[cig-commits] r18486 - mc/3D/CitcomCU/trunk/src
becker at geodynamics.org
becker at geodynamics.org
Sat May 28 20:13:35 PDT 2011
Author: becker
Date: 2011-05-28 20:13:35 -0700 (Sat, 28 May 2011)
New Revision: 18486
Modified:
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/Output_gzdir.c
mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.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
mc/3D/CitcomCU/trunk/src/prototypes.h
mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h
Log:
Added tracer flavors to determine, say, material properties. For now,
only USE_GGRD initialization of a single flavor is supported which allows
selection of Byerlee or not.
Modified: mc/3D/CitcomCU/trunk/src/Composition_adv.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Composition_adv.c 2011-05-29 00:17:47 UTC (rev 18485)
+++ mc/3D/CitcomCU/trunk/src/Composition_adv.c 2011-05-29 03:13:35 UTC (rev 18486)
@@ -93,8 +93,13 @@
transfer_markers_processors(E, on_off);
+ /* assign element numbers to markers */
element_markers(E, on_off);
+
+ /* get nodal values from markers */
get_C_from_markers(E, C);
+ if(E->tracers_add_flavors)
+ get_CF_from_markers(E,E->CF);
E->advection.markers_g = -1;
if(E->advection.timesteps % 10 == 0)
@@ -136,7 +141,12 @@
transfer_markers_processors(E, on_off);
/* predicted compositional field at t+dt */
element_markers(E, on_off);
+ /* update nodal values */
get_C_from_markers(E, C);
+ if(E->tracers_add_flavors)
+ get_CF_from_markers(E,E->CF);
+
+
return;
}
@@ -318,7 +328,7 @@
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->tflavors[ii][k] = E->RINS[neighbor][j * rioff + 2 + k];
E->traces_leave[ii] = 0;
}
}
@@ -347,7 +357,7 @@
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->tflavors[ii][k] = E->RINS[neighbor][j * rioff + 2 + k];
E->traces_leave[ii] = 0;
}
}
@@ -380,8 +390,7 @@
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->tflavors[ii][k] = E->tflavors[i][k];
E->traces_leave[ii] = 0;
}
}
@@ -431,7 +440,7 @@
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];
+ E->PINS[neighbor][k3++] = E->tflavors[part][k];
}
//if((E->parallel.me==0) && (E->monitor.solution_cycles>199))fprintf(stderr,"ta inner loop ok\n");
@@ -545,7 +554,105 @@
return;
}
+/*
+obtain nodal flavor from markers
+
+*/
+void get_CF_from_markers(struct All_variables *E, int **CF)
+{
+ int el, itracer, imark, jnode, node,k,i,j;
+ float *dmin,dist;
+ static int method = 0,been_here = 0;
+ float temp3, temp1, temp0;
+ static int *element[3];
+
+ const int nno = E->lmesh.nno;
+ const int nel = E->lmesh.nel;
+ const int dims = E->mesh.nsd;
+ const int ends = enodes[dims];
+ const int lev = E->mesh.levmax;
+ const int nno1 = nno + 1;
+
+ if(!been_here){
+ /* decide on method */
+ for(i=-1,k=0;k < E->tracers_add_flavors;k++)
+ if(E->tmaxflavor[k] > i)
+ i = E->tmaxflavor[k];
+ if(i > 1) /* any of the flavors is more than one */
+ method = 0;
+ else{
+ method = 1; /* for counting entries */
+ element[0] = (int *)malloc((nel + 1) * sizeof(int));
+ element[1] = (int *)malloc((nel + 1) * sizeof(int));
+ }
+ if(E->parallel.me==0)
+ if(method == 0)
+ fprintf(stderr,"get_CF_from_markers: using closest tracer method\n");
+ else
+ fprintf(stderr,"get_CF_from_markers: using mean flavor method for range 0..1\n");
+ been_here = 1;
+ }
+ dmin = (float *)malloc((nno1) * sizeof(float));
+ switch(method){
+ case 0:
+ /* determine nodal flavor based on the closest tracer */
+ for(k=0;k < nno1;k++)
+ dmin[k] = 1e10;
+ for(itracer=0;itracer < E->advection.markers;itracer++){
+ el = E->CElement[itracer] ;
+ for(jnode = 1; jnode <= ends; jnode++){ /* loop through the nodes within this element */
+ node = E->ien[el].node[jnode];
+ dist = (float) distance_to_node(E->XMC[1][itracer],E->XMC[2][itracer],E->XMC[3][itracer],E,node);
+ if(dmin[node] > dist){
+ dmin[node] = dist;
+ for(k=0;k < E->tracers_add_flavors;k++)
+ CF[k][node] = E->tflavors[itracer][k];
+ }
+ }
+ }
+ break;
+ case 1:
+ /*
+ mean flavor method : compute based on the mean flavor of each element
+ */
+ for(k=0;k < E->tracers_add_flavors;k++){
+ for(el = 1; el <= nel; el++){
+ element[0][el] = element[1][el] = 0;
+ }
+ for(i = 1; i <= nno; i++){
+ dmin[i] = 0.0;
+ }
+ for(imark = 1; imark <= E->advection.markers; imark++){
+ element[E->tflavors[imark][k]][E->CElement[imark]]++;
+ }
+ for(el = 1; el <= nel; el++){
+ temp0 = element[0][el];
+ temp1 = element[1][el];
+ if(element[0][el] || element[1][el])
+ temp3 = temp1 / (temp0 + temp1); /* mean elemental CF */
+ else
+ temp3 = 0; /* no tracers */
+ for(j = 1; j <= ends; j++){
+ node = E->ien[el].node[j];
+ dmin[node] += E->TWW[lev][el].node[j] * temp3;
+ }
+ }
+ exchange_node_f20(E, dmin, E->mesh.levmax);
+ for(node = 1; node <= nno; node++){
+ CF[k][node] = (int) (dmin[node] * E->Mass[node] +.5);
+ }
+ }
+ break; /* */
+ default:
+ myerror("method for flavor interpolation undefined",E);
+ break;
+ }
+ free(dmin);
+}
+
+
+
/* ================================================ */
void element_markers(struct All_variables *E, int con)
{
Modified: mc/3D/CitcomCU/trunk/src/Convection.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Convection.c 2011-05-29 00:17:47 UTC (rev 18485)
+++ mc/3D/CitcomCU/trunk/src/Convection.c 2011-05-29 03:13:35 UTC (rev 18486)
@@ -170,7 +170,7 @@
void convection_initial_fields(struct All_variables *E)
{
- int i;
+ int i,j;
if(E->control.composition)
{
@@ -188,8 +188,18 @@
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));
+ if(E->tracers_add_flavors){
+ /* nodal compotional flavors */
+ E->CF = (int **)malloc( E->tracers_add_flavors * sizeof(int *));
+ /* maximum flavor value */
+ E->tmaxflavor = (int *)calloc(E->tracers_add_flavors,sizeof(int ));
+ for(j=0;j<E->tracers_add_flavors;j++)
+ E->CF[j] = (int *)calloc((E->lmesh.nno + 1), sizeof(int));
+ /* tracer flavor value */
+ E->tflavors = (int **)malloc((E->advection.markers_uplimit + 1) * sizeof(int *));
+ for(j = 0; j < E->advection.markers_uplimit + 1;j++)
+ E->tflavors[j] = (int *)calloc(E->tracers_add_flavors,sizeof(int));
+ }
}
}
}
@@ -220,9 +230,9 @@
NOTE:
- THERE'S A DIFFERENT ROUTINE FOR GGRD HANDLING IN GGRD_HANDLING
- WHICH HAS DIFFERNET LOGIC FOR MARKER INIT ETC. THE OTHER ONE GETS
- CALLED when USE_GGRD is compiled in
+ there's a different routine for ggrd handling in ggrd_handling
+ which has differnet logic for marker init etc. the other one gets
+ called when USE_GGRD is compiled in
=============================== */
@@ -247,6 +257,9 @@
p = 0;
+ if(E->tracers_add_flavors)
+ myerror("convection_initial_temperature: not set up for tracer flavors",E);
+
if(E->control.restart == 0)
{
if(E->control.CART3D)
@@ -596,8 +609,10 @@
*/
}
-
+ /* nodal values */
get_C_from_markers(E, E->C);
+ if(E->tracers_add_flavors)
+ get_CF_from_markers(E,E->CF);
return;
}
@@ -605,7 +620,7 @@
void convection_initial_markers(struct All_variables *E,int use_element_nodes_for_init_c)
{
//int el, i, j, k, p, node, ii, jj;
- int el, node,j;
+ int el, ntracer,j,i;
//double x, y, z, r, t, f, dX[4], dx, dr;
double x, y, z, r, t, f, dX[4];
//char input_s[100], output_file[255];
@@ -617,7 +632,7 @@
if(E->control.CART3D)
{
- node = 0;
+ ntracer = 0;
do
{
x = drand48() * (E->XG2[1] - E->XG1[1]);
@@ -626,33 +641,33 @@
if((x >= E->XP[1][1] && x <= E->XP[1][E->lmesh.nox]) && (y >= E->XP[2][1] && y <= E->XP[2][E->lmesh.noy]) && (z >= E->XP[3][1] && z <= E->XP[3][E->lmesh.noz]))
{
- node++;
- E->XMC[1][node] = x;
- E->XMC[2][node] = y;
- E->XMC[3][node] = z;
+ ntracer++;
+ E->XMC[1][ntracer] = x;
+ E->XMC[2][ntracer] = y;
+ E->XMC[3][ntracer] = z;
- el = get_element(E, E->XMC[1][node], E->XMC[2][node], E->XMC[3][node], dX);
- E->CElement[node] = el;
+ el = get_element(E, E->XMC[1][ntracer], E->XMC[2][ntracer], E->XMC[3][ntracer], dX);
+ E->CElement[ntracer] = el;
if(use_element_nodes_for_init_c){
for(temp=0.0,j = 1; j <= ends; j++)
temp += E->C[E->ien[el].node[j]];
temp /= ends;
if(temp >0.5)
- E->C12[node] = 1;
+ E->C12[ntracer] = 1;
else
- E->C12[node] = 0;
+ E->C12[ntracer] = 0;
}else{ /* use depth */
- if(E->XMC[3][node] > E->viscosity.zcomp)
- E->C12[node] = 0;
+ if(E->XMC[3][ntracer] > E->viscosity.zcomp)
+ E->C12[ntracer] = 0;
else
- E->C12[node] = 1;
+ E->C12[ntracer] = 1;
}
}
- } while(node < E->advection.markers);
+ } while(ntracer < E->advection.markers);
}
else if(E->control.Rsphere)
{
- node = 0;
+ ntracer = 0;
do
{
x = (drand48() - 0.5) * 2.0;
@@ -664,34 +679,82 @@
f = myatan(y, x);
if((t >= E->XP[1][1] && t <= E->XP[1][E->lmesh.nox]) && (f >= E->XP[2][1] && f <= E->XP[2][E->lmesh.noy]) && (r >= E->XP[3][1] && r <= E->XP[3][E->lmesh.noz]))
{
- node++;
- E->XMC[1][node] = t;
- E->XMC[2][node] = f;
- E->XMC[3][node] = r;
- el = get_element(E, E->XMC[1][node], E->XMC[2][node], E->XMC[3][node], dX);
- E->CElement[node] = el;
+ ntracer++;
+ E->XMC[1][ntracer] = t;
+ E->XMC[2][ntracer] = f;
+ E->XMC[3][ntracer] = r;
+ el = get_element(E, E->XMC[1][ntracer], E->XMC[2][ntracer], E->XMC[3][ntracer], dX);
+ E->CElement[ntracer] = el;
if(use_element_nodes_for_init_c){
for(temp=0.0,j = 1; j <= ends; j++)
temp += E->C[E->ien[el].node[j]];
temp /= ends;
if(temp >0.5)
- E->C12[node] = 1;
+ E->C12[ntracer] = 1;
else
- E->C12[node] = 0;
+ E->C12[ntracer] = 0;
}else{ /* use depth */
if(r > E->viscosity.zcomp)
- E->C12[node] = 0;
+ E->C12[ntracer] = 0;
else
- E->C12[node] = 1;
+ E->C12[ntracer] = 1;
}
+
+
}
- } while(node < E->advection.markers);
+ } while(ntracer < E->advection.markers);
}
+
+ /* assign tracers based on ggrd */
+ if(E->tracers_add_flavors){
+#ifdef USE_GGRD
+ assign_flavor_to_tracer_from_grd(E);
+#else
+ myerror("convection_initial_markers: flavor init requires GGRD compilation",E);
+#endif
+ }
+ /*
+ get nodal values
+ */
+ get_C_from_markers(E, E->C);
+ if(E->tracers_add_flavors)
+ get_CF_from_markers(E,E->CF);
+ return;
+}
+/*
+assign a flavor to this tracer based on proximity to a node
- get_C_from_markers(E, E->C);
+*/
+void assign_flavor_to_tracer_based_on_node(struct All_variables *E, int ntracer, int use_element_nodes_for_init_c)
+{
+ int j,el,i;
+ const int dims = E->mesh.nsd;
+ const int ends = enodes[dims];
+ float dmin,dist;
+ int nmin,lnode;
- return;
+ el = E->CElement[ntracer]; /* element we're in */
+ if(use_element_nodes_for_init_c){
+ for(dmin = 1e20,j = 1; j <= ends; j++){
+ lnode = E->ien[el].node[j];
+ dist = distance_to_node(E->XMC[1][ntracer],E->XMC[2][ntracer],E->XMC[3][ntracer],E,lnode);
+ if(dist < dmin){ /* pick the node that is closest to the tracer */
+ dmin = dist;
+ nmin = lnode;
+ }
+ }
+ for(i=0;i < E->tracers_add_flavors;i++)
+ E->tflavors[ntracer][i] = E->CF[i][lnode];
+ }else{
+ for(i=0;i < E->tracers_add_flavors;i++){
+ if(E->XMC[3][ntracer] > E->viscosity.zcomp)
+ E->tflavors[ntracer][i] = 0;
+ else
+ E->tflavors[ntracer][i] = 1;
+ }
+ }
+
}
void setup_plume_problem(struct All_variables *E)
Modified: mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2011-05-29 00:17:47 UTC (rev 18485)
+++ mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2011-05-29 03:13:35 UTC (rev 18486)
@@ -137,7 +137,7 @@
FALSE,FALSE))
myerror("grdtrack x-z init error",E);
}else{ /* several slab slices */
- for(slice=0;slice<E->control.ggrd_slab_slice;slice++){
+ for(slice=0;slice < E->control.ggrd_slab_slice;slice++){
sprintf(pfile,"%s.%i.grd",E->control.ggrd.temp.gfile,slice+1);
if(ggrd_grdtrack_init_general(FALSE,pfile,
char_dummy,"",(E->control.ggrd_ss_grd+slice),
@@ -194,7 +194,7 @@
for(slice=hit=0;(!hit) && (slice < E->control.ggrd_slab_slice);slice++){
if(E->control.Rsphere) {
if(in_slab_slice(E->SX[1][node],slice,E)){
- /* spherical interpolation */
+ /* spherical interpolation for a slice ! */
ggrd_grdtrack_interpolate_xy((double)E->SX[2][node] * ONEEIGHTYOVERPI,
(double)E->SX[1][node],
(E->control.ggrd_ss_grd+slice),&tadd,
@@ -376,97 +376,16 @@
if(E->control.ggrd.use_comp){ /* read composition init from grid */
if(!E->control.composition)
myerror("comp grd init but no composition control set!?",E);
- if(E->parallel.me==0)
- fprintf(stderr,"convection_initial_temperature: using GMT grd files for composition\n");
- if(E->parallel.me > 0){
- mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), mpi_tag,
- MPI_COMM_WORLD, &mpi_stat);
- }
- if(E->control.ggrd_slab_slice){
- if(E->control.ggrd_slab_slice == 1){
- if(ggrd_grdtrack_init_general(FALSE,E->control.ggrd.comp.gfile,
- char_dummy,"",E->control.ggrd_ss_grd,
- (E->parallel.me==0),
- FALSE,FALSE))
- myerror("grdtrack x-z init error",E);
- }else{ /* several slab slices */
- for(slice=0;slice<E->control.ggrd_slab_slice;slice++){
- sprintf(pfile,"%s.%i.grd",E->control.ggrd.comp.gfile,slice+1);
- if(ggrd_grdtrack_init_general(FALSE,pfile,
- char_dummy,"",(E->control.ggrd_ss_grd+slice),
- (E->parallel.me==0),
- FALSE,FALSE))
- myerror("grdtrack x-z init error",E);
- }
- }
- }else{
- /* 3-D */
- if(ggrd_grdtrack_init_general(TRUE,E->control.ggrd.comp.gfile,
- E->control.ggrd.comp.dfile,"",
- E->control.ggrd.comp.d,
- /* (E->parallel.me==0)*/
- FALSE,FALSE,FALSE))
- myerror("grdtrack init error",E);
- }
- if(E->parallel.me < E->parallel.nproc-1){
- /*
- tell the next processor to go ahead
- */
- mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT,
- (E->parallel.me+1), mpi_tag, MPI_COMM_WORLD);
- }else{
- fprintf(stderr,"convection_initial_temperature: last processor done comp with grd init\n");
- }
- for(i=1;i<=noy;i++)
- for(j=1;j<=nox;j++)
- for(k=1;k<=noz;k++) {
- node=k+(j-1)*noz+(i-1)*nox*noz; /* offset */
- if(E->control.ggrd_slab_slice){
- /* slab */
- for(hit = slice=0;(!hit) && (slice < E->control.ggrd_slab_slice);slice++){
- if(E->control.Rsphere) {
- if(in_slab_slice(E->SX[1][node],slice,E)){
- /* spherical interpolation */
- ggrd_grdtrack_interpolate_xy((double)E->SX[2][node] * ONEEIGHTYOVERPI,
- (double)E->SX[1][node],
- (E->control.ggrd_ss_grd+slice),&tadd,
- FALSE);
- hit = 1;
- }
- }else{ /* cartesian interpolation */
- if(in_slab_slice(E->X[2][node],slice,E)){
- ggrd_grdtrack_interpolate_xy((double)E->X[1][node],
- (double)E->X[3][node],
- (E->control.ggrd_ss_grd+slice),&tadd,
- FALSE);
- hit = 1 ;
- }
- }
- } /* end slice loop */
- if(!hit)
- tadd = 0;
- /* end slab slice */
- }else{
- /* 3-D */
- if(E->control.Rsphere) /* spherical interpolation */
- ggrd_grdtrack_interpolate_rtp((double)E->SX[3][node],
- (double)E->SX[1][node],
- (double)E->SX[2][node],
- E->control.ggrd.comp.d,&tadd,
- FALSE,FALSE);
- else /* cartesian interpolation */
- ggrd_grdtrack_interpolate_xyz((double)E->X[1][node],
- (double)E->X[2][node],
- (double)E->X[3][node],
- E->control.ggrd.comp.d,&tadd,
- FALSE);
- }
- E->C[node] = E->control.ggrd.comp.offset + tadd *
- E->control.ggrd.comp.scale;
- }
- /* free the structure, not needed anymore */
- ggrd_grdtrack_free_gstruc(E->control.ggrd.comp.d);
+ ggrd_deal_with_composition_input(E, 1); /* compositional (density) */
+
+ /* flavors are dealt with based on markers, and then assigned
+ to nodes */
+
+ //if(E->tracers_add_flavors) /* flavors */
+ //ggrd_deal_with_composition_input(E, 0);
+
+
} /* end grid cinit */
/*
@@ -511,6 +430,257 @@
return;
}
+
+/*
+
+assign composition or flavors to nodes from grd files
+
+
+*/
+void ggrd_deal_with_composition_input(struct All_variables *E, int assign_composition)
+{
+ MPI_Status mpi_stat;
+ int mpi_rc, mpi_tag=1;
+ int mpi_inmsg, mpi_success_message = 1;
+ char ingfile[1000],indfile[1000],pfile[1000];
+ int use_nearneighbor,i,j,nox, noy, noz,slice,node,hit,k;
+ char *char_dummy="";
+ double tadd;
+ int *tmaxflavor;
+
+
+ noy = E->lmesh.noy;
+ noz = E->lmesh.noz;
+ nox = E->lmesh.nox;
+
+ if(E->parallel.me == 0){
+ if(assign_composition)
+ fprintf(stderr,"ggrd_deal_with_composition_input: using GMT grd files for composition\n");
+ else{
+ fprintf(stderr,"grd_deal_with_composition_input: using GMT grd files for flavors\n");
+
+ }
+ }
+ /* filenames */
+ if(assign_composition){
+ /* composition */
+ use_nearneighbor = FALSE;
+ strncpy(ingfile,E->control.ggrd.comp.gfile,1000);
+ strncpy(indfile,E->control.ggrd.comp.dfile,1000);
+ }else{
+ /* flavor */
+ use_nearneighbor = TRUE;
+ tmaxflavor = (int *)calloc(E->tracers_add_flavors,sizeof(int));
+ strncpy(ingfile,E->control.ggrd_flavor_gfile,1000);
+ strncpy(indfile,E->control.ggrd_flavor_dfile,1000);
+ if(E->tracers_add_flavors != 1)
+ myerror("ggrd_deal_with_composition_input: only set up for one flavor",E);
+
+ }
+ if(E->parallel.me > 0){
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), mpi_tag,
+ MPI_COMM_WORLD, &mpi_stat);
+ }
+ if(E->control.ggrd_slab_slice){
+ if(E->control.ggrd_slab_slice == 1){
+ if(ggrd_grdtrack_init_general(FALSE,ingfile, char_dummy,"",E->control.ggrd_ss_grd,(E->parallel.me==0),
+ FALSE,use_nearneighbor))
+ myerror("grdtrack x-z init error",E);
+ }else{ /* several slab slices */
+ for(slice=0;slice < E->control.ggrd_slab_slice;slice++){
+ sprintf(pfile,"%s.%i.grd",ingfile,slice+1);
+ if(ggrd_grdtrack_init_general(FALSE,pfile, char_dummy,"",(E->control.ggrd_ss_grd+slice),(E->parallel.me==0),
+ FALSE,use_nearneighbor))
+ myerror("grdtrack x-z init error",E);
+ }
+ }
+ }else{
+ /* 3-D */
+ if(ggrd_grdtrack_init_general(TRUE,ingfile,indfile,"",E->control.ggrd.comp.d,FALSE,FALSE,use_nearneighbor))
+ myerror("grdtrack init error",E);
+ }
+ if(E->parallel.me < E->parallel.nproc-1){
+ /*
+ tell the next processor to go ahead
+ */
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT,
+ (E->parallel.me+1), mpi_tag, MPI_COMM_WORLD);
+ }else{
+ fprintf(stderr,"ggrd_deal_with_composition_input: last processor done comp with grd init\n");
+ }
+ for(i=1;i<=noy;i++) {
+ for(j=1;j<=nox;j++){
+ for(k=1;k<=noz;k++) {
+ node=k+(j-1)*noz+(i-1)*nox*noz; /* offset */
+ if(E->control.ggrd_slab_slice){
+ /* slab */
+ for(hit = slice=0;(!hit) && (slice < E->control.ggrd_slab_slice);slice++){
+ if(E->control.Rsphere) {
+ if(in_slab_slice(E->SX[1][node],slice,E)){
+ /* spherical interpolation for a slice */
+ ggrd_grdtrack_interpolate_xy((double)E->SX[2][node] * ONEEIGHTYOVERPI,
+ (double)E->SX[1][node],
+ (E->control.ggrd_ss_grd+slice),&tadd,
+ FALSE);
+ hit = 1;
+ }
+ }else{ /* cartesian interpolation */
+ if(in_slab_slice(E->X[2][node],slice,E)){
+
+ ggrd_grdtrack_interpolate_xy((double)E->X[1][node],(double)E->X[3][node],
+ (E->control.ggrd_ss_grd+slice),&tadd,
+ FALSE);
+ hit = 1 ;
+ }
+ }
+ } /* end slice loop */
+ if(!hit)
+ tadd = 0;
+ /* end slab slice */
+ }else{
+ /* 3-D */
+ if(E->control.Rsphere) /* spherical interpolation */
+ ggrd_grdtrack_interpolate_rtp((double)E->SX[3][node],(double)E->SX[1][node],(double)E->SX[2][node],
+ E->control.ggrd.comp.d,&tadd,
+ FALSE,FALSE);
+ else /* cartesian interpolation */
+ ggrd_grdtrack_interpolate_xyz((double)E->X[1][node],(double)E->X[2][node],(double)E->X[3][node],
+ E->control.ggrd.comp.d,&tadd,
+ FALSE);
+ }
+ if(assign_composition){
+ E->C[node] = E->control.ggrd.comp.offset + tadd * E->control.ggrd.comp.scale;
+ }else{ /* flavor */
+
+ E->CF[0][node] = (int)(tadd+.5);
+ if(E->CF[0][node] < 0)
+ E->CF[0][node]=0;
+ if(tmaxflavor[0] < E->CF[0][node]) /* maximum flavor count */
+ tmaxflavor[0] = E->CF[0][node];
+ }
+ }
+ }
+ }
+ if(!assign_composition){
+ MPI_Allreduce(tmaxflavor,E->tmaxflavor,E->tracers_add_flavors,
+ MPI_INT,MPI_MAX,MPI_COMM_WORLD);
+ if(E->parallel.me == 0)
+ fprintf(stderr,"ggrd_deal_with_composition_input: assigned maximum flavor %i\n",
+ E->tmaxflavor[0]);
+ free(tmaxflavor);
+ }
+
+
+ /* free the structure, not needed anymore */
+ ggrd_grdtrack_free_gstruc(E->control.ggrd.comp.d);
+}
+
+
+/*
+
+assign composition or flavors to nodes from grd files
+
+
+*/
+void assign_flavor_to_tracer_from_grd(struct All_variables *E)
+{
+ MPI_Status mpi_stat;
+ int mpi_rc, mpi_tag=1;
+ int mpi_inmsg, mpi_success_message = 1;
+ char pfile[1000];
+ int use_nearneighbor = TRUE,i,j,slice,hit,k;
+ char *char_dummy="";
+ double tadd;
+ struct ggrd_gt grd[4];
+ int *tmaxflavor;
+ tmaxflavor = (int *)calloc(E->tracers_add_flavors,sizeof(int));
+ if(E->parallel.me == 0)
+ fprintf(stderr,"assign_flavor_to_tracer_from_grd: using GMT grd files for flavors\n");
+
+ if(E->parallel.me > 0){
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), mpi_tag,
+ MPI_COMM_WORLD, &mpi_stat);
+ }
+ if(E->control.ggrd_slab_slice){
+ if(E->control.ggrd_slab_slice == 1){
+ if(ggrd_grdtrack_init_general(FALSE,E->control.ggrd_flavor_gfile, char_dummy,"",
+ grd,(E->parallel.me==0),FALSE,use_nearneighbor))
+ myerror("grdtrack x-z init error",E);
+ }else{ /* several slab slices */
+ for(slice=0;slice < E->control.ggrd_slab_slice;slice++){
+ sprintf(pfile,"%s.%i.grd",E->control.ggrd_flavor_gfile,slice+1);
+ if(ggrd_grdtrack_init_general(FALSE,pfile, char_dummy,"",(grd+slice),(E->parallel.me==0),
+ FALSE,use_nearneighbor))
+ myerror("grdtrack x-z init error",E);
+ }
+ }
+ }else{
+ /* 3-D */
+ if(ggrd_grdtrack_init_general(TRUE,E->control.ggrd_flavor_gfile,E->control.ggrd_flavor_dfile,"",
+ grd,FALSE,FALSE,use_nearneighbor))
+ myerror("grdtrack init error",E);
+ }
+ if(E->parallel.me < E->parallel.nproc-1){
+ /*
+ tell the next processor to go ahead
+ */
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT,
+ (E->parallel.me+1), mpi_tag, MPI_COMM_WORLD);
+ }else{
+ fprintf(stderr,"ggrd_deal_with_composition_input: last processor done comp with grd init\n");
+ }
+ for(i=0;i < E->advection.markers;i++) {
+ if(E->control.ggrd_slab_slice){
+ /* slab */
+ for(hit = slice=0;(!hit) && (slice < E->control.ggrd_slab_slice);slice++){
+ if(E->control.Rsphere) {
+ if(in_slab_slice(E->XMC[1][i],slice,E)){
+ /* spherical interpolation for a slice */
+ ggrd_grdtrack_interpolate_xy((double)E->XMC[2][i]* ONEEIGHTYOVERPI,(double)E->XMC[1][i],
+ (grd+slice),&tadd,
+ FALSE);
+ hit = 1;
+ }
+ }else{ /* cartesian interpolation */
+ if(in_slab_slice(E->XMC[2][i],slice,E)){
+ ggrd_grdtrack_interpolate_xy((double)E->XMC[1][i],(double)E->XMC[3][i],
+ (grd+slice),&tadd,FALSE);
+ hit = 1 ;
+ }
+ }
+ } /* end slice loop */
+ if(!hit)
+ tadd = 0;
+ /* end slab slice */
+ }else{
+ /* 3-D */
+ if(E->control.Rsphere) /* spherical interpolation */
+ ggrd_grdtrack_interpolate_rtp((double)E->XMC[3][i],(double)E->XMC[1][i],(double)E->XMC[2][i],
+ grd,&tadd,FALSE,FALSE);
+ else /* cartesian interpolation */
+ ggrd_grdtrack_interpolate_xyz((double)E->XMC[1][i],(double)E->XMC[2][i],(double)E->XMC[3][i],
+ grd,&tadd,FALSE);
+ }
+
+ /* assign to tracer */
+ E->tflavors[i][0] = (int)(tadd+.5);
+ if(tmaxflavor[0] < E->tflavors[i][0]) /* maximum flavor count */
+ tmaxflavor[0] = E->tflavors[i][0];
+ } /* end tracer loop */
+
+ MPI_Allreduce(tmaxflavor,E->tmaxflavor,E->tracers_add_flavors,
+ MPI_INT,MPI_MAX,MPI_COMM_WORLD);
+ if(E->parallel.me == 0){
+ fprintf(stderr,"ggrd_deal_with_composition_input: assigned maximum flavor %i\n",
+ E->tmaxflavor[0]);
+
+ }
+ free(tmaxflavor);
+
+}
+
+
+
int in_slab_slice(float coord, int slice, struct All_variables *E)
{
if((slice < 0)||(slice > E->control.ggrd_slab_slice-1))
Modified: mc/3D/CitcomCU/trunk/src/Instructions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Instructions.c 2011-05-29 00:17:47 UTC (rev 18485)
+++ mc/3D/CitcomCU/trunk/src/Instructions.c 2011-05-29 03:13:35 UTC (rev 18486)
@@ -848,13 +848,19 @@
/* comp */
input_boolean("ggrd_cinit",&(E->control.ggrd.use_comp),"off", m);
input_double("ggrd_cinit_scale",&(E->control.ggrd.comp.scale),"1.0", m);
+
input_string("ggrd_cinit_gfile",E->control.ggrd.comp.gfile,"", m);
input_string("ggrd_cinit_dfile",E->control.ggrd.comp.dfile,"", m);
+
+ input_string("ggrd_cflavor_init_gfile",E->control.ggrd_flavor_gfile,"", m);
+ input_string("ggrd_cflavor_init_dfile",E->control.ggrd_flavor_dfile,"", m);
+
+
input_double("ggrd_cinit_offset",&(E->control.ggrd.comp.offset),"0.0", m);
/* slab slice handling */
input_int("slab_slice",&(E->control.ggrd_slab_slice),"0", m);
if(E->control.ggrd_slab_slice > 3)
- myerror("too many slab slices",E);
+ myerror("too many slab slices, only four allowed",E);
input_float_vector("slab_theta_bound",E->control.ggrd_slab_slice,(E->control.ggrd_slab_theta_bound), m);
/*
@@ -1038,8 +1044,8 @@
E->monitor.tau_scale = (E->data.ref_viscosity/E->monitor.time_scale); /* scaling stress */
+ (E->problem_settings) (E);
- (E->problem_settings) (E);
viscosity_parameters(E);
return;
}
Modified: mc/3D/CitcomCU/trunk/src/Output_gzdir.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2011-05-29 00:17:47 UTC (rev 18485)
+++ mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2011-05-29 03:13:35 UTC (rev 18486)
@@ -490,6 +490,22 @@
for(i=1;i<=E->lmesh.nno;i++)
gzprintf(gzout,"%12.4e\n",E->C[i]);
gzclose(gzout);
+
+
+
+ if(E->tracers_add_flavors){
+ sprintf(output_file,"%s/%d/cf.%d.%d.gz",
+ E->control.data_file2,file_number, E->parallel.me,file_number);
+ gzout=safe_gzopen(output_file,"w");
+ gzprintf(gzout,"%3d %7d\n",1,E->lmesh.nno);
+ for(i=1;i <= E->lmesh.nno;i++) {
+ for(j=0;j < E->tracers_add_flavors;j++)
+ gzprintf(gzout,"%3i ",E->CF[j][i]);
+ gzprintf(gzout,"\n");
+ }
+ gzclose(gzout);
+ }
+
}
/*
end vtk out
Modified: mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c 2011-05-29 00:17:47 UTC (rev 18485)
+++ mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c 2011-05-29 03:13:35 UTC (rev 18486)
@@ -1015,3 +1015,36 @@
a[1][1] -= trace;
a[2][2] -= trace;
}
+/*
+
+compute the distance between the location x1,x2,x3 and the ndoe lnode
+where those coordinates are theta,phi,r for spherical and x,y,z for Cartesian
+
+ still need to test this properly
+
+*/
+double distance_to_node(double x1, double x2, double x3,struct All_variables *E, int lnode)
+{
+ double tmp1,tmp2,tmp3,dr;
+ if(E->control.Rsphere){
+ tmp1 = (E->SX[1][lnode] - x1)/2.0;
+ tmp1 = sin(tmp1);
+ tmp1 = tmp1 * tmp1;
+ tmp2 = (x2-E->SX[2][lnode])/2.0;
+ tmp2 = sin(tmp2);
+ tmp2 = tmp2 * tmp2;
+ tmp3 = tmp1;
+ tmp3 += sin(x1) * sin(E->SX[1][lnode]) * tmp2;
+ tmp3 = sqrt(tmp3);
+
+ tmp3 = 2.0*asin(tmp3); /* horizontal distance */
+ dr = E->SX[3][lnode]-x3; /* vertical distance */
+ return sqrt(tmp3*tmp3 + dr*dr);
+ }else{
+ tmp1 = x1 - E->X[1][lnode];tmp2 = tmp1*tmp1;
+ tmp1 = x2 - E->X[2][lnode];tmp2 += tmp1*tmp1;
+ tmp1 = x3 - E->X[3][lnode];tmp2 += tmp1*tmp1;
+ return sqrt(tmp2);
+ }
+
+}
Modified: mc/3D/CitcomCU/trunk/src/Parallel_related.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Parallel_related.c 2011-05-29 00:17:47 UTC (rev 18485)
+++ mc/3D/CitcomCU/trunk/src/Parallel_related.c 2011-05-29 03:13:35 UTC (rev 18486)
@@ -1238,6 +1238,82 @@
return;
}
+void exchange_node_int(struct All_variables *E, int *U, int lev)
+{
+ int target_proc, kk, i, j, k, idb;
+ static int *S[27], *R[27];
+
+ static int been_here = 0;
+ static int sizeofk;
+ const int levmax = E->mesh.levmax;
+
+ MPI_Status status[100];
+ MPI_Request request[100];
+
+ if(been_here == 0){
+ sizeofk = 0;
+ for(i = 1; i <= E->mesh.nsd; i++)
+ for(k = 1; k <= 2; k++){
+ if(E->parallel.NUM_PASS[levmax].bound[i][k])
+ sizeofk = max(sizeofk, (1 + E->parallel.NUM_NODE[levmax].pass[i][k]) * sizeof(int));
+ }
+ for(k = 1; k <= 2; k++){
+ S[k] = (int *)malloc(sizeofk);
+ R[k] = (int *)malloc(sizeofk);
+ }
+ been_here++;
+ }
+ for(i = 1; i <= E->mesh.nsd; i++)
+ {
+ idb = 0;
+ for(k = 1; k <= 2; k++)
+ if(E->parallel.NUM_PASS[levmax].bound[i][k])
+ {
+ for(j = 1; j <= E->parallel.NUM_NODE[lev].pass[i][k]; j++)
+ S[k][j - 1] = U[E->parallel.EXCHANGE_NODE[lev][j].pass[i][k]];
+
+ target_proc = E->parallel.PROCESSOR[lev].pass[i][k];
+ if(target_proc != E->parallel.me)
+ {
+ idb++;
+ MPI_Isend(S[k], E->parallel.NUM_NODE[lev].pass[i][k], MPI_INT, target_proc, 1, MPI_COMM_WORLD, &request[idb - 1]);
+ }
+ } /* for k */
+ for(k = 1; k <= 2; k++)
+ if(E->parallel.NUM_PASS[levmax].bound[i][k])
+ {
+
+ target_proc = E->parallel.PROCESSOR[lev].pass[i][k];
+ if(target_proc != E->parallel.me)
+ {
+ idb++;
+ MPI_Irecv(R[k], E->parallel.NUM_NODE[lev].pass[i][k], MPI_INT, target_proc, 1, MPI_COMM_WORLD, &request[idb - 1]);
+ }
+ } /* for k */
+ MPI_Waitall(idb, request, status);
+
+ for(k = 1; k <= 2; k++){
+ if(E->parallel.NUM_PASS[levmax].bound[i][k])
+ {
+ target_proc = E->parallel.PROCESSOR[lev].pass[i][k];
+ if(target_proc != E->parallel.me)
+ {
+ for(j = 1; j <= E->parallel.NUM_NODE[lev].pass[i][k]; j++)
+ U[E->parallel.EXCHANGE_NODE[lev][j].pass[i][k]] += R[k][j - 1];
+ }
+ else
+ {
+ kk = 1;
+ if(k == 1)
+ kk = 2;
+ for(j = 1; j <= E->parallel.NUM_NODE[lev].pass[i][k]; j++)
+ U[E->parallel.EXCHANGE_NODE[lev][j].pass[i][k]] += S[kk][j - 1];
+ }
+ } /* for k */
+ }
+ } /* for dim */
+ return;
+}
/* ========================== */
Modified: mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2011-05-29 00:17:47 UTC (rev 18485)
+++ mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2011-05-29 03:13:35 UTC (rev 18486)
@@ -156,6 +156,8 @@
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_int("pdepv_for_flavor",&(E->viscosity.pdepv_for_flavor),"0",m); /* byerlee only for certain flavor */
+
input_boolean("CDEPV",&(E->viscosity.CDEPV),"off",m);
input_boolean("cdepv_absolute",&(E->viscosity.cdepv_absolute),"off",m); /* make comp dep viscosity absolute */
@@ -337,7 +339,6 @@
if(E->viscosity.CDEPV)
visc_from_C(E, visc, evisc, propogate);
- /* two different kinds of implementations of powerlaw */
if(E->viscosity.SDEPV)
visc_from_S(E, visc, evisc, propogate);
@@ -1230,9 +1231,10 @@
static int visited = 0;
float scale,stress_magnitude,depth,exponent1,eta_old,eta_old2,eta_new;
float *eedot;
- float zzz,zz[9];
+ float zzz,zz[9],wmax,weight;
+ int point_flavor,tfn[9];
float tau,tau2,ettby,ettnew;
- int m,l,z,jj,kk,i;
+ int m,l,z,jj,kk,i,node,tepoints,epoints,npc,tnpc;
static float ndz_to_m;
#ifdef DEBUG
FILE *out;
@@ -1240,9 +1242,11 @@
const int vpts = vpoints[E->mesh.nsd];
const int nel = E->lmesh.nel;
const int ends = enodes[E->mesh.nsd];
+ epoints = nel * vpts;
eedot = (float *) safe_malloc((2+nel)*sizeof(float));
-
+
+
#ifdef DEBUG
out=fopen("tmp.visc","w");
#endif
@@ -1269,6 +1273,11 @@
E->viscosity.plasticity_trans,
E->viscosity.plasticity_viscosity_offset);
fprintf(stderr,"\tpsrw: %i\n",E->viscosity.psrw);
+
+ if(E->viscosity.pdepv_for_flavor){
+ if((E->tracers_add_flavors < 1) || (E->control.composition == 0))
+ myerror("Byerlee is to apply only to certain flavor, but no flavor is set",E);
+ }
}
/*
get strain rates for all elements
@@ -1282,133 +1291,175 @@
else
strain_rate_2_inv(E,eedot,1);
}
+
if(E->viscosity.psrw){
/* strain-rate weakening */
+ npc=0;
for(i=1;i <= nel;i++) {
l = E->mat[i] - 1; /* material of element */
for(kk=1;kk <= ends;kk++){/* loop through integration points*/
+ node = E->ien[i].node[kk];
if(E->control.Rsphere)
- zz[kk] = (1.0 - E->SX[3][E->ien[i].node[kk]]);
+ zz[kk] = (1.0 - E->SX[3][node]);
else
- zz[kk] = (1.0 - E->X[3][E->ien[i].node[kk]]);
+ zz[kk] = (1.0 - E->X[3][node]);
if(E->viscosity.plasticity_dimensional)
zz[kk] *= ndz_to_m; /* scale to meters */
+ if(E->viscosity.pdepv_for_flavor)
+ tfn[kk]= E->CF[0][node];
}
for(jj=1;jj <= vpts;jj++) {
- zzz = 0.0;
+ zzz = wmax = 0.0;
for(kk=1;kk <= ends;kk++) {
- zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+ weight = E->N.vpt[GNVINDEX(kk,jj)];
+ zzz += zz[kk] * weight;
+ if(E->viscosity.pdepv_for_flavor){
+ if(weight > wmax){
+ wmax = weight;
+ point_flavor = tfn[kk]; /* use the most important nodal flavor */
+ }
+ }
}
- if(E->viscosity.plasticity_dimensional){
- tau = (E->viscosity.abyerlee[l] * zzz +
- E->viscosity.bbyerlee[l]) * E->viscosity.lbyerlee[l];
- tau /= E->monitor.tau_scale;
- }else{
- tau = E->viscosity.abyerlee[l] * zzz + E->viscosity.bbyerlee[l];
- tau = min(tau, E->viscosity.lbyerlee[l]);
+ if((E->viscosity.pdepv_for_flavor==0)||(point_flavor == E->viscosity.pdepv_for_flavor)){
+ npc++;
+ /* plasticity applies */
+ if(E->viscosity.plasticity_dimensional){
+ tau = (E->viscosity.abyerlee[l] * zzz +
+ E->viscosity.bbyerlee[l]) * E->viscosity.lbyerlee[l];
+ tau /= E->monitor.tau_scale;
+ }else{
+ tau = E->viscosity.abyerlee[l] * zzz + E->viscosity.bbyerlee[l];
+ tau = min(tau, E->viscosity.lbyerlee[l]);
+ }
+ if((visited > 1) && (tau < 1e15)){
+ tau2 = tau * tau;
+ eta_old = EEta[ (i-1)*vpts + jj ] ;
+ eta_old2 = eta_old * eta_old;
+ eta_new = (tau2 * eta_old)/(tau2 + 2.0 * eta_old2 * eedot[i]);
+ EEta[ (i-1)*vpts + jj ] = ettnew;
+ //if(E->parallel.me==0)fprintf(stderr,"tau: %11.4e eII: %11.4e eta_old: %11.4e eta_new: %11.4e\n",tau, eedot[i],eta_old,eta_new);
+ }
+
}
- if((visited > 1) && (tau < 1e15)){
- tau2 = tau * tau;
- eta_old = EEta[ (i-1)*vpts + jj ] ;
- eta_old2 = eta_old * eta_old;
- eta_new = (tau2 * eta_old)/(tau2 + 2.0 * eta_old2 * eedot[i]);
- EEta[ (i-1)*vpts + jj ] = ettnew;
- //if(E->parallel.me==0)fprintf(stderr,"tau: %11.4e eII: %11.4e eta_old: %11.4e eta_new: %11.4e\n",tau, eedot[i],eta_old,eta_new);
- }
}
}
}else{
/* regular plasticity */
-
+ npc=0;
for(i=1;i <= nel;i++) {
/*
loop through all elements
*/
l = E->mat[i] - 1; /* material of element */
for(kk=1;kk <= ends;kk++){/* loop through integration points*/
+ node = E->ien[i].node[kk];
/* depth in meters */
if(E->control.Rsphere)
- zz[kk] = (1.0 - E->SX[3][E->ien[i].node[kk]]);
+ zz[kk] = (1.0 - E->SX[3][node]);
else
- zz[kk] = (1.0 - E->X[3][E->ien[i].node[kk]]);
+ zz[kk] = (1.0 - E->X[3][node]);
if(E->viscosity.plasticity_dimensional)
zz[kk] *= ndz_to_m; /* scale to meters */
+ if(E->viscosity.pdepv_for_flavor) /* nodal flavors */
+ tfn[kk]= E->CF[0][node];
}
for(jj=1;jj <= vpts;jj++) {
/* loop over nodes in element */
- zzz = 0.0;
- for(kk=1;kk <= ends;kk++) {
+ zzz = wmax = 0.0;
+ for(kk=1;kk <= ends;kk++) { /* add the nodal contribution */
/*
depth [m]
*/
- zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+ weight = E->N.vpt[GNVINDEX(kk,jj)];
+ zzz += zz[kk] * weight;
+ if(E->viscosity.pdepv_for_flavor){
+ if(weight > wmax){
+ wmax = weight;
+ point_flavor = tfn[kk]; /* use the most important nodal flavor */
+ }
+ }
}
- if(E->viscosity.plasticity_dimensional){
- /* byerlee type */
-
+ if((E->viscosity.pdepv_for_flavor==0)||(point_flavor == E->viscosity.pdepv_for_flavor)){
+ npc++;
/*
- yield stress in [Pa]
+ apply plasticity for this integration point if
+ pdepv_for_flavor is zero, or if the flavor matches
+
*/
- tau=(E->viscosity.abyerlee[l] * zzz + E->viscosity.bbyerlee[l]) * E->viscosity.lbyerlee[l];
+
+ if(E->viscosity.plasticity_dimensional){
+ /* byerlee type */
+
+ /*
+ yield stress in [Pa]
+ */
+ tau=(E->viscosity.abyerlee[l] * zzz + E->viscosity.bbyerlee[l]) * E->viscosity.lbyerlee[l];
+ /*
+ scaled stress
+ */
+ tau /= E->monitor.tau_scale;
+ }else{
+ tau = E->viscosity.abyerlee[l] * zzz + E->viscosity.bbyerlee[l];
+ tau = min(tau, E->viscosity.lbyerlee[l]);
+ }
/*
- scaled stress
+
+ `byerlee viscosity' : tau = 2 eps eta, this is non-dim
+ plus some offset as in Stein et al.
+
*/
- tau /= E->monitor.tau_scale;
- }else{
+ ettby = tau/(2.0 * (eedot[i]+1e-7)) + E->viscosity.plasticity_viscosity_offset;
+ /*
+
+ decide on the plasticity transition
- tau = E->viscosity.abyerlee[l] * zzz + E->viscosity.bbyerlee[l];
- tau = min(tau, E->viscosity.lbyerlee[l]);
- }
- /*
-
- `byerlee viscosity' : tau = 2 eps eta, this is non-dim
- plus some offset as in Stein et al.
+ */
+ if(E->viscosity.plasticity_trans){
+ /*
+ eta = 1./(1./eta(k)+1./ettby)
+ */
+
+ ettnew = 1.0/(1.0/EEta[ (i-1)*vpts + jj ] + 1.0/ettby);
+ //fprintf(stderr,"a: %g %g %g\n",EEta[ (i-1)*vpts + jj ],ettby,ettnew);
+ }else{
+ /*
+ min(\eta_p, \eta_visc )
+ */
+ ettnew = min(EEta[ (i-1)*vpts + jj ],ettby);
+ //fprintf(stderr,"m: %g %g %g\n",EEta[ (i-1)*vpts + jj ],ettby,ettnew);
+ }
- */
- ettby = tau/(2.0 * (eedot[i]+1e-7)) + E->viscosity.plasticity_viscosity_offset;
- /*
-
- decide on the plasticity transition
+#ifdef DEBUG
+ /* output format
+
+ z[km] tau[MPa] eps[s^-1] eta_b[Pas] eta_T[Pas] eta_c[Pas]
-
- */
- if(E->viscosity.plasticity_trans){
- /*
- eta = 1./(1./eta(k)+1./ettby)
*/
-
- ettnew = 1.0/(1.0/EEta[ (i-1)*vpts + jj ] + 1.0/ettby);
- //fprintf(stderr,"a: %g %g %g\n",EEta[ (i-1)*vpts + jj ],ettby,ettnew);
- }else{
- /*
- min(\eta_p, \eta_visc )
- */
- ettnew = min(EEta[ (i-1)*vpts + jj ],ettby);
- //fprintf(stderr,"m: %g %g %g\n",EEta[ (i-1)*vpts + jj ],ettby,ettnew);
- }
-#ifdef DEBUG
- /* output format
-
- z[km] tau[MPa] eps[s^-1] eta_b[Pas] eta_T[Pas] eta_c[Pas]
-
- */
- if(visited)
- fprintf(out,"%10.2f %17.4e %17.4e %17.4e %17.4e %17.4e\n",
- zzz/1e3,tau*E->monitor.tau_scale/1e6,
- eedot[i]/E->monitor.time_scale,
- ettby*E->data.ref_viscosity,
- EEta[ (i-1)*vpts + jj ]*E->data.ref_viscosity,
- ettnew*E->data.ref_viscosity);
+ if(visited)
+ fprintf(out,"%10.2f %17.4e %17.4e %17.4e %17.4e %17.4e\n",
+ zzz/1e3,tau*E->monitor.tau_scale/1e6,
+ eedot[i]/E->monitor.time_scale,
+ ettby*E->data.ref_viscosity,
+ EEta[ (i-1)*vpts + jj ]*E->data.ref_viscosity,
+ ettnew*E->data.ref_viscosity);
#endif
- // if(visited)
- // fprintf(stderr,"%11.4e %11.4e %11.4e %11.4e\n",ettnew,EEta[ (i-1)*vpts + jj ] ,ettby,eedot[i]);
- EEta[ (i-1)*vpts + jj ] = ettnew;
- }
+ // if(visited)
+ // fprintf(stderr,"%11.4e %11.4e %11.4e %11.4e\n",ettnew,EEta[ (i-1)*vpts + jj ] ,ettby,eedot[i]);
+ EEta[ (i-1)*vpts + jj ] = ettnew;
+ } /* plasticity applies, for this flavor branch */
+ } /* end integration point loop */
} /* end regular plasticity */
}
+ if(E->viscosity.pdepv_for_flavor){
+ MPI_Allreduce(&epoints,&tepoints,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
+ MPI_Allreduce(&npc,&tnpc,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
+ if(E->parallel.me==0)fprintf(stderr,"applied plasticity for flavor %i to %10i out of %10i integration points (%5.2f%%)\n",
+ E->viscosity.pdepv_for_flavor,tnpc,tepoints,(float)tnpc/(float)tepoints*100);
+
+ }
#ifdef DEBUG
fclose(out);
#endif
Modified: mc/3D/CitcomCU/trunk/src/global_defs.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/global_defs.h 2011-05-29 00:17:47 UTC (rev 18485)
+++ mc/3D/CitcomCU/trunk/src/global_defs.h 2011-05-29 03:13:35 UTC (rev 18486)
@@ -762,7 +762,7 @@
float ggrd_slab_theta_bound[4];
struct ggrd_gt ggrd_ss_grd[4];
int ggrd_mat_limit_prefactor,ggrd_mat_is_3d;
- char ggrd_mat_depth_file[1000];
+ char ggrd_mat_depth_file[1000],ggrd_flavor_gfile[1000],ggrd_flavor_dfile[1000];
#endif
int adi_heating, visc_heating;
@@ -946,8 +946,9 @@
int *traces_leave_index;
int *CElement;
- int *tflavors; /* this should also be unsigned short */
+ int **tflavors; /* this should also be unsigned short */
int tracers_add_flavors;
+ int *tmaxflavor;
int *RG[4];
double *XRG[4];
@@ -978,6 +979,7 @@
float *Vi, *EVi;
float *diffusivity, *expansivity;
float *T, *C, *CE, *buoyancy;
+ int **CF; /* tracer flavors */
float *Tdot;
float *VI[MAX_LEVELS]; /* viscosity has to soak down to all levels */
float *EVI[MAX_LEVELS]; /* element viscosity has to soak down to all levels */
Modified: mc/3D/CitcomCU/trunk/src/prototypes.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/prototypes.h 2011-05-29 00:17:47 UTC (rev 18485)
+++ mc/3D/CitcomCU/trunk/src/prototypes.h 2011-05-29 03:13:35 UTC (rev 18486)
@@ -170,6 +170,8 @@
void ggrd_read_mat_from_file(struct All_variables *);
void ggrd_solve_eigen3x3(double [3][3], double [3], double [3][3], struct All_variables *);
void ggrd_read_anivisc_from_file(struct All_variables *);
+void ggrd_deal_with_composition_input(struct All_variables *, int );
+
/* Global_operations.c */
void return_horiz_sum(struct All_variables *, float *, float *, int);
void return_horiz_ave(struct All_variables *, float *, float *);
@@ -373,3 +375,8 @@
int layers(struct All_variables *, float);
int weak_zones(struct All_variables *, int, float);
float boundary_thickness(struct All_variables *, float *);
+double distance_to_node(double , double , double ,struct All_variables *, int );
+void assign_flavor_to_tracer_based_on_node(struct All_variables *, int , int );
+void exchange_node_int(struct All_variables *, int *, int );
+void get_CF_from_markers(struct All_variables *, int **);
+void assign_flavor_to_tracer_from_grd(struct All_variables *);
Modified: mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h 2011-05-29 00:17:47 UTC (rev 18485)
+++ mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h 2011-05-29 03:13:35 UTC (rev 18486)
@@ -144,11 +144,11 @@
int cdepv_absolute; /* make the composition an absolute viscosity,
not a prefactor */
float pre_comp[CITCOM_CU_VISC_MAXLAYER]; /* prefactors */
+ int pdepv_for_flavor;
-
int weak_blobs;
float weak_blobx[CITCOM_CU_VISC_MAXLAYER];
float weak_bloby[CITCOM_CU_VISC_MAXLAYER];
More information about the CIG-COMMITS
mailing list