[cig-commits] r16402 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Wed Mar 10 11:38:15 PST 2010
Author: becker
Date: 2010-03-10 11:38:15 -0800 (Wed, 10 Mar 2010)
New Revision: 16402
Modified:
mc/3D/CitcomS/trunk/lib/BC_util.c
mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c
mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c
mc/3D/CitcomS/trunk/lib/Topo_gravity.c
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Experimental implementation of a single node layer for internal
boundary condition assignment, assigned if toplayerbc < 0, to nodes at
layer noz+toplayerbc, but only for the top level multigrid, for now.
Modified: mc/3D/CitcomS/trunk/lib/BC_util.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/BC_util.c 2010-03-10 15:50:23 UTC (rev 16401)
+++ mc/3D/CitcomS/trunk/lib/BC_util.c 2010-03-10 19:38:15 UTC (rev 16402)
@@ -255,10 +255,13 @@
options:
-toplayerbc > 0: assign surface BC down to toplayerbc nd
+toplayerbc > 0: assign surface boundary condition down to all nodes within the the toplayerbc
+ layer, e.g. layer one for lithosphere
toplayerbc == 0: no action
+toplayerbc < 0: assign surface boundary condition within medium at node -toplayerbc depth, ie.
+ toplayerbc = -1 is one node underneath surface
- */
+*/
void assign_internal_bc(struct All_variables *E,int is_global)
{
@@ -307,8 +310,46 @@
if(E->control.ggrd.vtop_control)
ggrd_read_vtop_from_file(E, is_global, 1);
#endif
- } /* end toplayerbc > 0 branch */
-
+ /* end toplayerbc > 0 branch */
+ }else if(E->mesh.toplayerbc < 0){
+ /* internal node at noz-toplayerbc */
+ 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];
+ /* we're looping through all nodes for the possibility that
+ there are several internal processors which need BCs */
+ k = noz + E->mesh.toplayerbc;
+ if(k < 1)myerror(E,"out of bounds for noz and toplayerbc");
+ onbottom = ((k==1) && (E->parallel.me_loc[3]==0))?(1):(0);
+ if(!onbottom)
+ ncount++; /* not in top or bottom */
+ if(E->mesh.topvbc != 1) { /* free slip */
+ internal_horizontal_bc(E,E->sphere.cap[j].VB,k,1,0.0,VBX,0,lv,j);
+ if(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);
+ if(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);
+ if(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);
+ if(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);
+ }
+ }
+ /* read in velocities/stresses from grd file? */
+#ifdef USE_GGRD
+ if(E->control.ggrd.vtop_control)
+ 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"):
Modified: mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c 2010-03-10 15:50:23 UTC (rev 16401)
+++ mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c 2010-03-10 19:38:15 UTC (rev 16402)
@@ -120,7 +120,7 @@
/*
apply stress or velocity boundary conditions, read from file
settings are to be implemented in those routines (will only do
- anything at present, if E->mesh.toplayerbc > 0
+ anything at present, if E->mesh.toplayerbc != 0
*/
assign_internal_bc(E,1);
#ifdef USE_GGRD
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-03-10 15:50:23 UTC (rev 16401)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-03-10 19:38:15 UTC (rev 16402)
@@ -326,7 +326,8 @@
/*
-read in material, i.e. viscosity prefactor from ggrd file, this will get assigned if
+read in material, i.e. viscosity prefactor from ggrd file, this will
+get assigned for all nodes if their
layer <= E->control.ggrd.mat_control for E->control.ggrd.mat_control > 0
@@ -335,6 +336,9 @@
layer == -E->control.ggrd.mat_control for E->control.ggrd.mat_control < 0
+the grd model can be 2D (a layer in itself), or 3D (a model with
+several layers)
+
*/
void ggrd_read_mat_from_file(struct All_variables *E, int is_global)
{
@@ -679,7 +683,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,ontop;
+ int mpi_rc,interpolate,timedep,use_codes,code,assign,ontop,kinc;
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;
@@ -704,7 +708,7 @@
noy = E->lmesh.noy;
noxnoz = nox*noz;
- if(E->mesh.toplayerbc > 0)
+ if(E->mesh.toplayerbc != 0)
allow_internal = TRUE;
else
allow_internal = FALSE;
@@ -934,19 +938,35 @@
/* determine vertical nodes */
if(allow_internal){
/* internal */
- /* check for internal nodes in layers */
- for(k=nozl;k >= 1;k--){
- if(layers(E,m,k) > E->mesh.toplayerbc) /* assume regular mesh structure */
- break;
- }
- if(k == nozl){ /* */
- assign = FALSE;
+ if(E->mesh.toplayerbc > 0){
+ /* check for internal nodes in layers */
+ for(k=nozl;k >= 1;k--){
+ if(layers(E,m,k) > E->mesh.toplayerbc) /* assume regular mesh structure */
+ break;
+ }
+ if(k == nozl){ /* */
+ assign = FALSE;
+ }else{
+ assign = TRUE;topnode = nozl;botnode = k+1;
+ }
+ kinc = 1;
}else{
- assign = TRUE;topnode = nozl;botnode = k+1;
+ /* only one internal node */
+ assign = TRUE;
+ topnode = nozl;
+ if(level == E->mesh.gridmax){
+ botnode = nozl + E->mesh.toplayerbc;
+ kinc = topnode - botnode;
+ }else{
+ if(E->parallel.me == 0)
+ fprintf(stderr,"WARNING: assigning single layer internal boundary condition only to top level multigrid\n");
+ botnode = nozl;kinc = 1;
+ }
}
}else{ /* just top node */
assign = TRUE;
topnode = botnode = nozl;
+ kinc = 1;
}
if(verbose)
fprintf(stderr,"ggrd_read_vtop_from_file: mixed: internal: %i assign: %i k: %i to %i (%i)\n",
@@ -1007,7 +1027,7 @@
*/
- for(k = botnode;k <= topnode;k++){
+ for(k = botnode;k <= topnode;k += kinc){
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; */
@@ -1153,7 +1173,7 @@
XXX
*/
- for(k = botnode;k <= topnode;k++){
+ for(k = botnode;k <= topnode;k += kinc){
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){
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2010-03-10 15:50:23 UTC (rev 16401)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2010-03-10 19:38:15 UTC (rev 16402)
@@ -434,8 +434,17 @@
input_int("topvbc",&(E->mesh.topvbc),"0",m);
input_int("botvbc",&(E->mesh.botvbc),"0",m);
- input_int("toplayerbc",&(E->mesh.toplayerbc),"0",m); /* apply surface BC throughout all layer nodes */
+ input_int("toplayerbc",&(E->mesh.toplayerbc),"0",m); /* > 0: apply
+ surface BC
+ throughout
+ all layer
+ nodes of layers toplayerbc (i.e. = 1 for top layer)
+
+ < 0: apply to single node layer noz+toplayerbc
+
+ */
+
input_float("topvbxval",&(E->control.VBXtopval),"0.0",m);
input_float("botvbxval",&(E->control.VBXbotval),"0.0",m);
input_float("topvbyval",&(E->control.VBYtopval),"0.0",m);
@@ -511,6 +520,9 @@
30 60
+
+ in the example above, the input grid is a layer. if it's a 3D model, provide
+ ggrd_mat_depth_file, akin to temperature input
*/
Modified: mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c 2010-03-10 15:50:23 UTC (rev 16401)
+++ mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c 2010-03-10 19:38:15 UTC (rev 16402)
@@ -119,7 +119,7 @@
apply stress or velocity boundary conditions, read from file
settings are to be implemented in those routines (will only do
- anything at present, if E->mesh.toplayerbc > 0
+ anything at present, if E->mesh.toplayerbc != 0
*/
assign_internal_bc(E,0);
Modified: mc/3D/CitcomS/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2010-03-10 15:50:23 UTC (rev 16401)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2010-03-10 19:38:15 UTC (rev 16402)
@@ -492,11 +492,11 @@
} else {
/*
- regular case
+ no side boundary conditions
*/
if(E->mesh.toplayerbc > 0){
- /* internal BCs */
+ /* internal BCs for top toplayerbc layers */
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++)
@@ -510,6 +510,19 @@
}
}
+ }else if(E->mesh.toplayerbc < 0){
+ /* internal BCs for a single node layer noz+toplayerbc down */
+ 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++){
+ k = E->lmesh.noz + E->mesh.toplayerbc;
+ 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 */
+ 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++)
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2010-03-10 15:50:23 UTC (rev 16401)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2010-03-10 19:38:15 UTC (rev 16402)
@@ -346,7 +346,7 @@
int botvbc;
int sidevbc;
- int toplayerbc; /* apply surface BC throughout layer */
+ int toplayerbc; /* apply surface BC throughout top layers, or for a single internal node */
int periodic_x;
int periodic_y;
More information about the CIG-COMMITS
mailing list