[cig-commits] r5923 - mc/3D/CitcomS/trunk/lib
gurnis at geodynamics.org
gurnis at geodynamics.org
Mon Jan 29 12:46:29 PST 2007
Author: gurnis
Date: 2007-01-29 12:46:28 -0800 (Mon, 29 Jan 2007)
New Revision: 5923
Modified:
mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
mc/3D/CitcomS/trunk/lib/BC_util.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Problem_related.c
Log:
Changed the seequence placement in the time step when the age grids are assimilated into the temperature field. Previously, going back to the original code contributed by CC, the assimilation occurred avery time the temperature BCs were set (meaning multiple times in a PC loop). Now assimilation happens only once at the end of the time step. This is placed at the end of the routine PG_time_step_solve. Once we have worked with this some more, we will want to expose this code at a higher level so that assimulation will work properly from either plain C code or the Pyre version. Currently only tested with Full.
Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-01-29 20:37:00 UTC (rev 5922)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-01-29 20:46:28 UTC (rev 5923)
@@ -168,6 +168,13 @@
free((void *) T1[m] );
free((void *) Tdot1[m] );
}
+
+ if(E->control.lith_age) {
+ lith_age_conform_tbc(E);
+ assimilate_lith_conform_bcs(E);
+ }
+
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/BC_util.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/BC_util.c 2007-01-29 20:37:00 UTC (rev 5922)
+++ mc/3D/CitcomS/trunk/lib/BC_util.c 2007-01-29 20:46:28 UTC (rev 5923)
@@ -52,8 +52,11 @@
void assimilate_lith_conform_bcs2(struct All_variables *);
if(E->control.lith_age) {
+ /*
+ This sequence now moved to end of PG_time_step_solve
lith_age_conform_tbc(E);
assimilate_lith_conform_bcs(E);
+ */
}
else
temperatures_conform_bcs2(E);
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-01-29 20:37:00 UTC (rev 5922)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-01-29 20:46:28 UTC (rev 5923)
@@ -141,8 +141,11 @@
set_sphere_harmonics (E);
- if(E->control.mat_control)
+ if(E->control.mat_control) {
+ if(E->parallel.me ==0) fprintf(stderr,"IN Instructions.c\n");
+ fflush(stderr);
read_mat_from_file(E);
+ }
else
construct_mat_group(E);
Modified: mc/3D/CitcomS/trunk/lib/Problem_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Problem_related.c 2007-01-29 20:37:00 UTC (rev 5922)
+++ mc/3D/CitcomS/trunk/lib/Problem_related.c 2007-01-29 20:46:28 UTC (rev 5923)
@@ -48,104 +48,7 @@
void read_mat_from_file(E)
struct All_variables *E;
{
- float find_age_in_MY();
-
- int nn,m,i,j,k,kk,el,lev,els,cap;
- int elx,ely,elz,e,elg,emax,gmax;
- float *VIP1,*VIP2;
-
- float age1,newage1,newage2;
- int nodea,nage;
-
- int llayer;
- int layers();
- FILE *fp,*fp1,*fp2,*fp3,*fp4;
- char output_file[255];
-
- const int dims=E->mesh.nsd,dofs=E->mesh.dof;
- const int ends=enodes[dims];
-
- elx=E->lmesh.elx;
- elz=E->lmesh.elz;
- ely=E->lmesh.ely;
-
- emax=E->mesh.elx*E->mesh.elz*E->mesh.ely;
- gmax=E->mesh.elx*E->mesh.ely;
-
- VIP1 = (float*) malloc ((gmax+1)*sizeof(float));
- VIP2 = (float*) malloc ((gmax+1)*sizeof(float));
-
-
-
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for (el=1; el<=elx*ely*elz; el++) {
- nodea = E->ien[m][el].node[2];
- llayer = layers(E,m,nodea);
- if (llayer) { /* for layers:1-lithosphere,2-upper, 3-trans, and 4-lower mantle */
- E->mat[m][el] = llayer;
- }
- }
-
-
- if(E->control.mat_control==1) {
-
- age1 = find_age_in_MY(E);
-
- nage=age1/1.;
- newage1=1.*nage;
-
- for(m=1;m<=E->sphere.caps_per_proc;m++) {
- cap = E->sphere.capid[m] - 1; /* capid: 1-12 */
-
- sprintf(output_file,"%s%0.0f.%d",E->control.mat_file,newage1,cap);
- fprintf(E->fp,"%s %f %s\n","open mat file newage1",newage1,output_file);
- fp1=fopen(output_file,"r");
- if (fp1 == NULL) {
- fprintf(E->fp,"(Problem_related #1) Cannot open %s\n",output_file);
- exit(8);
- }
-
- newage2=newage1+1.;
- sprintf(output_file,"%s%0.0f.%d",E->control.mat_file,newage2,cap);
- fprintf(E->fp,"%s %f %s\n","open mat file newage2",newage2,output_file);
- fp2=fopen(output_file,"r");
- if (fp2 == NULL) {
- fprintf(E->fp,"(Problem_related #2) Cannot open %s\n",output_file);
- exit(8);
- }
- } /*end m */
-
- for(i=1;i<=gmax;i++) {
- fscanf(fp1,"%d %f", &nn,&(VIP1[i]));
- fscanf(fp2,"%d %f", &nn,&(VIP2[i]));
- }
-
- fclose(fp1);
- fclose(fp2);
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for (k=1;k<=ely;k++)
- for (i=1;i<=elx;i++) {
- elg = E->lmesh.exs+i + (E->lmesh.eys+k-1)*E->mesh.elx;
-
- for (j=1;j<=elz;j++) {
- el = j + (i-1)*E->lmesh.elz + (k-1)*E->lmesh.elz*E->lmesh.elx;
-
- if(E->sx[m][3][E->ien[m][el].node[2]]>=E->sphere.ro-E->viscosity.zlith)
- E->VIP[m][el] = VIP1[elg]+(VIP2[elg]-VIP1[elg])/(newage2-newage1)*(age1-newage1);
-
- } /* end for j */
-
- } /* end for m */
-
- } /* end for E->control.mat==1 */
-
-
- /* mat output moved to Output.c */
-
- free ((void *) VIP1);
- free ((void *) VIP2);
-
+ (E->solver.read_input_files_for_timesteps)(E,3,1); /* read element material(3) and output(1) */
return;
}
More information about the cig-commits
mailing list