[cig-commits] r19332 - mc/3D/CitcomCU/trunk/src
becker at geodynamics.org
becker at geodynamics.org
Sat Jan 7 12:05:13 PST 2012
Author: becker
Date: 2012-01-07 12:05:13 -0800 (Sat, 07 Jan 2012)
New Revision: 19332
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/global_defs.h
Log:
Implemented a simplified version of the absolute tracer method,
just for testing purposes.
Modified: mc/3D/CitcomCU/trunk/src/Composition_adv.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Composition_adv.c 2012-01-07 18:11:40 UTC (rev 19331)
+++ mc/3D/CitcomCU/trunk/src/Composition_adv.c 2012-01-07 20:05:13 UTC (rev 19332)
@@ -549,28 +549,53 @@
{
element[E->C12[imark]][E->CElement[imark]]++;
}
+ if(E->tracers_assign_dense_only){
+ /*
+ generally, a bad idea, counting dense only
+ */
+ for(el = 1; el <= nel; el++)
+ {
+ temp1 = element[1][el];
+ if(element[1][el]){
+ temp3 = temp1 / (float)E->advection.markers_per_ele; /* elemental C */
+ if(temp3 > 1)
+ temp3 = 1.0; /* against tracer clumping */
+ }else
+ temp3 = 0.0;
+ for(j = 1; j <= ends; j++)
+ {
+ node = E->ien[el].node[j];
+ C[node] += E->TWW[lev][el].node[j] * temp3;
+ }
+ E->CE[el] = temp3;
+ }
+ }else{
- 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); /* elemental C */
- else
- temp3 = E->CE[el]; /* elemental C */
- for(j = 1; j <= ends; j++)
+ /*
+ ratio method
+ */
+ 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); /* elemental C */
+ else
+ temp3 = E->CE[el]; /* elemental C */
+ for(j = 1; j <= ends; j++)
{
- node = E->ien[el].node[j];
- C[node] += E->TWW[lev][el].node[j] * temp3;
+ node = E->ien[el].node[j];
+ C[node] += E->TWW[lev][el].node[j] * temp3;
}
- E->CE[el] = temp3;
+ E->CE[el] = temp3;
+ }
}
exchange_node_f20(E, C, E->mesh.levmax);
for(node = 1; node <= nno; node++)
{
- C[node] = C[node] * E->Mass[node];
+ C[node] = C[node] * E->Mass[node];
}
return;
Modified: mc/3D/CitcomCU/trunk/src/Convection.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Convection.c 2012-01-07 18:11:40 UTC (rev 19331)
+++ mc/3D/CitcomCU/trunk/src/Convection.c 2012-01-07 20:05:13 UTC (rev 19332)
@@ -53,9 +53,16 @@
input_boolean("composition", &(E->control.composition), "0", E->parallel.me);
input_int("tracers_add_flavors", &(E->tracers_add_flavors), "0", E->parallel.me);
+
+ input_float("tracers_assign_dense_fraction",&(E->tracers_assign_dense_fraction),"1.0",E->parallel.me);
+ if(E->tracers_assign_dense_fraction < 1.0){
+ force_report(E,"WARNING: only assigning tracers to elements which are predominantly dense");
+ if(E->parallel.me == 0)fprintf(stderr,"assuming a fraction of %g%%\n",E->tracers_assign_dense_fraction * 100.0);
+ E->tracers_assign_dense_only = 1;
+ }else{
+ E->tracers_assign_dense_only = 0;
+ }
- input_boolean("tracers_assign_dense_only", &(E->tracers_assign_dense_only),"0",E->parallel.me);
-
if(E->control.composition)
E->next_buoyancy_field = PG_timestep_particle;
else
@@ -178,8 +185,11 @@
{
if(E->tracers_assign_dense_only){
- force_report(E,"WARNING: assigning only dense tracers, generally not a good idea!");
- E->advection.markers = E->advection.markers_per_ele * E->mesh.nel;
+ if(E->parallel.me == 0)
+ fprintf(stderr,"WARNING: assigning only to %g%% of elements, not a good idea!\n",
+ E->tracers_assign_dense_fraction*100);
+ E->advection.markers = (int)((float)E->advection.markers_per_ele * (float) E->mesh.nel *
+ E->tracers_assign_dense_fraction);
}else{
E->advection.markers = E->advection.markers_per_ele * E->mesh.nel;
@@ -647,112 +657,173 @@
return;
}
-void convection_initial_markers(struct All_variables *E,int use_element_nodes_for_init_c)
+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, 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];
- //FILE *fp;
- float temp;
+ //int el, i, j, k, p, node, ii, jj;
+ int el, ntracer,j,i,node,k;
+ //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];
+ //FILE *fp;
+ float temp,frac,locx[3],locp[3];
- const int dims = E->mesh.nsd;
- const int ends = enodes[dims];
+ const int dims = E->mesh.nsd;
+ const int ends = enodes[dims];
- if(E->control.CART3D)
- {
- ntracer = 0;
- do
- {
- x = drand48() * (E->XG2[1] - E->XG1[1]);
- y = drand48() * (E->XG2[2] - E->XG1[2]);
- z = drand48() * (E->XG2[3] - E->XG1[3]);
+ if(E->tracers_assign_dense_only){
+ /*
+ assign only tracers within elements that are on average dense
+
+ this is, in general, not a good idea, but i added it for
+ testing purposes
- 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]))
- {
- ntracer++;
- E->XMC[1][ntracer] = x;
- E->XMC[2][ntracer] = y;
- E->XMC[3][ntracer] = z;
+ */
+ ntracer=0;
+ for(el = 1; el <= E->lmesh.nel; el++){
+ /* get mean composition */
+ for(temp=0.0,j = 1; j <= ends; j++)
+ temp += E->C[E->ien[el].node[j]];
+ temp /= ends;
+ if(temp > 0.5){
+ for(j=0;j < E->advection.markers_per_ele;j++){
+ ntracer++;
- 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[ntracer] = 1;
- else
- E->C12[ntracer] = 0;
- }else{ /* use depth */
- if(E->XMC[3][ntracer] > E->viscosity.zcomp)
- E->C12[ntracer] = 0;
- else
- E->C12[ntracer] = 1;
- }
- }
- } while(ntracer < E->advection.markers);
+ if(ntracer > E->advection.markers)
+ myerror("only dense assign: out of markers",E);
+ /* get a random location within the element */
+ if(E->control.CART3D){ /* cartesian */
+ x = y = z = 0;
+ for(k=1;k <= ends;k++){
+ node = E->ien[el].node[k];
+ frac = drand48();
+ x += frac * E->X[1][node];
+ y += frac * E->X[2][node];
+ z += frac * E->X[3][node];
+ }
+ E->XMC[1][ntracer] = x / (float)ends;
+ E->XMC[2][ntracer] = y / (float)ends;
+ E->XMC[3][ntracer] = z / (float)ends;
+ }else{
+ x = y = z = 0;
+ for(k=1;k <= ends;k++){
+ node = E->ien[el].node[k];
+ frac = drand48();
+ rtp2xyz((float)E->SX[3][node],(float)E->SX[1][node],(float)E->SX[2][node],locx); /* convert to cart */
+ x += frac * locx[0];
+ y += frac * locx[1];
+ z += frac * locx[2];
+ }
+ locx[0] = x / (float)ends; locx[1] = y /(float)ends; locx[2] = z/(float)ends;
+ xyz2rtp(locx[0],locx[1],locx[2],locp); /* back to spherical */
+ E->XMC[1][ntracer] = locp[1]; /* theta */
+ E->XMC[2][ntracer] = locp[2]; /* phi */
+ E->XMC[3][ntracer] = locp[0]; /* r */
+ }
+ E->C12[ntracer] = 1;
+ E->CElement[ntracer] = el;
}
- else if(E->control.Rsphere)
- {
- ntracer = 0;
- do
- {
- x = (drand48() - 0.5) * 2.0;
- y = drand48();
-// y = (drand48()-0.5)*2.0;
- z = (drand48() - 0.5) * 2.0;
- r = sqrt(x * x + y * y + z * z);
- t = acos(z / r);
- 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]))
- {
- 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[ntracer] = 1;
- else
- E->C12[ntracer] = 0;
- }else{ /* use depth */
- if(r > E->viscosity.zcomp)
- E->C12[ntracer] = 0;
- else
- E->C12[ntracer] = 1;
- }
+ }
+ }
-
- }
- } while(ntracer < E->advection.markers);
- }
+ }else{
+ /* regular operation */
+
+
+ if(E->control.CART3D)
+ {
+ ntracer = 0;
+ do
+ {
+ x = drand48() * (E->XG2[1] - E->XG1[1]);
+ y = drand48() * (E->XG2[2] - E->XG1[2]);
+ z = drand48() * (E->XG2[3] - E->XG1[3]);
+
+ 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]))
+ {
+ ntracer++;
+ E->XMC[1][ntracer] = x;
+ E->XMC[2][ntracer] = y;
+ E->XMC[3][ntracer] = z;
+
+ 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[ntracer] = 1;
+ else
+ E->C12[ntracer] = 0;
+ }else{ /* use depth */
+ if(E->XMC[3][ntracer] > E->viscosity.zcomp)
+ E->C12[ntracer] = 0;
+ else
+ E->C12[ntracer] = 1;
+ }
+ }
+ } while(ntracer < E->advection.markers);
+ }
+ else if(E->control.Rsphere)
+ {
+ ntracer = 0;
+ do
+ {
+ x = (drand48() - 0.5) * 2.0;
+ y = drand48();
+ // y = (drand48()-0.5)*2.0;
+ z = (drand48() - 0.5) * 2.0;
+ r = sqrt(x * x + y * y + z * z);
+ t = acos(z / r);
+ 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]))
+ {
+ 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[ntracer] = 1;
+ else
+ E->C12[ntracer] = 0;
+ }else{ /* use depth */
+ if(r > E->viscosity.zcomp)
+ E->C12[ntracer] = 0;
+ else
+ E->C12[ntracer] = 1;
+ }
+
+
+ }
+ } while(ntracer < E->advection.markers);
+ }
+ } /* end regular operation */
- /* assign tracers based on ggrd */
- if(E->tracers_add_flavors){
+ /* assign tracers based on ggrd */
+ if(E->tracers_add_flavors){
#ifdef USE_GGRD
- assign_flavor_to_tracer_from_grd(E);
+ assign_flavor_to_tracer_from_grd(E);
#else
- myerror("convection_initial_markers: flavor init requires GGRD compilation",E);
+ 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;
+ }
+ /*
+ 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
*/
Modified: mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2012-01-07 18:11:40 UTC (rev 19331)
+++ mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2012-01-07 20:05:13 UTC (rev 19332)
@@ -64,6 +64,8 @@
char *char_dummy="";
+ int c1_local,c1_total,n_total,int_unity = 1;
+
/* twb additions */
double rho_prem;
char pfile[1000];
@@ -405,10 +407,30 @@
if(E->C[i]<0)E->C[i]=0;
}
+ /* compute how many dense nodes were assigned */
+ c1_local = c1_total = n_total = 0;
+ for(i=1;i <= E->lmesh.nno;i++){
+ if(E->C[i] > 0.5)
+ c1_local++;
+ }
+ MPI_Allreduce(&(E->lmesh.nno), &n_total,1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+ MPI_Allreduce(&c1_local, &c1_total,1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+ E->tracers_dense_frac = (float)c1_total/(float)n_total;
+ /* check if we restrict assignment */
+ if(E->parallel.me==0)
+ fprintf(stderr,"assigned C>0.5 %i/%i times out of %i/%i, %.1f%%\n",
+ c1_local,c1_total,E->lmesh.nno,n_total,E->tracers_dense_frac*100.);
+ if(E->tracers_assign_dense_only){
+ if(E->parallel.me == 0)
+ fprintf(stderr,"compares with restricted set dense fraction estimate of %g%%\n",
+ E->tracers_assign_dense_fraction * 100);
+ if(E->tracers_assign_dense_fraction < E->tracers_dense_frac)
+ myerror("increase the dense fraction for assignment, too small",E);
+ }
+
if(E->control.composition)
convection_initial_markers(E,1);
} // end for restart==0
-
else if(E->control.restart)
{
#ifdef USE_GZDIR
Modified: mc/3D/CitcomCU/trunk/src/global_defs.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/global_defs.h 2012-01-07 18:11:40 UTC (rev 19331)
+++ mc/3D/CitcomCU/trunk/src/global_defs.h 2012-01-07 20:05:13 UTC (rev 19332)
@@ -961,6 +961,8 @@
int *tmaxflavor;
int tracers_assign_dense_only;
+ float tracers_assign_dense_fraction,tracers_dense_frac;
+
int *RG[4];
double *XRG[4];
double *XP[4], XG1[4], XG2[4];
More information about the CIG-COMMITS
mailing list