[cig-commits] r16086 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Tue Dec 8 09:19:49 PST 2009
Author: becker
Date: 2009-12-08 09:19:49 -0800 (Tue, 08 Dec 2009)
New Revision: 16086
Modified:
mc/3D/CitcomS/trunk/lib/BC_util.c
mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
mc/3D/CitcomS/trunk/lib/Topo_gravity.c
Log:
Truly apply internal boundary conditions on horizontal surfaces for
E->mesh.toplayerbc > 0 (default, of course, 0)
Modified: mc/3D/CitcomS/trunk/lib/BC_util.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/BC_util.c 2009-12-08 16:01:03 UTC (rev 16085)
+++ mc/3D/CitcomS/trunk/lib/BC_util.c 2009-12-08 17:19:49 UTC (rev 16086)
@@ -129,7 +129,7 @@
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 */
+ if(level == E->mesh.levmax) /* NB */
BC[dirn][node] = value;
} /* end for loop i & j */
}
@@ -260,7 +260,7 @@
void assign_internal_bc(struct All_variables *E,int is_global)
{
- int lv, j, noz, k,lay,ncount;
+ int lv, j, noz, k,lay,ncount,ontop,onbottom;
/* stress or vel BC within a layer */
ncount = 0;
@@ -271,25 +271,31 @@
/* 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 */
+ ontop = ((k==noz) && (E->parallel.me_loc[3]==E->parallel.nprocz-1))?(TRUE):(FALSE);
+ onbottom = ((k==1) && (E->parallel.me_loc[3]==0))?(TRUE):(FALSE);
/* 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))))
+
+ if((!ontop)&&(!onbottom))
ncount++; /* not in top or bottom */
- if(E->mesh.topvbc != 1) { /* free slip top */
+ if(E->mesh.topvbc != 1) { /* free slip */
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);
+ if(ontop || onbottom)
+ 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);
+ if(ontop || onbottom)
+ 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 */
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);
+ if(ontop || onbottom)
+ 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);
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,1,0.0, SBX,0,lv,j);
+ if(ontop || onbottom)
+ 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);
}
}
}
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2009-12-08 16:01:03 UTC (rev 16085)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2009-12-08 17:19:49 UTC (rev 16086)
@@ -679,7 +679,7 @@
void ggrd_read_vtop_from_file(struct All_variables *E, int is_global)
{
MPI_Status mpi_stat;
- int mpi_rc,interpolate,timedep,use_codes,code,assign;
+ int mpi_rc,interpolate,timedep,use_codes,code,assign,ontop;
int mpi_inmsg, mpi_success_message = 1;
int m,el,i,k,i1,i2,ind,nodel,j,level, verbose;
int nox,noz,noy,noxl,noyl,nozl,lselect,idim,noxnoz,noxlnozl,save_codes,topnode,botnode;
@@ -1008,6 +1008,7 @@
*/
for(k = botnode;k <= topnode;k++){
+ ontop = ((k==nozl) && (E->parallel.me_loc[3]==E->parallel.nprocz-1))?(TRUE):(FALSE);
/* depth loop */
nodel = k + (j-1) * nozl + (i-1)*noxlnozl; /* top node = nozl + (j-1) * nozl + (i-1)*noxlnozl; */
if(fabs(v[2]) > cutoff){
@@ -1153,6 +1154,7 @@
*/
for(k = botnode;k <= topnode;k++){
+ ontop = ((k==noz) && (E->parallel.me_loc[3]==E->parallel.nprocz-1))?(TRUE):(FALSE);
nodel = k + (j-1) * noz + (i-1)*noxnoz ; /* node = k + (j-1) * nozg + (i-1)*noxgnozg; */
if(use_codes){
/* find code from v[1], theta */
@@ -1187,7 +1189,8 @@
E->sphere.cap[m].VB[1][nodel] = v[1]; /* theta */
E->sphere.cap[m].VB[2][nodel] = v[2]; /* phi */
}
- E->sphere.cap[m].VB[3][nodel] = 0.0; /* r */
+ if(use_vel && ontop)
+ E->sphere.cap[m].VB[3][nodel] = 0.0; /* r */
} /* end z */
} /* end x */
} /* end y */
Modified: mc/3D/CitcomS/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2009-12-08 16:01:03 UTC (rev 16085)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2009-12-08 17:19:49 UTC (rev 16086)
@@ -450,19 +450,27 @@
{
int m, i, j, k, n, d;
const unsigned sbc_flag[4] = {0, SBX, SBY, SBZ};
- const int stress_index[4][4] = { {0, 0, 0, 0},
- {0, 1, 4, 5}, /* N-S sides */
- {0, 4, 2, 6}, /* E-W sides */
- {0, 5, 6, 3} }; /* U-D sides */
+ /*
+ stress tensor is sorted like so: 1: xx 2: yy 3: zz 4: xy 5: xz 6: yz
+ tt pp rr tp tr pr
+ */
+ const int stress_index[4][4] = { {0, 0, 0, 0}, /* traction to stress tensor conversion */
+ {0, 1, 4, 5}, /* N-S sides, xx xy xz */
+ {0, 4, 2, 6}, /* E-W sides yx yy yz */
+ {0, 5, 6, 3} }; /* U-D sides zx zy zz */
- if(E->control.side_sbcs) {
+ int noxnoz;
+ noxnoz = E->lmesh.nox*E->lmesh.noz;
+
+ if(E->control.side_sbcs) { /* side boundary conditions */
+
for(m=1; m<=E->sphere.caps_per_proc; m++)
for(i=1; i<=E->lmesh.noy; i++)
for(j=1; j<=E->lmesh.nox; j++)
for(k=1; k<=E->lmesh.noz; k++) {
- n = k+(j-1)*E->lmesh.noz+(i-1)*E->lmesh.nox*E->lmesh.noz;
+ n = k+(j-1)*E->lmesh.noz+(i-1)*noxnoz;
for(d=1; d<=E->mesh.nsd; d++)
@@ -483,22 +491,43 @@
}
} else {
+ /*
+ regular case
- for(m=1; m<=E->sphere.caps_per_proc; m++)
- for(i=1; i<=E->lmesh.noy; i++)
- for(j=1; j<=E->lmesh.nox; j++)
- for(k=1; k<=E->lmesh.noz; k++) {
- n = k+(j-1)*E->lmesh.noz+(i-1)*E->lmesh.nox*E->lmesh.noz;
- for(d=1; d<=E->mesh.nsd; d++)
- if(E->node[m][n] & sbc_flag[d]) {
- if(i==1 || i==E->lmesh.noy)
- E->gstress[m][(n-1)*6+stress_index[d][2]] = E->sphere.cap[m].VB[d][n];
- if(j==1 || j==E->lmesh.nox)
- E->gstress[m][(n-1)*6+stress_index[d][1]] = E->sphere.cap[m].VB[d][n];
- if(k==1 || k==E->lmesh.noz)
- E->gstress[m][(n-1)*6+stress_index[d][3]] = E->sphere.cap[m].VB[d][n];
- }
- }
+ */
+ if(E->mesh.toplayerbc > 0){
+ /* internal BCs */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(i=1; i<=E->lmesh.noy; i++)
+ for(j=1; j<=E->lmesh.nox; j++)
+ for(k=1; k<=E->lmesh.noz; k++) {
+ n = k+(j-1)*E->lmesh.noz+(i-1)*noxnoz;
+ for(d=1; d<=E->mesh.nsd; d++)
+ if(E->node[m][n] & sbc_flag[d]) {
+ /* apply internal traction vector on horizontal surface */
+ if(layers(E,m,n) <= E->mesh.toplayerbc)
+ E->gstress[m][(n-1)*6+stress_index[d][3]] = E->sphere.cap[m].VB[d][n];
+ }
+ }
+
+ }else{
+ /* default */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(i=1; i<=E->lmesh.noy; i++)
+ for(j=1; j<=E->lmesh.nox; j++)
+ for(k=1; k<=E->lmesh.noz; k++) {
+ n = k+(j-1)*E->lmesh.noz+(i-1)*noxnoz;
+ for(d=1; d<=E->mesh.nsd; d++)
+ if(E->node[m][n] & sbc_flag[d]) {
+ if(i==1 || i==E->lmesh.noy)
+ E->gstress[m][(n-1)*6+stress_index[d][2]] = E->sphere.cap[m].VB[d][n];
+ if(j==1 || j==E->lmesh.nox)
+ E->gstress[m][(n-1)*6+stress_index[d][1]] = E->sphere.cap[m].VB[d][n];
+ if(k==1 || k==E->lmesh.noz)
+ E->gstress[m][(n-1)*6+stress_index[d][3]] = E->sphere.cap[m].VB[d][n];
+ }
+ }
+ }
}
}
More information about the CIG-COMMITS
mailing list