[cig-commits] r16067 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Fri Dec 4 09:44:24 PST 2009
Author: becker
Date: 2009-12-04 09:44:24 -0800 (Fri, 04 Dec 2009)
New Revision: 16067
Modified:
mc/3D/CitcomS/trunk/lib/BC_util.c
mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c
mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c
mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c
mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
mc/3D/CitcomS/trunk/lib/ggrd_handling.h
Log:
Suggested unified treatment of horizontal_bc (moved to BC_util) for
review by Eh
Improved handling of grd-read velocity and traction boundary conditions.
Modified: mc/3D/CitcomS/trunk/lib/BC_util.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/BC_util.c 2009-12-04 04:42:31 UTC (rev 16066)
+++ mc/3D/CitcomS/trunk/lib/BC_util.c 2009-12-04 17:44:24 UTC (rev 16067)
@@ -28,8 +28,104 @@
#include "global_defs.h"
+/*
+assign boundary conditions to a horizontal layer of nodes
+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[];
+ int row;
+ int dirn;
+ float value;
+ unsigned int mask;
+ char onoff;
+ int level,m;
+
+{
+ 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");
+
+ /* 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 */
+ }
+ }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 */
+ }
+ return;
+}
+
+
+
void strip_bcs_from_residual(E,Res,level)
struct All_variables *E;
double **Res;
@@ -137,7 +233,58 @@
return;
}
+/*
+facility to apply internal velocity or stress conditions after
+top/bottom
+options:
+
+toplayerbc > 0: assign surface BC down to toplayerbc nd
+toplayerbc == 0: no action
+
+ */
+void assign_internal_bc(struct All_variables *E,int is_global)
+{
+
+ int lv, j, noz, k,node,lay;
+ /* stress or vel BC within a layer */
+
+
+ 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){
+ 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);
+ }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);
+ }
+ }
+ }
+ }
+ /* 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 */
+}
+
+
+
/* End of file */
Modified: mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c 2009-12-04 04:42:31 UTC (rev 16066)
+++ mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c 2009-12-04 17:44:24 UTC (rev 16067)
@@ -35,7 +35,8 @@
#endif
/* ========================================== */
-static void horizontal_bc();
+void horizontal_bc(struct All_variables *,float *[],int,int,float,unsigned int,char,int,int);
+void assign_internal_bc(struct All_variables *,int );
static void velocity_apply_periodic_bcs();
static void temperature_apply_periodic_bcs();
void read_temperature_boundary_from_file(struct All_variables *);
@@ -48,7 +49,6 @@
{
void velocity_imp_vert_bc();
void velocity_apply_periodicapply_periodic_bcs();
- void assign_internal_bc(struct All_variables * );
void apply_side_sbc();
@@ -64,11 +64,6 @@
horizontal_bc(E,E->sphere.cap[j].VB,noz,1,E->control.VBXtopval,SBX,1,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,noz,3,0.0,SBZ,0,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,noz,2,E->control.VBYtopval,SBY,1,lv,j);
-#ifdef USE_GGRD
- /* Ggrd traction control */
- if((lv==E->mesh.gridmax) && E->control.ggrd.vtop_control)
- ggrd_read_vtop_from_file(E, 1,0);
-#endif
}
if(E->mesh.botvbc != 1) { /* free slip bottom */
horizontal_bc(E,E->sphere.cap[j].VB,1,1,0.0,VBX,0,lv,j);
@@ -121,15 +116,17 @@
/* If any imposed internal velocity structure it goes here */
-
+
/*
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
*/
- assign_internal_bc(E);
-
-
+ assign_internal_bc(E,1);
+#ifdef USE_GGRD
+ if(E->control.ggrd.vtop_control) /* assign stress or velocity BCs */
+ ggrd_read_vtop_from_file(E,1);
+#endif
return; }
/* ========================================== */
@@ -183,94 +180,7 @@
/* ========================================================= */
-/*
-why don't we have only one horizontal_bc for both regional and full
-versions ?
-
-
- */
-static void 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;
- static short int warned = FALSE;
- /* safety feature */
- if(dirn > E->mesh.nsd)
- return;
-
-
- /*
- now warn if not assigning to top or bottom, but allow
- */
- if((row != 1) && (row != E->lmesh.NOZ[level])){
- if(!warned){
- 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 = 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 */
- }
-
- /* 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)*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 */
- }
- }else{
-
-
- 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)*E->lmesh.NOX[level]*E->lmesh.NOZ[level];
- 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)*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 */
- }
- return;
-}
-
-
static void velocity_apply_periodic_bcs(E)
struct All_variables *E;
{
@@ -288,58 +198,9 @@
}
-/*
-facility to apply internal velocity or stress conditions after
-top/bottom
-options:
-toplayerbc > 0: assign surface BC down to toplayerbc nd
-
- */
-void assign_internal_bc(struct All_variables *E)
-{
-
- int lv, j, noz, k,node,lay;
- /* stress or vel BC within a layer */
-
-
- 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){
- 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);
- }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);
- }
- }
- }
- }
-#ifdef USE_GGRD
- if(E->control.ggrd.vtop_control)
- ggrd_read_vtop_from_file(E, 1, 1);
-#endif
- }
-}
-
-
-
-
/* version */
/* $Id$ */
Modified: mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c 2009-12-04 04:42:31 UTC (rev 16066)
+++ mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c 2009-12-04 17:44:24 UTC (rev 16067)
@@ -222,9 +222,7 @@
case 1: /* velocity boundary conditions */
#ifdef USE_GGRD
- if(E->control.ggrd.vtop_control){
- ggrd_read_vtop_from_file(E, 1,0);
- }else{
+ if(!E->control.ggrd.vtop_control){ /* grd control is called from boundary conditions subroutine */
#endif
nnn=nox*noy;
for(i=1;i<=dims;i++) {
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2009-12-04 04:42:31 UTC (rev 16066)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2009-12-04 17:44:24 UTC (rev 16067)
@@ -674,18 +674,15 @@
if topvbc=1, will apply to velocities
if topvbc=0, will apply to tractions
-
-if is_inetrnal is set, will check for toplayerbc > 0, and assign
-internal boundary conditions to nodes within layer <= toplayerbc
*/
-void ggrd_read_vtop_from_file(struct All_variables *E, int is_global, int is_internal)
+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_inmsg, mpi_success_message = 1;
int m,el,i,k,i1,i2,ind,nodel,j,level, verbose;
- int noxg,nozg,noyg,noxl,noyl,nozl,lselect,idim,noxgnozg,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],
@@ -694,25 +691,28 @@
static pole_warned = FALSE, mode_warned = FALSE;
static ggrd_boolean shift_to_pos_lon = FALSE;
const int dims=E->mesh.nsd;
- int top_echo,nfree,nfixed,use_vel;
+ int top_proc,nfree,nfixed,use_vel,allow_internal;
#ifdef USE_GZDIR
gzFile *fp1;
#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;
+ /* number of nodes for this processor at highest MG level */
+ nox = E->lmesh.nox;
+ noz = E->lmesh.noz;
+ noy = E->lmesh.noy;
+ noxnoz = nox*noz;
+ if(E->mesh.toplayerbc > 0)
+ allow_internal = TRUE;
+ else
+ allow_internal = FALSE;
-
/* top processor check */
- top_echo = E->parallel.nprocz-1;
+ top_proc = E->parallel.nprocz-1;
/* output of warning messages */
- if((is_internal && (E->parallel.me == 0))||(E->parallel.me == top_echo))
+ if((allow_internal && (E->parallel.me == 0))||(E->parallel.me == top_proc))
verbose = TRUE;
else
verbose = FALSE;
@@ -729,7 +729,12 @@
break;
}
- /* read in plate code files? */
+ /*
+
+ read in plate code files? this will assign Euler vector
+ respective velocities
+
+ */
use_codes = (E->control.ggrd_vtop_omega[0] > 1e-7)?(1):(0);
save_codes = 0;
if(use_codes && (!use_vel)){
@@ -760,7 +765,7 @@
mode_warned = TRUE;
}
}
- if (is_internal || (E->parallel.me_loc[3] == top_echo)) {
+ if (allow_internal || (E->parallel.me_loc[3] == top_proc)) {
/*
top processors for regular operations, all for internal
@@ -845,7 +850,7 @@
snprintf(tfilename1,1000,"%s/codes.%d.gz", E->control.data_dir,E->parallel.me);
fp1 = gzdir_output_open(tfilename1,"w");
}
- if(E->parallel.me == top_echo)
+ if(verbose)
if(use_codes)
fprintf(stderr,"ggrd_read_vtop_from_file: assigning Euler vector %g, %g, %g to plates with code %i\n",
E->control.ggrd_vtop_omega[1],
@@ -898,11 +903,10 @@
fprintf(stderr,"ggrd_read_vtop_from_file: temporally constant velocity BC \n");
}
- if(verbose){
+ if(verbose)
fprintf(stderr,"ggrd_read_vtop_from_file: assigning %s BC, timedep: %i time: %g\n",
(use_vel)?("velocities"):("tractions"), timedep,age);
-
- }
+
/* if mixed BCs are allowed, need to reassign the boundary
condition */
if(E->control.ggrd_allow_mixed_vbcs){
@@ -914,15 +918,13 @@
*/
if(use_codes)
myerror(E,"cannot mix Euler velocities for plate codes and mixed vbcs");
- if((is_internal && (E->parallel.me==0)) || (E->parallel.me == top_echo))
+ if(verbose)
fprintf(stderr,"WARNING: allowing mixed velocity BCs\n");
-
-
/* velocities larger than the cutoff will be assigned as free
slip */
cutoff = E->control.ggrd.svp->fmaxlim[0] + 1e-5;
for(level=E->mesh.gridmax;level >= E->mesh.gridmin;level--){/* multigrid levels */
- /* all levels */
+ /* assign BCs to all levels */
noxl = E->lmesh.NOX[level];
noyl = E->lmesh.NOY[level];
nozl = E->lmesh.NOZ[level];
@@ -930,31 +932,28 @@
for (m=1;m <= E->sphere.caps_per_proc;m++) {
/* determine vertical nodes */
- if(is_internal){
+ if(allow_internal){
/* internal */
- if(E->mesh.toplayerbc == 0){ /* don't assign */
+ /* 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{ /* check for internal nodes in layers */
- for(k=nozl-1;k >= 1;k--){
- if(layers(E,m,k) > E->mesh.toplayerbc) /* assume regular mesh structure */
- break;
- }
- if(k == nozl-1){
- assign = FALSE;
- }else{
- assign = TRUE;
- topnode = nozl-1;
- botnode = k+1;
- }
+ }else{
+ assign = TRUE;topnode = nozl;botnode = k+1;
}
}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",
- is_internal,assign,botnode,topnode,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);
/*
- loop through all horizontal nodes
+ loop through all horizontal nodes and assign boundary
+ conditions for all required levels
*/
if(assign){
for(i=1;i <= noyl;i++){
@@ -988,15 +987,13 @@
if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],(E->control.ggrd.svp+i1),
vin1,FALSE,shift_to_pos_lon)){
fprintf(stderr,"ggrd_read_vtop_from_file: interpolation error at %g, %g\n",
- rout[1],rout[2]);
- parallel_process_termination();
+ rout[1],rout[2]);parallel_process_termination();
}
if(interpolate){ /* second time */
if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],(E->control.ggrd.svp+i2),vin2,
FALSE,shift_to_pos_lon)){
fprintf(stderr,"ggrd_read_mat_from_file: interpolation error at %g, %g\n",
- rout[1],rout[2]);
- parallel_process_termination();
+ rout[1],rout[2]);parallel_process_termination();
}
v[2] = (f1*vin1[0] + f2*vin2[0]); /* vphi unscaled! */
}else{
@@ -1051,40 +1048,38 @@
condition values
*/
+ /* scaled cutoff velocity */
+ if(!use_codes) /* else, is not defined */
+ cutoff = E->control.ggrd.svp->fmaxlim[0] * vscale + 1e-5;
+ else{
+ cutoff = 1e30;
+ if(save_codes) /* those will be surface nodes only */
+ gzprintf(fp1,"%3d %7d\n",m,E->lmesh.nsf);
+ }
+
for (m=1;m <= E->sphere.caps_per_proc;m++) {
- if(is_internal){
- if(E->mesh.toplayerbc == 0){ /* don't assign */
+ 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{ /* check for internal nodes in layers */
- for(k=nozg-1;k >= 1;k--){
- if(layers(E,m,k) > E->mesh.toplayerbc) /* assume regular mesh structure */
- break;
- }
- if(k == nozg-1){
- assign = FALSE;
- }else{
- assign = TRUE;topnode = nozg-1;botnode = k+1;
- }
+ }else{
+ assign = TRUE;topnode = noz;botnode = k+1;
}
}else{ /* just top node */
assign = TRUE;
- topnode = botnode = nozg;
+ topnode = botnode = noz;
}
- if(verbose)fprintf(stderr,"ggrd_read_vtop_from_file: internal: %i assign: %i k: %i to %i (%i)\n",
- is_internal,assign,botnode,topnode,nozg);
-
- /* scaled cutoff velocity */
- if(!use_codes) /* else, is not defined */
- cutoff = E->control.ggrd.svp->fmaxlim[0] * vscale + 1e-5;
- else{
- cutoff = 1e30;
- if(save_codes) /* those will be surface nodes only */
- gzprintf(fp1,"%3d %7d\n",m,E->lmesh.nsf);
- }
- if(assign){
- for(i=1;i <= noyg;i++) {/* loop through surface nodes */
- for(j=1;j <= noxg;j++) {
- nodel = nozg + (j-1) * nozg + (i-1)*noxgnozg; /* top node = nozg + (j-1) * nozg + (i-1)*noxgnozg; */
+ 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);
+ if(assign){
+ for(i=1;i <= noy;i++) {/* loop through surface nodes */
+ for(j=1;j <= nox;j++) {
+ nodel = noz + (j-1) * noz + (i-1)*noxnoz; /* 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];
@@ -1158,13 +1153,12 @@
*/
for(k = botnode;k <= topnode;k++){
- nodel = k + (j-1) * nozg + (i-1)*noxgnozg ; /* node = k + (j-1) * nozg + (i-1)*noxgnozg; */
+ 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 */
code = (int)(v[1] + 0.5);
if(save_codes) /* lon lat code */
- gzprintf(fp1, "%9.4f %9.4f %i\n",
- rout[2]/M_PI*180,90-rout[1]*180/M_PI,code);
+ gzprintf(fp1, "%9.4f %9.4f %i\n",rout[2]/M_PI*180,90-rout[1]*180/M_PI,code);
if((int)E->control.ggrd_vtop_omega[0] == code){
/* within plate */
sin_theta=sin(rout[1]);cos_theta=cos(rout[1]);
@@ -1183,7 +1177,6 @@
v[1] = v[2] = 0.0;
}
}
-
/* assign velociites */
if(fabs(v[2]) > cutoff){
/* huge velocitie - free slip */
@@ -1195,7 +1188,7 @@
E->sphere.cap[m].VB[2][nodel] = v[2]; /* phi */
}
E->sphere.cap[m].VB[3][nodel] = 0.0; /* r */
- }
+ } /* end z */
} /* end x */
} /* end y */
} /* end assign */
@@ -1211,7 +1204,9 @@
gzclose(fp1);
}
} /* end top proc branch */
- E->control.ggrd.vtop_control_init = 1;
+
+ /* */
+ E->control.ggrd.vtop_control_init = TRUE;
//if(E->parallel.me == 0)
//fprintf(stderr,"vtop from grd done: %i\n",lc++);
}
Modified: mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c 2009-12-04 04:42:31 UTC (rev 16066)
+++ mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c 2009-12-04 17:44:24 UTC (rev 16067)
@@ -35,8 +35,8 @@
#endif
/* ========================================== */
-
-static void horizontal_bc();
+void horizontal_bc(struct All_variables *,float *[],int,int,float,unsigned int,char,int,int);
+void assign_internal_bc(struct All_variables * ,int);
static void velocity_apply_periodic_bcs();
static void temperature_apply_periodic_bcs();
static void velocity_refl_vert_bc();
@@ -52,7 +52,7 @@
void velocity_imp_vert_bc();
void renew_top_velocity_boundary();
void apply_side_sbc();
- void assign_internal_bc(struct All_variables * );
+
int node,d,j,noz,lv,k;
for(lv=E->mesh.gridmax;lv>=E->mesh.gridmin;lv--)
@@ -66,11 +66,6 @@
horizontal_bc(E,E->sphere.cap[j].VB,noz,1,E->control.VBXtopval,SBX,1,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,noz,3,0.0,SBZ,0,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,noz,2,E->control.VBYtopval,SBY,1,lv,j);
-#ifdef USE_GGRD
- /* Ggrd traction control */
- if((lv==E->mesh.gridmax) && E->control.ggrd.vtop_control)
- ggrd_read_vtop_from_file(E, 0, 0);
-#endif
}
else if(E->mesh.topvbc == 1) {
horizontal_bc(E,E->sphere.cap[j].VB,noz,1,E->control.VBXtopval,VBX,1,lv,j);
@@ -127,8 +122,13 @@
anything at present, if E->mesh.toplayerbc > 0
*/
- assign_internal_bc(E);
+ assign_internal_bc(E,0);
+#ifdef USE_GGRD
+ if(E->control.ggrd.vtop_control)
+ ggrd_read_vtop_from_file(E,0);
+#endif
+
if(E->control.verbose) {
for (j=1;j<=E->sphere.caps_per_proc;j++)
for (node=1;node<=E->lmesh.nno;node++)
@@ -390,88 +390,6 @@
}
-/* ========================================================= */
-
-
-static void 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;
- static short int warned = FALSE;
-
- /* safety feature */
- if(dirn > E->mesh.nsd)
- return;
- 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, %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 = 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 */
- }
-
- /* 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)*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 */
- }
- }else{
- 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 = 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 */
- }
-
- /* 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)*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 */
- }
-
- return;
-}
-
-
static void velocity_apply_periodic_bcs(E)
struct All_variables *E;
{
@@ -495,57 +413,8 @@
}
-/*
-facility to apply internal velocity or stress conditions after
-top/bottom
-options:
-
-toplayerbc > 0: assign surface BC down to toplayerbc nd
-
- */
-void assign_internal_bc(struct All_variables *E)
-{
-
- int lv, j, noz, k,node,lay;
- /* stress or vel BC within a layer */
-
-
- 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){
- 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);
- }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);
- }
- }
- }
- }
-#ifdef USE_GGRD
- if(E->control.ggrd.vtop_control)
- ggrd_read_vtop_from_file(E, 0, 1);
-#endif
- }
-}
-
-
-
/* version */
/* $Id$ */
Modified: mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c 2009-12-04 04:42:31 UTC (rev 16066)
+++ mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c 2009-12-04 17:44:24 UTC (rev 16067)
@@ -227,9 +227,7 @@
case 1: /* velocity boundary conditions */
#ifdef USE_GGRD
- if(E->control.ggrd.vtop_control){
- ggrd_read_vtop_from_file(E, 0,0);
- }else{
+ if(!E->control.ggrd.vtop_control){ /* grd control is called from boundary condition file */
#endif
nnn=nox*noy;
for(i=1;i<=dims;i++) {
Modified: mc/3D/CitcomS/trunk/lib/ggrd_handling.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/ggrd_handling.h 2009-12-04 04:42:31 UTC (rev 16066)
+++ mc/3D/CitcomS/trunk/lib/ggrd_handling.h 2009-12-04 17:44:24 UTC (rev 16067)
@@ -39,7 +39,7 @@
void ggrd_temp_init_general(struct All_variables *,int);
void ggrd_read_mat_from_file(struct All_variables *, int );
void ggrd_read_ray_from_file(struct All_variables *, int );
-void ggrd_read_vtop_from_file(struct All_variables *, int ,int);
+void ggrd_read_vtop_from_file(struct All_variables *, int);
void ggrd_read_age_from_file(struct All_variables *, int );
void xyz2rtp(float ,float ,float ,float *);
void xyz2rtpd(float ,float ,float ,double *);
More information about the CIG-COMMITS
mailing list