[cig-commits] r20511 - mc/3D/CitcomCU/trunk/src
becker at geodynamics.org
becker at geodynamics.org
Tue Jul 10 03:09:17 PDT 2012
Author: becker
Date: 2012-07-10 03:09:17 -0700 (Tue, 10 Jul 2012)
New Revision: 20511
Modified:
mc/3D/CitcomCU/trunk/src/Boundary_conditions.c
mc/3D/CitcomCU/trunk/src/Composition_adv.c
mc/3D/CitcomCU/trunk/src/Convection.c
mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
mc/3D/CitcomCU/trunk/src/Instructions.c
mc/3D/CitcomCU/trunk/src/Process_velocity.c
mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
mc/3D/CitcomCU/trunk/src/global_defs.h
mc/3D/CitcomCU/trunk/src/prototypes.h
Log:
Fixed missing "break" in Viscosity_structures, and minor adjustments
Modified: mc/3D/CitcomCU/trunk/src/Boundary_conditions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Boundary_conditions.c 2012-07-09 10:34:28 UTC (rev 20510)
+++ mc/3D/CitcomCU/trunk/src/Boundary_conditions.c 2012-07-10 10:09:17 UTC (rev 20511)
@@ -116,7 +116,8 @@
if(lv == E->mesh.levmax){ /* NB */
E->VB[1][node] = 0; /* set to zero */
}
- }else{ /* assume only x or y periodic */
+ }
+ if(E->mesh.periodic_y){
E->NODE[lv][node] = E->NODE[lv][node] & (~SBY);
E->NODE[lv][node] = E->NODE[lv][node] | (VBY);
if(lv == E->mesh.levmax){ /* NB */
@@ -126,6 +127,11 @@
}
}
}
+
+ if(E->mesh.slab_influx_side_bc) /* slab in-flux on side */
+ velocity_apply_slab_influx_side_bc(E);
+
+
if(E->control.verbose)
{
for(node = 1; node <= E->lmesh.nno; node++)
@@ -155,7 +161,7 @@
int lev,top;
if(E->parallel.me == 0)
- fprintf(stderr,"WARNING: freezing surface boundary condition at time step %i\n",
+ fprintf(stderr,"WARNING: freezing surface boundary condition at time step %i (not working yet)\n",
E->monitor.solution_cycles);
/*
@@ -237,121 +243,186 @@
/* except one side with XOZ and y=0, all others are not reflecting BC*/
/* for two YOZ planes if 3-D, or two OZ side walls for 2-D */
-
- if (E->parallel.me_loc[1]==0 || E->parallel.me_loc[1]==E->parallel.nprocx-1)
+ if (((!E->mesh.slab_influx_side_bc)&&(E->parallel.me_loc[1]==0)) ||
+ (E->parallel.me_loc[1]==E->parallel.nprocx-1)){
for(j=1;j<=E->lmesh.noy;j++)
for(i=1;i<=E->lmesh.noz;i++) {
- node1 = i + (j-1)*E->lmesh.noz*E->lmesh.nox;
+ node1 = i + (j-1)*E->lmesh.noz*E->lmesh.nox;
node2 = node1 + (E->lmesh.nox-1)*E->lmesh.noz;
-
+
ii = i + E->lmesh.nzs - 1;
- if (E->parallel.me_loc[1]==0 ) {
+ if (E->parallel.me_loc[1]==0 ) {
E->VB[1][node1] = 0.0;
if((ii != 1) && (ii != E->mesh.noz))
- E->VB[3][node1] = 0.0;
- }
- if (E->parallel.me_loc[1]==E->parallel.nprocx-1) {
+ E->VB[3][node1] = 0.0;
+ }
+ if (E->parallel.me_loc[1]==E->parallel.nprocx-1) {
E->VB[1][node2] = 0.0;
if((ii != 1) && (ii != E->mesh.noz))
- E->VB[3][node2] = 0.0;
- }
- } /* end loop for i and j */
-
+ E->VB[3][node2] = 0.0;
+ }
+ } /* end loop for i and j */
+ }
/* for two XOZ planes if 3-D */
- if (E->parallel.me_loc[2]==0 || E->parallel.me_loc[2]==E->parallel.nprocy-1)
- for(j=1;j<=E->lmesh.nox;j++)
- for(i=1;i<=E->lmesh.noz;i++) {
- node1 = i +(j-1)*E->lmesh.noz;
- node2 = node1+(E->lmesh.noy-1)*E->lmesh.noz*E->lmesh.nox;
- ii = i + E->lmesh.nzs - 1;
+ if ((E->parallel.me_loc[2]==0) || (E->parallel.me_loc[2]==E->parallel.nprocy-1)){
+ for(j=1;j<=E->lmesh.nox;j++)
+ for(i=1;i<=E->lmesh.noz;i++) {
+ node1 = i +(j-1)*E->lmesh.noz;
+ node2 = node1+(E->lmesh.noy-1)*E->lmesh.noz*E->lmesh.nox;
+ ii = i + E->lmesh.nzs - 1;
- if (E->parallel.me_loc[2]==0) {
- E->VB[2][node1] = 0.0;
- if((ii != 1) && (ii != E->mesh.noz))
- E->VB[3][node1] = 0.0;
- }
- if (E->parallel.me_loc[2]==E->parallel.nprocy-1) {
- E->VB[2][node2] = 0.0;
- if((ii != 1) && (ii != E->mesh.noz))
- E->VB[3][node2] = 0.0;
- }
- } /* end of loop i & j */
-
+ if (E->parallel.me_loc[2]==0) {
+ E->VB[2][node1] = 0.0;
+ if((ii != 1) && (ii != E->mesh.noz))
+ E->VB[3][node1] = 0.0;
+ }
+ if (E->parallel.me_loc[2]==E->parallel.nprocy-1) {
+ E->VB[2][node2] = 0.0;
+ if((ii != 1) && (ii != E->mesh.noz))
+ E->VB[3][node2] = 0.0;
+ }
+ } /* end of loop i & j */
+ }
/* all vbc's apply at all levels */
for(level=E->mesh.levmax;level>=E->mesh.levmin;level--) {
nox = E->lmesh.NOX[level] ;
noz = E->lmesh.NOZ[level] ;
noy = E->lmesh.NOY[level] ;
-
- if (E->parallel.me_loc[1]==0 || E->parallel.me_loc[1]==E->parallel.nprocx-1)
+ if (((!E->mesh.slab_influx_side_bc)&&(E->parallel.me_loc[1]==0)) ||
+ (E->parallel.me_loc[1]==E->parallel.nprocx-1)){
for(j=1;j<=noy;j++)
- for(i=1;i<=noz;i++) {
+ for(i=1;i<=noz;i++) {
node1 = i + (j-1)*noz*nox ;
node2 = node1 + (nox-1) * noz ;
ii = i + E->lmesh.NZS[level] - 1;
- if (E->parallel.me_loc[1]==0 ) {
- E->NODE[level][node1] = E->NODE[level][node1] & (~SBX);
- E->NODE[level][node1] = E->NODE[level][node1] | (VBX);
- if((ii!=1) && (ii!=E->mesh.NOZ[level])) {
- E->NODE[level][node1] = E->NODE[level][node1] & (~VBY);
- E->NODE[level][node1] = E->NODE[level][node1] | SBY;
- E->NODE[level][node1] = E->NODE[level][node1] & (~ VBZ);
- E->NODE[level][node1] = E->NODE[level][node1] | SBZ;
- }
- }
- if (E->parallel.me_loc[1]==E->parallel.nprocx-1) {
- E->NODE[level][node2] = E->NODE[level][node2] & (~SBX);
- E->NODE[level][node2] = E->NODE[level][node2] | (VBX);
- if((ii!=1) && (ii!=E->mesh.NOZ[level])) {
- E->NODE[level][node2] = E->NODE[level][node2] & (~VBY);
- E->NODE[level][node2] = E->NODE[level][node2] | SBY;
- E->NODE[level][node2] = E->NODE[level][node2] & (~ VBZ);
- E->NODE[level][node2] = E->NODE[level][node2] | SBZ;
- }
- }
- } /* end for loop i & j */
-
- if (E->parallel.me_loc[2]==0 || E->parallel.me_loc[2]==E->parallel.nprocy-1)
- for(j=1;j<=nox;j++)
- for(i=1;i<=noz;i++) {
- node1 = i + (j-1)*noz;
- node2 = node1+(noy-1)*noz*nox;
- ii = i + E->lmesh.NZS[level] - 1;
- jj = j + E->lmesh.NXS[level] - 1;
- if (E->parallel.me_loc[2]==0) {
- E->NODE[level][node1] = E->NODE[level][node1] | VBY;
- E->NODE[level][node1] = E->NODE[level][node1] & (~SBY);
- if((ii!= 1) && (ii != E->mesh.NOZ[level])) {
- E->NODE[level][node1] = E->NODE[level][node1] & (~VBZ);
- E->NODE[level][node1] = E->NODE[level][node1] | SBZ;
- }
- if((jj!=1) && (jj!=E->mesh.NOX[level]) && (ii!=1) && (ii!=E->mesh.NOZ[level])){
- E->NODE[level][node1] = E->NODE[level][node1] & (~VBX);
- E->NODE[level][node1] = E->NODE[level][node1] | SBX;
- }
- }
- if (E->parallel.me_loc[2]==E->parallel.nprocy-1) {
- E->NODE[level][node2] = E->NODE[level][node2] | VBY;
- E->NODE[level][node2] = E->NODE[level][node2] & (~SBY);
- if((ii!= 1) && (ii != E->mesh.NOZ[level])) {
- E->NODE[level][node2] = E->NODE[level][node2] & (~VBZ);
- E->NODE[level][node2] = E->NODE[level][node2] | SBZ;
- }
- if((jj!=1) && (jj!=E->mesh.NOX[level]) && (ii!=1) && (ii!=E->mesh.NOZ[level])){
- E->NODE[level][node2] = E->NODE[level][node2] & (~VBX);
- E->NODE[level][node2] = E->NODE[level][node2] | SBX;
- }
- }
-
- } /* end for loop i & j */
+ if (E->parallel.me_loc[1]==0 ) {
+ E->NODE[level][node1] = E->NODE[level][node1] & (~SBX);
+ E->NODE[level][node1] = E->NODE[level][node1] | (VBX);
+ if((ii!=1) && (ii!=E->mesh.NOZ[level])) {
+ E->NODE[level][node1] = E->NODE[level][node1] & (~VBY);
+ E->NODE[level][node1] = E->NODE[level][node1] | SBY;
+ E->NODE[level][node1] = E->NODE[level][node1] & (~ VBZ);
+ E->NODE[level][node1] = E->NODE[level][node1] | SBZ;
+ }
+ }
+ if (E->parallel.me_loc[1]==E->parallel.nprocx-1) {
+ E->NODE[level][node2] = E->NODE[level][node2] & (~SBX);
+ E->NODE[level][node2] = E->NODE[level][node2] | (VBX);
+ if((ii!=1) && (ii!=E->mesh.NOZ[level])) {
+ E->NODE[level][node2] = E->NODE[level][node2] & (~VBY);
+ E->NODE[level][node2] = E->NODE[level][node2] | SBY;
+ E->NODE[level][node2] = E->NODE[level][node2] & (~ VBZ);
+ E->NODE[level][node2] = E->NODE[level][node2] | SBZ;
+ }
+ }
+ } /* end for loop i & j */
+ }
+ if ((E->parallel.me_loc[2]==0) || (E->parallel.me_loc[2]==E->parallel.nprocy-1)){
+ for(j=1;j<=nox;j++)
+ for(i=1;i<=noz;i++) {
+ node1 = i + (j-1)*noz;
+ node2 = node1+(noy-1)*noz*nox;
+ ii = i + E->lmesh.NZS[level] - 1;
+ jj = j + E->lmesh.NXS[level] - 1;
+ if (E->parallel.me_loc[2]==0) {
+ E->NODE[level][node1] = E->NODE[level][node1] | VBY;
+ E->NODE[level][node1] = E->NODE[level][node1] & (~SBY);
+ if((ii!= 1) && (ii != E->mesh.NOZ[level])) {
+ E->NODE[level][node1] = E->NODE[level][node1] & (~VBZ);
+ E->NODE[level][node1] = E->NODE[level][node1] | SBZ;
+ }
+ if((jj!=1) && (jj!=E->mesh.NOX[level]) && (ii!=1) && (ii!=E->mesh.NOZ[level])){
+ E->NODE[level][node1] = E->NODE[level][node1] & (~VBX);
+ E->NODE[level][node1] = E->NODE[level][node1] | SBX;
+ }
+ }
+ if (E->parallel.me_loc[2]==E->parallel.nprocy-1) {
+ E->NODE[level][node2] = E->NODE[level][node2] | VBY;
+ E->NODE[level][node2] = E->NODE[level][node2] & (~SBY);
+ if((ii!= 1) && (ii != E->mesh.NOZ[level])) {
+ E->NODE[level][node2] = E->NODE[level][node2] & (~VBZ);
+ E->NODE[level][node2] = E->NODE[level][node2] | SBZ;
+ }
+ if((jj!=1) && (jj!=E->mesh.NOX[level]) && (ii!=1) && (ii!=E->mesh.NOZ[level])){
+ E->NODE[level][node2] = E->NODE[level][node2] & (~VBX);
+ E->NODE[level][node2] = E->NODE[level][node2] | SBX;
+ }
+ }
+
+ } /* end for loop i & j */
+ }
} /* end for loop level */
-
+
return;
}
+void velocity_apply_slab_influx_side_bc(struct All_variables *E)
+{
+
+ int i,j,node,level,nox,noy,noz;
+ float vminus,d1,d2;
+ d1 = 1.0 - E->mesh.slab_influx_z2;
+ d2 = E->mesh.slab_influx_z2 - E->mesh.slab_influx_z1;
+ vminus = -2*(d1+d2/2.)/(1-d1-d2)*E->control.sub_vel;
+ if(E->parallel.me == 0)
+ fprintf(stderr,"WARNING: using experimental side influx boundary condition, v0: %g v-: %g z1: %g z2: %g\n",
+ E->control.sub_vel,vminus, E->mesh.slab_influx_z1, E->mesh.slab_influx_z2);
+
+ if(E->mesh.periodic_x)
+ myerror("error, need reflective boundary conditions on X for slab influx",E);
+ if(E->control.Rsphere)
+ myerror("error, need Cartesian geoemetry for slab influx",E);
+ if (E->parallel.me_loc[1]==0){
+ /*
+ assign velocity values
+ */
+ for(j=1;j<=E->lmesh.noy;j++){
+ for(i=1;i<=E->lmesh.noz;i++) {
+ node = i + (j-1)*E->lmesh.noz*E->lmesh.nox;
+ //if(E->X[2][node] <= E->mesh.slab_influx_y1){
+ /* apply VBc for slab */
+ if(E->X[3][node] >= E->mesh.slab_influx_z2)
+ E->VB[1][node] = E->control.sub_vel; /* plate velocity above z2 */
+ else if(E->X[3][node] >= E->mesh.slab_influx_z1){
+ E->VB[1][node] = E->control.sub_vel * (E->X[3][node] - E->mesh.slab_influx_z1)/d2; /* tapered down to z1 */
+ }else{
+ E->VB[1][node] = vminus * (E->X[3][node]/E->mesh.slab_influx_z1);
+ }
+
+ /* y and z directions fixed */
+ E->VB[2][node] = 0.0;
+ E->VB[3][node] = 0.0;
+ }
+ }
+ /*
+
+ set velocity boundary conditions at all levels on the complete left hand size
+
+ */
+ for(level=E->mesh.levmax;level>=E->mesh.levmin;level--) {
+ nox = E->lmesh.NOX[level] ;
+ noz = E->lmesh.NOZ[level] ;
+ noy = E->lmesh.NOY[level] ;
+ for(j=1;j<=noy;j++)
+ for(i=1;i<=noz;i++) {
+ node = i + (j-1)*noz*nox ;
+ /* no slip all around */
+ E->NODE[level][node] = E->NODE[level][node] | (VBX);
+ E->NODE[level][node] = E->NODE[level][node] & (~SBX);
+ E->NODE[level][node] = E->NODE[level][node] | (VBY);
+ E->NODE[level][node] = E->NODE[level][node] & (~SBY);
+
+ E->NODE[level][node] = E->NODE[level][node] | (VBZ);
+ E->NODE[level][node] = E->NODE[level][node] & (~SBZ);
+
+ }
+ }
+ } /* end only left-most processors branch */
+}
+
void temperature_refl_vert_bc(struct All_variables *E)
{
int i, j;
@@ -548,7 +619,7 @@
}
/* all vbc's apply at all levels */
-for(level=E->mesh.levmax;level>=E->mesh.levmin;level--) {
+ for(level=E->mesh.levmax;level>=E->mesh.levmin;level--) {
nox = E->lmesh.NOX[level] ;
noz = E->lmesh.NOZ[level] ;
noy = E->lmesh.NOY[level] ;
Modified: mc/3D/CitcomCU/trunk/src/Composition_adv.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Composition_adv.c 2012-07-09 10:34:28 UTC (rev 20510)
+++ mc/3D/CitcomCU/trunk/src/Composition_adv.c 2012-07-10 10:09:17 UTC (rev 20511)
@@ -543,6 +543,11 @@
{
C[i] = 0.0;
}
+
+ /* this will assign C=1 on top left side */
+ if(E->mesh.slab_influx_side_bc)
+ composition_apply_slab_influx_side_bc(E);
+
/* for each element, count dense and regular marks */
for(imark = 1; imark <= E->advection.markers; imark++)
Modified: mc/3D/CitcomCU/trunk/src/Convection.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Convection.c 2012-07-09 10:34:28 UTC (rev 20510)
+++ mc/3D/CitcomCU/trunk/src/Convection.c 2012-07-10 10:09:17 UTC (rev 20511)
@@ -505,61 +505,67 @@
// for composition
- if(E->control.composition && (E->control.restart == 1 || E->control.restart == 2))
- {
- convection_initial_markers(E,0);
- }
-
- else if(E->control.composition && (E->control.restart == 3))
- {
- sprintf(output_file, "%s.comp.%d.%d", E->convection.old_T_file, E->parallel.me, E->monitor.solution_cycles);
- fp = fopen(output_file, "r");
- fgets(input_s, 200, fp);
- sscanf(input_s, "%d %d %g", &i, &j, &E->monitor.elapsed_time);
-
- for(node = 1; node <= E->lmesh.nno; node++)
+ if(E->control.composition){
+
+ if(E->control.restart == 1 || E->control.restart == 2)
+ {
+ convection_initial_markers(E,0);
+ }
+
+ else if(E->control.restart == 3)
+ {
+ sprintf(output_file, "%s.comp.%d.%d", E->convection.old_T_file, E->parallel.me, E->monitor.solution_cycles);
+ fp = fopen(output_file, "r");
+ fgets(input_s, 200, fp);
+ sscanf(input_s, "%d %d %g", &i, &j, &E->monitor.elapsed_time);
+
+ for(node = 1; node <= E->lmesh.nno; node++)
{
- fgets(input_s, 200, fp);
- sscanf(input_s, "%g", &E->C[node]);
+ fgets(input_s, 200, fp);
+ sscanf(input_s, "%g", &E->C[node]);
}
- fclose(fp);
-
+ fclose(fp);
+
convection_initial_markers1(E);
- }
-
- else if(E->control.composition && E->control.restart == -1)
- {
-
- sprintf(output_file, "%s.traces.%d", E->convection.old_T_file, E->parallel.me);
- fp = fopen(output_file, "r");
- fgets(input_s, 200, fp);
- sscanf(input_s, "%d %d %g", &E->advection.markers, &j, &temp1);
- for(i = 1; i <= E->advection.markers; i++)
+ }
+
+ else if(E->control.restart == -1)
+ {
+
+ sprintf(output_file, "%s.traces.%d", E->convection.old_T_file, E->parallel.me);
+ fp = fopen(output_file, "r");
+ fgets(input_s, 200, fp);
+ sscanf(input_s, "%d %d %g", &E->advection.markers, &j, &temp1);
+ for(i = 1; i <= E->advection.markers; i++)
{
- fgets(input_s, 200, fp);
- sscanf(input_s, "%lf %lf %lf %d %d", &E->XMC[1][i], &E->XMC[2][i], &E->XMC[3][i], &E->CElement[i], &E->C12[i]);
- if(E->XMC[3][i] < E->XP[3][1])
- E->XMC[3][i] = E->XP[3][1];
- else if(E->XMC[3][i] > E->XP[3][E->lmesh.noz])
- E->XMC[3][i] = E->XP[3][E->lmesh.noz];
- if(E->XMC[2][i] < E->XP[2][1])
- E->XMC[2][i] = E->XP[2][1];
- else if(E->XMC[2][i] > E->XP[2][E->lmesh.noy])
- E->XMC[2][i] = E->XP[2][E->lmesh.noy];
- if(E->XMC[1][i] < E->XP[1][1])
- E->XMC[1][i] = E->XP[1][1];
- else if(E->XMC[1][i] > E->XP[1][E->lmesh.nox])
- E->XMC[1][i] = E->XP[1][E->lmesh.nox];
+ fgets(input_s, 200, fp);
+ sscanf(input_s, "%lf %lf %lf %d %d", &E->XMC[1][i], &E->XMC[2][i], &E->XMC[3][i], &E->CElement[i], &E->C12[i]);
+ if(E->XMC[3][i] < E->XP[3][1])
+ E->XMC[3][i] = E->XP[3][1];
+ else if(E->XMC[3][i] > E->XP[3][E->lmesh.noz])
+ E->XMC[3][i] = E->XP[3][E->lmesh.noz];
+ if(E->XMC[2][i] < E->XP[2][1])
+ E->XMC[2][i] = E->XP[2][1];
+ else if(E->XMC[2][i] > E->XP[2][E->lmesh.noy])
+ E->XMC[2][i] = E->XP[2][E->lmesh.noy];
+ if(E->XMC[1][i] < E->XP[1][1])
+ E->XMC[1][i] = E->XP[1][1];
+ else if(E->XMC[1][i] > E->XP[1][E->lmesh.nox])
+ E->XMC[1][i] = E->XP[1][E->lmesh.nox];
}
- for(i = 1; i <= E->lmesh.nel; i++)
+ for(i = 1; i <= E->lmesh.nel; i++)
{
- fgets(input_s, 200, fp);
- sscanf(input_s, "%g", &E->CE[i]);
+ fgets(input_s, 200, fp);
+ sscanf(input_s, "%g", &E->CE[i]);
}
- fclose(fp);
+ fclose(fp);
-// get_C_from_markers(E,E->C);
-
+ // get_C_from_markers(E,E->C);
+
+ }
+ /* this will assign C=1 on top left side */
+ if(E->mesh.slab_influx_side_bc)
+ composition_apply_slab_influx_side_bc(E);
}
E->advection.timesteps = E->monitor.solution_cycles;
@@ -1001,3 +1007,21 @@
free((void *)P2);
return;
}
+
+/*
+ assign C=1 composition at left half of box, on top
+*/
+void composition_apply_slab_influx_side_bc(struct All_variables *E)
+{
+ int i;
+ if(E->control.composition){
+ for(i=1;i<=E->advection.markers;i++){
+ if((E->XMC[1][i] <= 0.2)&&
+ (E->XMC[2][i] <= E->mesh.slab_influx_y2)&&
+ (E->XMC[3][i] >= E->mesh.slab_influx_z2)){
+ E->C12[i] = 1;
+ }
+ }
+ }
+
+}
Modified: mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2012-07-09 10:34:28 UTC (rev 20510)
+++ mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2012-07-10 10:09:17 UTC (rev 20511)
@@ -1029,7 +1029,7 @@
nr.grd, nt.grd, np.grd for the directors
*/
-void ggrd_read_anivisc_from_file(struct All_variables *E)
+void ggrd_read_anivisc_from_file_cu(struct All_variables *E)
{
MPI_Status mpi_stat;
int mpi_rc;
@@ -1074,7 +1074,11 @@
}
}
}
- sprintf(gmt_string,""); /* regional */
+ if(E->control.Rsphere)
+ sprintf(gmt_string,GGRD_GMT_GEOGRAPHIC_STRING);
+ else
+ sprintf(gmt_string,"");
+
/*
Modified: mc/3D/CitcomCU/trunk/src/Instructions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Instructions.c 2012-07-09 10:34:28 UTC (rev 20510)
+++ mc/3D/CitcomCU/trunk/src/Instructions.c 2012-07-10 10:09:17 UTC (rev 20511)
@@ -494,10 +494,19 @@
E->mesh.topvbc = 0; /* stress */
E->mesh.botvbc = 0;
E->mesh.sidevbc = 0;
+
E->mesh.periodic_x = 0; /* reflection is default */
E->mesh.periodic_y = 0;
E->mesh.periodic_pin_or_filter = 0; /* filter by default */
-
+
+ E->mesh.slab_influx_side_bc = 0;
+
+ E->mesh.slab_influx_z1=0.6; /* above z1, below z2 --> tapered inflow */
+ E->mesh.slab_influx_z2=0.9; /* above z2 --> lithosphere */
+
+ E->mesh.slab_influx_y2=2.37; /* y-extent of continental BC */
+
+
E->control.VBXtopval = 0.0;
E->control.VBYtopval = 0.0;
E->control.VBXbotval = 0.0;
@@ -930,7 +939,9 @@
input_boolean("periodicx", &(E->mesh.periodic_x), "off", m);
input_boolean("periodicy", &(E->mesh.periodic_y), "off", m);
- input_int("periodic_pin_or_filter", &(E->mesh.periodic_pin_or_filter), "0", m); /* 1: pin 0: filter */
+ input_int("periodic_pin_or_filter", &(E->mesh.periodic_pin_or_filter), "0", m); /* 1: pin
+ 0: filter */
+ input_boolean("slab_influx_side_bc",&(E->mesh.slab_influx_side_bc),"off",m);
input_boolean("depthdominated", &(E->control.depth_dominated), "off", m);
input_boolean("eqnzigzag", &(E->control.eqn_zigzag), "off", m);
@@ -952,6 +963,7 @@
if((E->parallel.me == 0) && (fabs(E->control.plate_vel) > 5e-7))
fprintf(stderr,"WARNING: plate velocity is overriding VBx \n");
+ input_float("sub_vel", &(E->control.sub_vel), "0.0", m);
input_float("plate_age", &(E->control.plate_age), "0.0", m);
input_float("plume_radius", &(E->segment.plume_radius), "0.0", m);
Modified: mc/3D/CitcomCU/trunk/src/Process_velocity.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Process_velocity.c 2012-07-09 10:34:28 UTC (rev 20510)
+++ mc/3D/CitcomCU/trunk/src/Process_velocity.c 2012-07-10 10:09:17 UTC (rev 20511)
@@ -221,25 +221,24 @@
void remove_net_drift(struct All_variables *E)
{
- float net_drift;
- int coord,i;
+ float net_drift[3]={0.,0.,0.};
+ int i;
- if(E->mesh.periodic_x && E->mesh.periodic_y)
- myerror("remove_net_drift: can't deal with both x and y periodicity",E);
+
if(E->mesh.periodic_x)
- coord = 1;
- else if(E->mesh.periodic_y)
- coord = 2;
+ net_drift[1] = return_bulk_value(E, E->V[1], -1,1);
+ if(E->mesh.periodic_y)
+ net_drift[2] = return_bulk_value(E, E->V[2], -1,1);
+
- net_drift = return_bulk_value(E, E->V[coord], -1,1);
-
if(E->parallel.me == 0)
- fprintf(stderr,"remove_net_drift: removing net drift of %g in %i coord\n",
- net_drift,coord);
+ fprintf(stderr,"remove_net_drift: removing net drift of x: %g y: %g periodic: x: %i y: %i\n",
+ net_drift[1],net_drift[2],E->mesh.periodic_x,E->mesh.periodic_y);
+ for(i = 1; i <= E->lmesh.nno; i++){
+ E->V[1][i] -= net_drift[1];
+ E->V[2][i] -= net_drift[2];
+ }
- for(i = 1; i <= E->lmesh.nno; i++)
- E->V[coord][i] -= net_drift;
-
}
Modified: mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2012-07-09 10:34:28 UTC (rev 20510)
+++ mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2012-07-10 10:09:17 UTC (rev 20511)
@@ -546,6 +546,7 @@
EEta[(i - 1) * vpts + jj] = tempa * exp((E->viscosity.E[l]) / (temp + E->viscosity.T[l]));
}
}
+ break;
case 2:
/*
eta = eta0 * exp(E + (1-z)*Z0/(T+T0))
@@ -677,7 +678,7 @@
default:
myerror("RHEOL option undefined",E);
break;
- } /* end swith */
+ } /* end switch */
visits++;
Modified: mc/3D/CitcomCU/trunk/src/global_defs.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/global_defs.h 2012-07-09 10:34:28 UTC (rev 20510)
+++ mc/3D/CitcomCU/trunk/src/global_defs.h 2012-07-10 10:09:17 UTC (rev 20511)
@@ -578,6 +578,9 @@
int periodic_x;
int periodic_y;
+ int slab_influx_side_bc;
+ float slab_influx_z1,slab_influx_z2,slab_influx_y2;
+
int periodic_pin_or_filter; /* */
double volume;
float layer[4]; /* dimensionless dimensions */
@@ -794,6 +797,7 @@
int comparison;
int crust;
float plate_vel;
+ float sub_vel;
float plate_age;
//float tole_comp;
Modified: mc/3D/CitcomCU/trunk/src/prototypes.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/prototypes.h 2012-07-09 10:34:28 UTC (rev 20510)
+++ mc/3D/CitcomCU/trunk/src/prototypes.h 2012-07-10 10:09:17 UTC (rev 20511)
@@ -398,3 +398,6 @@
int weak_zones(struct All_variables *, int, float);
float boundary_thickness(struct All_variables *, float *);
void remove_net_drift(struct All_variables *);
+void velocity_apply_slab_influx_side_bc(struct All_variables *);
+void composition_apply_slab_influx_side_bc(struct All_variables *);
+
More information about the CIG-COMMITS
mailing list