[cig-commits] r16065 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Thu Dec 3 18:33:08 PST 2009
Author: becker
Date: 2009-12-03 18:33:07 -0800 (Thu, 03 Dec 2009)
New Revision: 16065
Modified:
mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c
mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c
Log:
Temporary fix for regular operations, for Eh's review.
Ggrd assignment still needs tending to.
Modified: mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c 2009-12-04 00:00:13 UTC (rev 16064)
+++ mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c 2009-12-04 02:33:07 UTC (rev 16065)
@@ -56,7 +56,7 @@
for(lv=E->mesh.gridmax;lv>=E->mesh.gridmin;lv--)
for (j=1;j<=E->sphere.caps_per_proc;j++) {
- noz = E->mesh.NOZ[lv];
+ noz = E->lmesh.NOZ[lv];
if(E->mesh.topvbc != 1) { /* free slip top */
horizontal_bc(E,E->sphere.cap[j].VB,noz,1,0.0,VBX,0,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,noz,3,0.0,VBZ,1,lv,j);
@@ -143,7 +143,7 @@
lev = E->mesh.levmax;
for (j=1;j<=E->sphere.caps_per_proc;j++) {
- noz = E->mesh.noz;
+ noz = E->lmesh.noz;
if(E->mesh.toptbc == 1) {
horizontal_bc(E,E->sphere.cap[j].TB,noz,3,E->control.TBCtopval,TBZ,1,lev,j);
horizontal_bc(E,E->sphere.cap[j].TB,noz,3,E->control.TBCtopval,FBZ,0,lev,j);
@@ -186,14 +186,14 @@
/*
why don't we have only one horizontal_bc for both regional and full
-versions ? TWB
+versions ?
*/
-static void horizontal_bc(E,BC,ROW,dirn,value,mask,onoff,level,m)
+static void horizontal_bc(E,BC,row,dirn,value,mask,onoff,level,m)
struct All_variables *E;
float *BC[];
- int ROW;
+ int row;
int dirn;
float value;
unsigned int mask;
@@ -201,7 +201,7 @@
int level,m;
{
- int i,j,node,rowl;
+ int i,j,node;
static short int warned = FALSE;
/* safety feature */
if(dirn > E->mesh.nsd)
@@ -209,31 +209,22 @@
/*
- I commented this out since the regular calls are always either
- ROW == 1 or ROW = E->lmesh.NOZ[level], and I want to allow
- assigning internal BCs TWB
-
- */
- /* if (ROW==1) */
- /* rowl = 1; */
- /* else */
- /* rowl = E->lmesh.NOZ[level]; */
- rowl = ROW;
- /*
now warn if not assigning to top or bottom, but allow
*/
- if((ROW != 1) && (ROW != E->lmesh.NOZ[level])){
+ if((row != 1) && (row != E->lmesh.NOZ[level])){
if(!warned){
- fprintf(stderr,"horizontal_bc: CPU %i: assigning internal BC\n",
- E->parallel.me);
+ fprintf(stderr,"horizontal_bc: CPU %i: assigning internal BC, %i %i %i\n",
+ E->parallel.me,row, E->lmesh.NOZ[level], E->mesh.noz);
warned = TRUE;
}
+ if(row > E->lmesh.NOZ[level])
+ myerror(E,"horizontal_bc: error, 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 = rowl+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
+ node = row+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
E->NODE[level][m][node] = E->NODE[level][m][node] & (~ mask);
} /* end for loop i & j */
}
@@ -242,7 +233,7 @@
else {
for(j=1;j<=E->lmesh.NOY[level];j++)
for(i=1;i<=E->lmesh.NOX[level];i++) {
- node = rowl+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
+ node = row+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
E->NODE[level][m][node] = E->NODE[level][m][node] | (mask);
if(level==E->mesh.levmax) /* NB */
BC[dirn][node] = value;
@@ -251,14 +242,14 @@
}else{
- if ( ( (ROW==1) && (E->parallel.me_loc[3]==0) ) ||
- ( (ROW==E->mesh.NOZ[level]) && (E->parallel.me_loc[3]==E->parallel.nprocz-1) ) ) {
+ 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 = rowl+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
+ node = row+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
E->NODE[level][m][node] = E->NODE[level][m][node] & (~ mask);
} /* end for loop i & j */
}
@@ -267,14 +258,14 @@
else {
for(j=1;j<=E->lmesh.NOY[level];j++)
for(i=1;i<=E->lmesh.NOX[level];i++) {
- node = rowl+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
+ node = row+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
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;
}
@@ -317,7 +308,7 @@
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->mesh.NOZ[lv];
+ 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){
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2009-12-04 00:00:13 UTC (rev 16064)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2009-12-04 02:33:07 UTC (rev 16065)
@@ -700,6 +700,15 @@
#else
myerror(E,"ggrd_read_vtop_from_file needs to use GZDIR (set USE_GZDIR flag) because of code output");
#endif
+
+ /* top level number of nodes */
+ noxg = E->lmesh.nox;
+ nozg = E->lmesh.noz;
+ noyg = E->lmesh.noy;
+ noxgnozg = noxg*nozg;
+
+
+
/* top processor check */
top_echo = E->parallel.nprocz-1;
/* output of warning messages */
@@ -732,12 +741,7 @@
fprintf(stderr,"ggrd_read_vtop_from_file: init stage, assigning %s, mixed mode: %i\n",
((E->mesh.topvbc)?("velocities"):("tractions")),E->control.ggrd_allow_mixed_vbcs);
- /* global, top level number of nodes */
- noxg = E->lmesh.nox;
- nozg = E->lmesh.noz;
- noyg = E->lmesh.noy;
- noxgnozg = noxg*nozg;
-
+
if(use_vel){
/*
velocity scaling, assuming input is cm/yr
@@ -918,10 +922,12 @@
slip */
cutoff = E->control.ggrd.svp->fmaxlim[0] + 1e-5;
for(level=E->mesh.gridmax;level >= E->mesh.gridmin;level--){/* multigrid levels */
+ /* all levels */
noxl = E->lmesh.NOX[level];
noyl = E->lmesh.NOY[level];
nozl = E->lmesh.NOZ[level];
noxlnozl = noxl*nozl;
+
for (m=1;m <= E->sphere.caps_per_proc;m++) {
/* determine vertical nodes */
if(is_internal){
@@ -951,9 +957,9 @@
loop through all horizontal nodes
*/
if(assign){
- for(i=1;i<=noyl;i++){
- for(j=1;j<=noxl;j++) {
- nodel = j * nozl + (i-1)*noxlnozl; /* top node = nozl + (j-1) * nozl + (i-1)*noxlnozl; */
+ for(i=1;i <= noyl;i++){
+ for(j=1;j <= noxl;j++) {
+ nodel = nozl + (j-1) * nozl + (i-1)*noxlnozl; /* top node = nozl + (j-1) * nozl + (i-1)*noxlnozl; */
/* node location */
rout[1] = E->SX[level][m][1][nodel]; /* theta,phi */
rout[2] = E->SX[level][m][2][nodel];
@@ -1078,7 +1084,7 @@
if(assign){
for(i=1;i <= noyg;i++) {/* loop through surface nodes */
for(j=1;j <= noxg;j++) {
- nodel = (i-1)*noxgnozg + j * nozg; /* top node = nozg + (j-1) * nozg + (i-1)*noxgnozg; */
+ nodel = nozg + (j-1) * nozg + (i-1)*noxgnozg; /* top node = nozg + (j-1) * nozg + (i-1)*noxgnozg; */
/* */
rout[1] = E->sx[m][1][nodel]; /* theta,phi coordinates */
rout[2] = E->sx[m][2][nodel];
@@ -1152,7 +1158,7 @@
*/
for(k = botnode;k <= topnode;k++){
- nodel = k + (i-1)*noxgnozg + (j-1) * nozg; /* node = k + (j-1) * nozg + (i-1)*noxgnozg; */
+ nodel = k + (j-1) * nozg + (i-1)*noxgnozg ; /* node = k + (j-1) * nozg + (i-1)*noxgnozg; */
if(use_codes){
/* find code from v[1], theta */
code = (int)(v[1] + 0.5);
Modified: mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c 2009-12-04 00:00:13 UTC (rev 16064)
+++ mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c 2009-12-04 02:33:07 UTC (rev 16065)
@@ -393,10 +393,10 @@
/* ========================================================= */
-static void horizontal_bc(E,BC,ROW,dirn,value,mask,onoff,level,m)
+static void horizontal_bc(E,BC,row,dirn,value,mask,onoff,level,m)
struct All_variables *E;
float *BC[];
- int ROW;
+ int row;
int dirn;
float value;
unsigned int mask;
@@ -404,38 +404,27 @@
int level,m;
{
- int i,j,node,rowl;
+ int i,j,node;
static short int warned = FALSE;
/* safety feature */
if(dirn > E->mesh.nsd)
return;
- /*
- I commented this out since the regular calls are always either
- ROW == 1 or ROW = E->lmesh.NOZ[level], and I want to allow
- assigning internal BCs TWB
-
- */
- /* if (ROW==1) */
- /* rowl = 1; */
- /* else */
- /* rowl = E->lmesh.NOZ[level]; */
- rowl = ROW;
- /*
- now warn if not assigning to top or bottom, but allow
- */
- if((ROW != 1) && (ROW != E->lmesh.NOZ[level])){
+ if((row != 1) && (row != E->lmesh.NOZ[level])){
+ /* not in top or bottom row of this processor */
if(!warned){
- fprintf(stderr,"horizontal_bc: CPU %i: assigning internal BC\n",
- E->parallel.me);
+ fprintf(stderr,"horizontal_bc: CPU %i: assigning internal BC, %i %i %i\n",
+ E->parallel.me,row,E->lmesh.NOZ[level],E->mesh.noz);
warned = TRUE;
}
-
+ if(row > E->lmesh.NOZ[level])
+ myerror(E,"horizontal_bc: error, 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 = rowl+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
+ node = row+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
E->NODE[level][m][node] = E->NODE[level][m][node] & (~ mask);
} /* end for loop i & j */
}
@@ -444,7 +433,7 @@
else {
for(j=1;j<=E->lmesh.NOY[level];j++)
for(i=1;i<=E->lmesh.NOX[level];i++) {
- node = rowl+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
+ node = row+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
E->NODE[level][m][node] = E->NODE[level][m][node] | (mask);
if(level==E->mesh.levmax) /* NB */
@@ -452,14 +441,14 @@
} /* end for loop i & j */
}
}else{
- if ( (ROW==1 && E->parallel.me_loc[3]==0) ||
- (ROW==E->lmesh.NOZ[level] && E->parallel.me_loc[3]==E->parallel.nprocz-1) ) {
+ if ( (row==1 && E->parallel.me_loc[3]==0) ||
+ (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 = rowl+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
+ node = row+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
E->NODE[level][m][node] = E->NODE[level][m][node] & (~ mask);
} /* end for loop i & j */
}
@@ -468,7 +457,7 @@
else {
for(j=1;j<=E->lmesh.NOY[level];j++)
for(i=1;i<=E->lmesh.NOX[level];i++) {
- node = rowl+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
+ node = row+(i-1)*E->lmesh.NOZ[level]+(j-1)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
E->NODE[level][m][node] = E->NODE[level][m][node] | (mask);
if(level==E->mesh.levmax) /* NB */
@@ -476,7 +465,7 @@
} /* end for loop i & j */
}
- } /* end for if ROW */
+ } /* end for if row */
}
return;
@@ -526,7 +515,7 @@
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->mesh.NOZ[lv];
+ 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){
More information about the CIG-COMMITS
mailing list