[cig-commits] r16404 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Thu Mar 11 16:41:47 PST 2010
Author: becker
Date: 2010-03-11 16:41:46 -0800 (Thu, 11 Mar 2010)
New Revision: 16404
Modified:
mc/3D/CitcomS/trunk/lib/BC_util.c
mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
Log:
Improvements for internal BCs.
Modified: mc/3D/CitcomS/trunk/lib/BC_util.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/BC_util.c 2010-03-10 19:55:02 UTC (rev 16403)
+++ mc/3D/CitcomS/trunk/lib/BC_util.c 2010-03-12 00:41:46 UTC (rev 16404)
@@ -281,7 +281,7 @@
/* node number is k, assuming no dependence on x and y */
if((lay = layers(E,j,k)) <= E->mesh.toplayerbc){
- if((!ontop)&&(!onbottom))
+ if((!ontop)&&(!onbottom)&&(lv==E->mesh.gridmax))
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);
@@ -318,18 +318,35 @@
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");
- ncount++; /* not in top or bottom */
+ if(lv == E->mesh.gridmax)
+ k = noz + E->mesh.toplayerbc;
+ else{
+ k = noz + (int)((float)E->mesh.toplayerbc / pow(2.,(float)(E->mesh.gridmax-lv)));
+ }
+ //fprintf(stderr,"BC_util: inner node: CPU: %i lv %i noz %i k %i\n",E->parallel.me,lv,noz,k);
+ if(k <= 1)
+ myerror(E,"out of bounds for noz and toplayerbc");
+ ontop = ((k==noz) && (E->parallel.me_loc[3]==E->parallel.nprocz-1))?(1):(0);
+ onbottom = ((k==1) && (E->parallel.me_loc[3]==0))?(1):(0);
+ if((!ontop)&&(!onbottom)&&(lv==E->mesh.gridmax))
+ 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(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);
+ 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);
+ 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);
+ 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 2010-03-10 19:55:02 UTC (rev 16403)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-03-12 00:41:46 UTC (rev 16404)
@@ -54,6 +54,8 @@
void construct_mat_group(struct All_variables *);
void temperatures_conform_bcs(struct All_variables *);
int layers(struct All_variables *,int ,int );
+void ggrd_vtop_helper_decide_on_internal_nodes(struct All_variables *, /* input */
+ int ,int ,int, int,int ,int *, int *,int *);
/*
@@ -686,7 +688,8 @@
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;
+ int nox,noz,noy,noxl,noyl,nozl,lselect,idim,noxnoz,
+ noxlnozl,save_codes,topnode,botnode;
char gmt_string[10],char_dummy;
static int lc =0; /* only for debugging */
double vin1[2],vin2[2],age,f1,f2,vscale,rout[3],cutoff,v[3],sin_theta,vx[4],
@@ -734,7 +737,7 @@
}
/*
-
+
read in plate code files? this will assign Euler vector
respective velocities
@@ -936,40 +939,8 @@
for (m=1;m <= E->sphere.caps_per_proc;m++) {
/* determine vertical nodes */
- if(allow_internal){
- /* internal */
- 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;
- }
- }else{
- /* only one internal node */
-
- if(level == E->mesh.gridmax){
- assign = TRUE;
- botnode = nozl + E->mesh.toplayerbc;
- topnode = botnode;
- }else{
- if(E->parallel.me == 0)
- fprintf(stderr,"WARNING: assigning single layer internal boundary condition only to top level multigrid\n");
- assign = FALSE;
- botnode = topnode = nozl;
- }
- }
- }else{ /* just top node */
- assign = TRUE;
- topnode = botnode = nozl;
- }
- if(verbose)
- fprintf(stderr,"ggrd_read_vtop_from_file: mixed: internal: %i assign: %i k: %i to %i (%i)\n",
- allow_internal,assign,botnode,topnode,nozl);
+ ggrd_vtop_helper_decide_on_internal_nodes(E,allow_internal,nozl,level,m,verbose,
+ &assign,&botnode,&topnode);
/*
loop through all horizontal nodes and assign boundary
conditions for all required levels
@@ -988,7 +959,7 @@
*/
if((is_global)&&(rout[1] > theta_max)){
if(!pole_warned){
- fprintf(stderr,"WARNING: shifting theta from %g (%g) to max theta %g (%g)\n",
+ fprintf(stderr,"ggrd_read_vtop_from_file: WARNING: shifting theta from %g (%g) to max theta %g (%g)\n",
rout[1],90-180/M_PI*rout[1],theta_max,90-180/M_PI*theta_max);
pole_warned = TRUE;
}
@@ -996,7 +967,7 @@
}
if((is_global)&&(rout[1] < theta_min)){
if(!pole_warned){
- fprintf(stderr,"WARNING: shifting theta from %g (%g) to min theta %g (%g)\n",
+ fprintf(stderr,"ggrd_read_vtop_from_file: WARNING: shifting theta from %g (%g) to min theta %g (%g)\n",
rout[1],90-180/M_PI*rout[1],theta_min,90-180/M_PI*theta_min);
pole_warned = TRUE;
}
@@ -1045,7 +1016,8 @@
E->NODE[level][m][nodel] = E->NODE[level][m][nodel] & (~SBX);
E->NODE[level][m][nodel] = E->NODE[level][m][nodel] | VBY;
E->NODE[level][m][nodel] = E->NODE[level][m][nodel] & (~SBY);
- }else{
+ }else{ fprintf(stderr,"t %i %i\n",level,nodel);
+
/* prescribed tractions */
E->NODE[level][m][nodel] = E->NODE[level][m][nodel] & (~VBX);
E->NODE[level][m][nodel] = E->NODE[level][m][nodel] | SBX;
@@ -1058,8 +1030,8 @@
} /* end y loop */
} /* actually assign */
} /* cap */
- } /* level */
- fprintf(stderr,"mixed_bc: %i free %i fixed for CPU %i\n",nfree,nfixed,E->parallel.me);
+ } /* MG level */
+ fprintf(stderr,"ggrd_read_vtop_from_file: mixed_bc: %i free %i fixed for CPU %i\n",nfree,nfixed,E->parallel.me);
} /* end mixed BC assign */
/*
@@ -1076,26 +1048,13 @@
if(save_codes) /* those will be surface nodes only */
gzprintf(fp1,"%3d %7d\n",m,E->lmesh.nsf);
}
-
+
+ /* top leevl */
+
for (m=1;m <= E->sphere.caps_per_proc;m++) {
- if(allow_internal){
- /* check for internal nodes in layers */
- for(k=noz;k >= 1;k--){
- if(layers(E,m,k) > E->mesh.toplayerbc) /* assumes regular mesh structure ! */
- break;
- }
- if(k == noz){
- assign = FALSE;
- }else{
- assign = TRUE;topnode = noz;botnode = k+1;
- }
- }else{ /* just top node */
- assign = TRUE;
- topnode = botnode = noz;
- }
- if(verbose)
- fprintf(stderr,"ggrd_read_vtop_from_file: internal: %i assign: %i k: %i to %i (%i)\n",
- allow_internal,assign,botnode,topnode,noz);
+ /* top level only */
+ ggrd_vtop_helper_decide_on_internal_nodes(E,allow_internal,E->lmesh.NOZ[E->mesh.gridmax],E->mesh.gridmax,m,verbose,
+ &assign,&botnode,&topnode);
if(assign){
for(i=1;i <= noy;i++) {/* loop through surface nodes */
for(j=1;j <= nox;j++) {
@@ -1204,7 +1163,7 @@
E->sphere.cap[m].VB[1][nodel] = 0; /* theta */
E->sphere.cap[m].VB[2][nodel] = 0; /* phi */
}else{
- /* regular no slip , assign velocities/tractuibs as BCs */
+ /* regular no slip , assign velocities/tractions as BCs */
E->sphere.cap[m].VB[1][nodel] = v[1]; /* theta */
E->sphere.cap[m].VB[2][nodel] = v[2]; /* phi */
}
@@ -1234,6 +1193,49 @@
}
+void ggrd_vtop_helper_decide_on_internal_nodes(struct All_variables *E, /* input */
+ int allow_internal,
+ int nozl,int level,int m,int verbose,
+ int *assign, /* output */
+ int *botnode,int *topnode)
+{
+ int k;
+ /* default is assign to top */
+ *assign = TRUE;
+ *botnode = *topnode = nozl;
+ /* determine vertical nodes */
+ if(allow_internal){
+ /* internal */
+ 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;
+ }
+ }else if(E->mesh.toplayerbc < 0){
+ /* only one internal node */
+ if(level == E->mesh.gridmax)
+ *botnode = nozl + E->mesh.toplayerbc;
+ else{
+ *botnode = nozl + (int)((float)E->mesh.toplayerbc / pow(2.,(float)(E->mesh.gridmax-level)));
+ }
+ *assign = TRUE;
+ *topnode = *botnode;
+ }else{
+ fprintf(stderr,"ggrd_vtop_helper_decide_on_internal_nodes: toplayerbc %i\n", E->mesh.toplayerbc);
+ myerror(E,"ggrd_vtop_helper_decide_on_internal_nodes: logic error");
+ }
+ }
+ if(verbose)
+ fprintf(stderr,"ggrd_vtop_helper_decide_on_internal_nodes: mixed: internal: %i assign: %i k: %i to %i (%i)\n",
+ allow_internal,*assign,*botnode,*topnode,nozl);
+
+}
void ggrd_read_age_from_file(struct All_variables *E, int is_global)
{
More information about the CIG-COMMITS
mailing list