[cig-commits] r16073 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Fri Dec 4 15:16:49 PST 2009
Author: becker
Date: 2009-12-04 15:16:49 -0800 (Fri, 04 Dec 2009)
New Revision: 16073
Modified:
mc/3D/CitcomS/trunk/lib/BC_util.c
Log:
Modified handling of internal BCs, as per Eh's suggestion, now handled
in a separate function.
Modified: mc/3D/CitcomS/trunk/lib/BC_util.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/BC_util.c 2009-12-04 22:48:42 UTC (rev 16072)
+++ mc/3D/CitcomS/trunk/lib/BC_util.c 2009-12-04 23:16:49 UTC (rev 16073)
@@ -27,20 +27,15 @@
*/
#include "global_defs.h"
+void horizontal_bc(struct All_variables *,float *[],int,int,float,unsigned int,char,int,int);
+void internal_horizontal_bc(struct All_variables *,float *[],int,int,float,unsigned int,char,int,int);
/*
-assign boundary conditions to a horizontal layer of nodes
+assign boundary conditions to a horizontal layer of nodes on top and
+bottom, only applying to those processors
-row > 0: regular operation, only assign BCs to bottom (row = 1) and
- top (row = E->lmesh.NOZ[level]) processors, as in the
- previous version where horizontal_bc was a function in both
- Full and Regional BC subroutines
-
-row < 0: assign BCs to -row, no matter what processor, to allow
- setting internal boundary conditions
-
- */
+*/
void horizontal_bc(E,BC,row,dirn,value,mask,onoff,level,m)
struct All_variables *E;
float *BC[];
@@ -53,26 +48,18 @@
{
int i,j,node,noxnoz;
- static short int warned = FALSE;
/* safety feature */
if(dirn > E->mesh.nsd)
return;
noxnoz = E->lmesh.NOX[level]*E->lmesh.NOZ[level];
- if(row < 0){
- /*
- assignment to any row
- */
- row = -row;
- if((row != E->lmesh.NOZ[level]) && (row != 1) && (!warned)){
- fprintf(stderr,"horizontal_bc: CPU %4i: assigning internal BC, row: %4i nozl: %4i noz: %4i\n",
- E->parallel.me,row, E->lmesh.NOZ[level], E->mesh.noz);
- warned = TRUE;
- }
- if((row > E->lmesh.NOZ[level])||(row<1))
- myerror(E,"horizontal_bc: error, out of bounds");
-
+ /* regular operation, assign only if in top
+ (row==E->lmesh.NOZ[level]) or bottom (row==1) processor */
+
+ if ( ( (row==1) && (E->parallel.me_loc[3]==0) ) || /* bottom or top */
+ ( (row==E->lmesh.NOZ[level]) && (E->parallel.me_loc[3]==E->parallel.nprocz-1) ) ) {
+
/* turn bc marker to zero */
if (onoff == 0) {
for(j=1;j<=E->lmesh.NOY[level];j++)
@@ -92,40 +79,66 @@
BC[dirn][node] = value;
} /* end for loop i & j */
}
- }else{
- /* regular operation, assign only if in top
- (row==E->lmesh.NOZ[level]) or bottom (row==1) processor */
- if ( ( (row==1) && (E->parallel.me_loc[3]==0) ) || /* bottom or top */
- ( (row==E->lmesh.NOZ[level]) && (E->parallel.me_loc[3]==E->parallel.nprocz-1) ) ) {
-
- /* turn bc marker to zero */
- if (onoff == 0) {
- for(j=1;j<=E->lmesh.NOY[level];j++)
- for(i=1;i<=E->lmesh.NOX[level];i++) {
- node = row+(i-1)*E->lmesh.NOZ[level]+(j-1)*noxnoz;
- E->NODE[level][m][node] = E->NODE[level][m][node] & (~ mask);
- } /* end for loop i & j */
- }
-
- /* turn bc marker to one */
- else {
- for(j=1;j<=E->lmesh.NOY[level];j++)
- for(i=1;i<=E->lmesh.NOX[level];i++) {
- node = row+(i-1)*E->lmesh.NOZ[level]+(j-1)*noxnoz;
- E->NODE[level][m][node] = E->NODE[level][m][node] | (mask);
- if(level==E->mesh.levmax) /* NB */
- BC[dirn][node] = value;
- } /* end for loop i & j */
- }
-
- } /* end for if row */
+ } /* end for if row */
+ return;}
+
+
+
+/*
+
+assign boundary conditions to a horizontal layer of nodes within mesh,
+without consideration of being in top or bottom processor
+
+*/
+void internal_horizontal_bc(E,BC,row,dirn,value,mask,onoff,level,m)
+ struct All_variables *E;
+ float *BC[];
+ int row;
+ int dirn;
+ float value;
+ unsigned int mask;
+ char onoff;
+ int level,m;
+
+{
+ int i,j,node,noxnoz;
+ /* safety feature */
+ if(dirn > E->mesh.nsd)
+ return;
+
+ noxnoz = E->lmesh.NOX[level]*E->lmesh.NOZ[level];
+
+ /*
+ assignment to any row, any processor
+ */
+
+ if((row > E->lmesh.NOZ[level])||(row < 1))
+ myerror(E,"internal_horizontal_bc: error, row out of bounds");
+
+ /* turn bc marker to zero */
+ if (onoff == 0) {
+ for(j=1;j<=E->lmesh.NOY[level];j++)
+ for(i=1;i<=E->lmesh.NOX[level];i++) {
+ node = row+(i-1)*E->lmesh.NOZ[level]+(j-1)*noxnoz;
+ E->NODE[level][m][node] = E->NODE[level][m][node] & (~ mask);
+ } /* end for loop i & j */
+ }else {
+ /* turn bc marker to one */
+ for(j=1;j<=E->lmesh.NOY[level];j++)
+ for(i=1;i<=E->lmesh.NOX[level];i++) {
+ node = row+(i-1)*E->lmesh.NOZ[level]+(j-1)*noxnoz;
+ E->NODE[level][m][node] = E->NODE[level][m][node] | (mask);
+ if(level==E->mesh.levmax) /* NB */
+ BC[dirn][node] = value;
+ } /* end for loop i & j */
}
+
+
return;
}
-
void strip_bcs_from_residual(E,Res,level)
struct All_variables *E;
double **Res;
@@ -247,31 +260,36 @@
void assign_internal_bc(struct All_variables *E,int is_global)
{
- int lv, j, noz, k,node,lay;
+ int lv, j, noz, k,lay,ncount;
/* stress or vel BC within a layer */
-
+ ncount = 0;
if(E->mesh.toplayerbc > 0){
for(lv=E->mesh.gridmax;lv>=E->mesh.gridmin;lv--)
for (j=1;j<=E->sphere.caps_per_proc;j++) {
noz = E->lmesh.NOZ[lv];
- for(k=noz-1;k >= 1;k--){ /* assumes regular grid */
- node = k; /* global node number */
- if((lay = layers(E,j,node)) <= E->mesh.toplayerbc){
+ /* we're looping through all nodes for the possibility that
+ there are several internal processors which need BCs */
+ for(k=noz;k >= 1;k--){ /* assumes regular grid */
+ /* node number is k, assuming no dependence on x and y */
+ if((lay = layers(E,j,k)) <= E->mesh.toplayerbc){
+ if((!((k==1) && (E->parallel.me_loc[3]==0))) &&
+ (!((k==noz) && (E->parallel.me_loc[3]==E->parallel.nprocz-1))))
+ ncount++; /* not in top or bottom */
if(E->mesh.topvbc != 1) { /* free slip top */
- horizontal_bc(E,E->sphere.cap[j].VB,-k,1,0.0,VBX,0,lv,j);
- horizontal_bc(E,E->sphere.cap[j].VB,-k,3,0.0,VBZ,1,lv,j);
- horizontal_bc(E,E->sphere.cap[j].VB,-k,2,0.0,VBY,0,lv,j);
- horizontal_bc(E,E->sphere.cap[j].VB,-k,1,E->control.VBXtopval,SBX,1,lv,j);
- horizontal_bc(E,E->sphere.cap[j].VB,-k,3,0.0,SBZ,0,lv,j);
- horizontal_bc(E,E->sphere.cap[j].VB,-k,2,E->control.VBYtopval,SBY,1,lv,j);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,1,0.0,VBX,0,lv,j);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,3,0.0,VBZ,1,lv,j);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,2,0.0,VBY,0,lv,j);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,1,E->control.VBXtopval,SBX,1,lv,j);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,3,0.0,SBZ,0,lv,j);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,2,E->control.VBYtopval,SBY,1,lv,j);
}else{ /* no slip */
- horizontal_bc(E,E->sphere.cap[j].VB,-k,1,E->control.VBXtopval,VBX,1,lv,j);
- horizontal_bc(E,E->sphere.cap[j].VB,-k,3,0.0,VBZ,1,lv,j);
- horizontal_bc(E,E->sphere.cap[j].VB,-k,2,E->control.VBYtopval,VBY,1,lv,j);
- horizontal_bc(E,E->sphere.cap[j].VB,-k,1,0.0,SBX,0,lv,j);
- horizontal_bc(E,E->sphere.cap[j].VB,-k,3,0.0,SBZ,0,lv,j);
- horizontal_bc(E,E->sphere.cap[j].VB,-k,2,0.0,SBY,0,lv,j);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,1,E->control.VBXtopval,VBX,1,lv,j);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,3,0.0,VBZ,1,lv,j);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,2,E->control.VBYtopval,VBY,1,lv,j);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,1,0.0,SBX,0,lv,j);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,3,0.0,SBZ,0,lv,j);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,2,0.0,SBY,0,lv,j);
}
}
}
@@ -282,6 +300,12 @@
ggrd_read_vtop_from_file(E, is_global, 1);
#endif
} /* end toplayerbc > 0 branch */
+
+ if(ncount)
+ fprintf(stderr,"assign_internal_bc: CPU %4i (%s): WARNING: assigned internal %s BCs to %6i nodes\n",
+ E->parallel.me,((E->parallel.me_loc[3]==0)&&(E->parallel.nprocz!=1))?("bottom"):
+ ((E->parallel.me_loc[3]==E->parallel.nprocz-1)?("top"):("interior")),
+ (E->mesh.topvbc!=1)?("stress"):("velocity"),ncount);
}
More information about the CIG-COMMITS
mailing list