[cig-commits] r18297 - mc/3D/CitcomCU/trunk/src
becker at geodynamics.org
becker at geodynamics.org
Wed Apr 27 05:19:44 PDT 2011
Author: becker
Date: 2011-04-27 05:19:44 -0700 (Wed, 27 Apr 2011)
New Revision: 18297
Modified:
mc/3D/CitcomCU/trunk/src/Composition_adv.c
mc/3D/CitcomCU/trunk/src/Convection.c
mc/3D/CitcomCU/trunk/src/Drive_solvers.c
mc/3D/CitcomCU/trunk/src/Parallel_related.c
mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
mc/3D/CitcomCU/trunk/src/global_defs.h
Log:
Added foundations for advection tracer flavors, not initialized or used for
anything useful yet
Added option to specify the number of maximum powerlaw iterations
Modified: mc/3D/CitcomCU/trunk/src/Composition_adv.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Composition_adv.c 2011-04-26 19:02:48 UTC (rev 18296)
+++ mc/3D/CitcomCU/trunk/src/Composition_adv.c 2011-04-27 12:19:44 UTC (rev 18297)
@@ -163,10 +163,10 @@
E->parallel.traces_transfer_index[neighbor] = (int *)malloc((markers + 1) * sizeof(int));
E->RVV[neighbor] = (float *)malloc(asize * sizeof(int));
E->RXX[neighbor] = (double *)malloc(asize * sizeof(double));
- E->RINS[neighbor] = (int *)malloc((markers + 1) * 2 * sizeof(int));
+ E->RINS[neighbor] = (int *)malloc((markers + 1) * (2 + E->tracers_add_flavors) * sizeof(int));
E->PVV[neighbor] = (float *)malloc(asize * sizeof(int));
E->PXX[neighbor] = (double *)malloc(asize * sizeof(double));
- E->PINS[neighbor] = (int *)malloc((markers + 1) * 2 * sizeof(int));
+ E->PINS[neighbor] = (int *)malloc((markers + 1) * (2 + E->tracers_add_flavors) * sizeof(int));
}
E->traces_leave_index = (int *)malloc((markers + 1) * sizeof(int));
been++;
@@ -266,11 +266,12 @@
void unify_markers_array(struct All_variables *E, int no_tran, int no_recv)
{
- int i, j;
+ int i, j,k;
int ii, jj, kk;
int nsd2, neighbor, no_trans1;
-
+ int rioff;
nsd2 = E->mesh.nsd * 2;
+ rioff = 2 + E->tracers_add_flavors;
ii = 0;
for(i = 1; i <= E->advection.markers; i++)
@@ -314,8 +315,10 @@
E->XMCpred[1][ii] = E->RXX[neighbor][j * nsd2 + 3];
E->XMCpred[2][ii] = E->RXX[neighbor][j * nsd2 + 4];
E->XMCpred[3][ii] = E->RXX[neighbor][j * nsd2 + 5];
- E->C12[ii] = E->RINS[neighbor][j * 2];
- E->CElement[ii] = E->RINS[neighbor][j * 2 + 1];
+ E->C12[ii] = E->RINS[neighbor][j * rioff];
+ E->CElement[ii] = E->RINS[neighbor][j * rioff + 1];
+ for(k=0;k < E->tracers_add_flavors;k++)
+ E->tflavors[ii*E->tracers_add_flavors+k] = E->RINS[neighbor][j * rioff + 2 + k];
E->traces_leave[ii] = 0;
}
}
@@ -341,8 +344,10 @@
E->XMCpred[1][ii] = E->RXX[neighbor][j * nsd2 + 3];
E->XMCpred[2][ii] = E->RXX[neighbor][j * nsd2 + 4];
E->XMCpred[3][ii] = E->RXX[neighbor][j * nsd2 + 5];
- E->C12[ii] = E->RINS[neighbor][j * 2];
- E->CElement[ii] = E->RINS[neighbor][j * 2 + 1];
+ E->C12[ii] = E->RINS[neighbor][j * rioff];
+ E->CElement[ii] = E->RINS[neighbor][j * rioff + 1];
+ for(k=0;k < E->tracers_add_flavors;k++)
+ E->tflavors[ii*E->tracers_add_flavors+k] = E->RINS[neighbor][j * rioff + 2 + k];
E->traces_leave[ii] = 0;
}
}
@@ -372,8 +377,11 @@
E->XMCpred[1][ii] = E->XMCpred[1][i];
E->XMCpred[2][ii] = E->XMCpred[2][i];
E->XMCpred[3][ii] = E->XMCpred[3][i];
- E->C12[ii] = E->C12[i];
+ E->C12[ii] = E->C12[i];
E->CElement[ii] = E->CElement[i];
+ for(k=0;k < E->tracers_add_flavors;k++)
+ E->tflavors[ii*E->tracers_add_flavors+k] = E->tflavors[i*E->tracers_add_flavors+k];
+
E->traces_leave[ii] = 0;
}
}
@@ -394,7 +402,7 @@
void prepare_transfer_arrays(struct All_variables *E)
{
- int j, part, neighbor, k1, k2, k3;
+ int j, part, neighbor, k1, k2, k3,k;
//if(E->parallel.me==0)fprintf(stderr,"ta 1 ok\n");
for(neighbor = 1; neighbor <= E->parallel.no_neighbors; neighbor++)
{
@@ -422,6 +430,9 @@
E->PINS[neighbor][k3++] = E->C12[part];
E->PINS[neighbor][k3++] = E->CElement[part];
+ for(k=0;k < E->tracers_add_flavors;k++)
+ E->PINS[neighbor][k3++] = E->tflavors[part*E->tracers_add_flavors+k];
+
}
//if((E->parallel.me==0) && (E->monitor.solution_cycles>199))fprintf(stderr,"ta inner loop ok\n");
}
Modified: mc/3D/CitcomCU/trunk/src/Convection.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Convection.c 2011-04-26 19:02:48 UTC (rev 18296)
+++ mc/3D/CitcomCU/trunk/src/Convection.c 2011-04-27 12:19:44 UTC (rev 18297)
@@ -52,7 +52,8 @@
{
input_int("composition", &(E->control.composition), "0", E->parallel.me);
-
+ input_int("tracers_add_flavors", &(E->tracers_add_flavors), "0", E->parallel.me);
+
if(E->control.composition)
E->next_buoyancy_field = PG_timestep_particle;
else
@@ -183,9 +184,13 @@
E->Vpred[i] = (float *)malloc((E->advection.markers_uplimit + 1) * sizeof(float));
E->XMCpred[i] = (double *)malloc((E->advection.markers_uplimit + 1) * sizeof(double));
E->XMC[i] = (double *)malloc((E->advection.markers_uplimit + 1) * sizeof(double));
- E->C12 = (int *)malloc((E->advection.markers_uplimit + 1) * sizeof(int));
- E->traces_leave = (int *)malloc((E->advection.markers_uplimit + 1) * sizeof(int));
- E->CElement = (int *)malloc((E->advection.markers_uplimit + 1) * sizeof(int));
+ if(i==1){ /* those should only get allocated once */
+ E->C12 = (int *)malloc((E->advection.markers_uplimit + 1) * sizeof(int));
+ E->traces_leave = (int *)malloc((E->advection.markers_uplimit + 1) * sizeof(int));
+ E->CElement = (int *)malloc((E->advection.markers_uplimit + 1) * sizeof(int));
+ if(E->tracers_add_flavors)
+ E->tflavors = (int *)(int *)malloc((E->advection.markers_uplimit + 1) * E->tracers_add_flavors * sizeof(int));
+ }
}
}
Modified: mc/3D/CitcomCU/trunk/src/Drive_solvers.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Drive_solvers.c 2011-04-26 19:02:48 UTC (rev 18296)
+++ mc/3D/CitcomCU/trunk/src/Drive_solvers.c 2011-04-27 12:19:44 UTC (rev 18297)
@@ -142,7 +142,7 @@
/* end for stress type iterations */
} while(iterate &&
(dUdot_mag > E->viscosity.sdepv_misfit) &&
- (E->monitor.visc_iter_count < 50) );
+ (E->monitor.visc_iter_count < E->monitor.max_sdep_visc_iter) );
if(iterate){ /* free the delta_U array */
free((void *)delta_U);
}
Modified: mc/3D/CitcomCU/trunk/src/Parallel_related.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Parallel_related.c 2011-04-26 19:02:48 UTC (rev 18296)
+++ mc/3D/CitcomCU/trunk/src/Parallel_related.c 2011-04-27 12:19:44 UTC (rev 18297)
@@ -1008,11 +1008,13 @@
{
//int target_proc, kk, e, node, i, ii, j, k, bound, type, idb, msginfo[8];
int target_proc, kk, k, idb;
+ int rioff;
-
MPI_Status status[100];
MPI_Request request[100];
+ rioff = 2 + E->tracers_add_flavors;
+
idb = 0;
for(k = 1; k <= E->parallel.no_neighbors; k++)
{
@@ -1020,7 +1022,7 @@
{
target_proc = E->parallel.neighbors[k];
idb++;
- kk = E->parallel.traces_transfer_number[k] * 2 + 1;
+ kk = E->parallel.traces_transfer_number[k] * rioff + 1;
MPI_Isend(E->PINS[k], kk, MPI_INT, target_proc, 1, MPI_COMM_WORLD, &request[idb - 1]);
idb++;
kk = E->parallel.traces_transfer_number[k] * 2 * E->mesh.nsd + 1;
@@ -1036,7 +1038,7 @@
{
target_proc = E->parallel.neighbors[k];
idb++;
- kk = E->parallel.traces_receive_number[k] * 2 + 1;
+ kk = E->parallel.traces_receive_number[k] * rioff + 1;
MPI_Irecv(E->RINS[k], kk, MPI_INT, target_proc, 1, MPI_COMM_WORLD, &request[idb - 1]);
idb++;
kk = E->parallel.traces_receive_number[k] * 2 * E->mesh.nsd + 1;
Modified: mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2011-04-26 19:02:48 UTC (rev 18296)
+++ mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2011-04-27 12:19:44 UTC (rev 18297)
@@ -153,6 +153,7 @@
input_boolean("sdepv_start_from_newtonian", &(E->viscosity.sdepv_start_from_newtonian), "off",m);
+ input_int("max_sdep_visc_iter",&(E->monitor.max_sdep_visc_iter),"50",m); /* max number of powerlaw iterations */
input_boolean("BDEPV",&(E->viscosity.BDEPV),"off",m);
input_boolean("CDEPV",&(E->viscosity.CDEPV),"off",m);
@@ -1443,6 +1444,9 @@
override all other viscosities (this is exact same as below but
for the viscosity assignment (repeated here to avoid if statments)
+
+ will make compositional viscosity an on/off kind of thing
+
*/
for(i = 1; i <= nel; i++){
for(kk = 1; kk <= ends; kk++){
@@ -1455,8 +1459,8 @@
cc_loc = 0.0;
for(kk = 1; kk <= ends; kk++)
cc_loc += CC[kk] * E->N.vpt[GNVINDEX(kk, jj)];
- vmean = exp(cc_loc * logv[1] + (1.0-cc_loc) * logv[0]);
- EEta[ (i-1)*vpts + jj ] = vmean;
+ /* mean between background viscosity and second material viscosity */
+ EEta[ (i-1)*vpts + jj ]= exp(cc_loc * logv[1] + (1.0-cc_loc) * log(EEta[ (i-1)*vpts + jj ]));
} /* end jj loop */
}
}else{
@@ -1481,7 +1485,7 @@
for(kk = 1; kk <= ends; kk++){/* the vpt takes care of averaging */
cc_loc += CC[kk] * E->N.vpt[GNVINDEX(kk, jj)];
}
- /* geometric mean of viscosity */
+ /* geometric mean of compositional viscosity prefactors */
vmean = exp(cc_loc * logv[1] + (1.0-cc_loc) * logv[0]);
EEta[ (i-1)*vpts + jj ] *= vmean;
} /* end jj loop */
Modified: mc/3D/CitcomCU/trunk/src/global_defs.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/global_defs.h 2011-04-26 19:02:48 UTC (rev 18296)
+++ mc/3D/CitcomCU/trunk/src/global_defs.h 2011-04-27 12:19:44 UTC (rev 18297)
@@ -650,7 +650,7 @@
int solution_cycles,solution_cycles_out;
int visc_iter_count;
-
+ int max_sdep_visc_iter;
float time_scale,time_scale_ma;
float length_scale;
float viscosity_scale;
@@ -941,11 +941,14 @@
float *Vpred[4];
double *XMC[4];
double *XMCpred[4];
- int *C12;
+ int *C12; /* this shold be unsigned short */
int *traces_leave;
int *traces_leave_index;
int *CElement;
+ int *tflavors; /* this should also be unsigned short */
+ int tracers_add_flavors;
+
int *RG[4];
double *XRG[4];
double *XP[4], XG1[4], XG2[4];
More information about the CIG-COMMITS
mailing list