[cig-commits] r5929 - in mc/3D/CitcomS/branches/compressible:
CitcomS/Solver lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Tue Jan 30 11:10:01 PST 2007
Author: tan2
Date: 2007-01-30 11:10:00 -0800 (Tue, 30 Jan 2007)
New Revision: 5929
Modified:
mc/3D/CitcomS/branches/compressible/CitcomS/Solver/Solver.py
mc/3D/CitcomS/branches/compressible/lib/Advection_diffusion.c
mc/3D/CitcomS/branches/compressible/lib/BC_util.c
mc/3D/CitcomS/branches/compressible/lib/Full_read_input_from_files.c
mc/3D/CitcomS/branches/compressible/lib/Instructions.c
mc/3D/CitcomS/branches/compressible/lib/Lith_age.c
mc/3D/CitcomS/branches/compressible/lib/Problem_related.c
mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c
Log:
Porting r5921-5928 from trunk
Modified: mc/3D/CitcomS/branches/compressible/CitcomS/Solver/Solver.py
===================================================================
--- mc/3D/CitcomS/branches/compressible/CitcomS/Solver/Solver.py 2007-01-30 19:04:35 UTC (rev 5928)
+++ mc/3D/CitcomS/branches/compressible/CitcomS/Solver/Solver.py 2007-01-30 19:10:00 UTC (rev 5929)
@@ -104,6 +104,12 @@
else:
stream = None
+ # XXX: This is a heck
+ # Write controller inventory to pid file
+ print >> stream, "[CitcomS.controller]"
+ print >> stream, "monitoringFrequency =", application.controller.inventory.monitoringFrequency
+ print >> stream
+
self.setProperties(stream)
self.restart = self.inventory.ic.inventory.restart
Modified: mc/3D/CitcomS/branches/compressible/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Advection_diffusion.c 2007-01-30 19:04:35 UTC (rev 5928)
+++ mc/3D/CitcomS/branches/compressible/lib/Advection_diffusion.c 2007-01-30 19:10:00 UTC (rev 5929)
@@ -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/branches/compressible/lib/BC_util.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/BC_util.c 2007-01-30 19:04:35 UTC (rev 5928)
+++ mc/3D/CitcomS/branches/compressible/lib/BC_util.c 2007-01-30 19:10:00 UTC (rev 5929)
@@ -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/branches/compressible/lib/Full_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Full_read_input_from_files.c 2007-01-30 19:04:35 UTC (rev 5928)
+++ mc/3D/CitcomS/branches/compressible/lib/Full_read_input_from_files.c 2007-01-30 19:10:00 UTC (rev 5929)
@@ -48,9 +48,18 @@
int nox,noz,noy,nnn,nox1,noz1,noy1,lev;
int i,ii,ll,m,mm,j,k,n,nodeg,nodel,node,cap;
int intage, pos_age;
+ int nodea;
+ int nn,el;
const int dims=E->mesh.nsd;
+ int elx,ely,elz,elg,emax;
+ float *VIP1,*VIP2;
+ int *LL1, *LL2;
+
+ int llayer;
+ int layers();
+
nox=E->mesh.nox;
noy=E->mesh.noy;
noz=E->mesh.noz;
@@ -59,6 +68,12 @@
noy1=E->lmesh.noy;
lev=E->mesh.levmax;
+ elx=E->lmesh.elx;
+ elz=E->lmesh.elz;
+ ely=E->lmesh.ely;
+
+ emax=E->mesh.elx*E->mesh.elz*E->mesh.ely;
+
age=find_age_in_MY(E);
if (age < 0.0) { /* age is negative -> use age=0 for input files */
@@ -103,7 +118,7 @@
}
break;
- case 2: /* read ages for lithosphere tempperature boundary conditions */
+ case 2: /* read ages for lithosphere tempperature assimilation */
sprintf(output_file1,"%s%0.0f.%d",E->control.lith_age_file,newage1,cap);
sprintf(output_file2,"%s%0.0f.%d",E->control.lith_age_file,newage2,cap);
fp1=fopen(output_file1,"r");
@@ -128,6 +143,33 @@
}
break;
+ case 3: /* read element materials */
+
+ sprintf(output_file1,"%s%0.0f.%d",E->control.mat_file,newage1,cap);
+ sprintf(output_file2,"%s%0.0f.%d",E->control.mat_file,newage2,cap);
+ fp1=fopen(output_file1,"r");
+ if (fp1 == NULL) {
+ fprintf(E->fp,"(Problem_related #8) Cannot open %s\n",output_file1);
+ exit(8);
+ }
+ if (pos_age) {
+ fp2=fopen(output_file2,"r");
+ if (fp2 == NULL) {
+ fprintf(E->fp,"(Problem_related #9) Cannot open %s\n",output_file2);
+ exit(8);
+ }
+ }
+ if((E->parallel.me==0) && (output==1)) {
+ fprintf(E->fp,"Mat: Starting Age = %g, Elapsed time = %g, Current Age = %g\n",E->control.start_age,E->monitor.elapsed_time,age);
+ fprintf(E->fp,"Mat: File1 = %s\n",output_file1);
+ if (pos_age)
+ fprintf(E->fp,"Mat: File2 = %s\n",output_file2);
+ else
+ fprintf(E->fp,"Mat: File2 = No file inputted (negative age)\n");
+ }
+
+ break;
+
} /* end switch */
@@ -178,7 +220,7 @@
}
break;
- case 2: /* ages for lithosphere temperature boundary conditions */
+ case 2: /* ages for lithosphere temperature assimilation */
for(i=1;i<=noy;i++)
for(j=1;j<=nox;j++) {
node=j+(i-1)*nox;
@@ -195,6 +237,53 @@
if (pos_age) fclose(fp2);
break;
+ case 3: /* read element materials */
+
+ VIP1 = (float*) malloc ((emax+1)*sizeof(float));
+ VIP2 = (float*) malloc ((emax+1)*sizeof(float));
+ LL1 = (int*) malloc ((emax+1)*sizeof(int));
+ LL2 = (int*) malloc ((emax+1)*sizeof(int));
+
+
+ /* probably can be safely removed */
+ 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;
+ }
+ }
+ for(i=1;i<=emax;i++) {
+ fscanf(fp1,"%d %d %f", &nn,&(LL1[i]),&(VIP1[i]));
+ fscanf(fp2,"%d %d %f", &nn,&(LL2[i]),&(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++) {
+ for (j=1;j<=elz;j++) {
+ el = j + (i-1)*E->lmesh.elz + (k-1)*E->lmesh.elz*E->lmesh.elx;
+ elg = E->lmesh.ezs+j + (E->lmesh.exs+i-1)*E->mesh.elz + (E->lmesh.eys+k-1)*E->mesh.elz*E->mesh.elx;
+
+ E->VIP[m][el] = VIP1[elg]+(VIP2[elg]-VIP1[elg])/(newage2-newage1)*(age-newage1);
+ E->mat[m][el] = LL1[elg];
+
+ } /* end for j */
+ } /* end for i */
+ } /* end for k */
+ } /* end for m */
+
+ free ((void *) VIP1);
+ free ((void *) VIP2);
+ free ((void *) LL1);
+ free ((void *) LL2);
+
+ break;
+
} /* end switch */
} /* end for m */
Modified: mc/3D/CitcomS/branches/compressible/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Instructions.c 2007-01-30 19:04:35 UTC (rev 5928)
+++ mc/3D/CitcomS/branches/compressible/lib/Instructions.c 2007-01-30 19:10:00 UTC (rev 5929)
@@ -144,8 +144,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/branches/compressible/lib/Lith_age.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Lith_age.c 2007-01-30 19:04:35 UTC (rev 5928)
+++ mc/3D/CitcomS/branches/compressible/lib/Lith_age.c 2007-01-30 19:10:00 UTC (rev 5929)
@@ -341,7 +341,7 @@
/* NOW SET THE TEMPERATURES IN THE LITHOSPHERE IF CHANGING EVERY TIME STEP */
- if(E->monitor.solution_cycles>1 && E->control.lith_age_time) {
+ if(E->monitor.solution_cycles>0 && E->control.lith_age_time) {
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=noy;i++)
for(j=1;j<=nox;j++)
Modified: mc/3D/CitcomS/branches/compressible/lib/Problem_related.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Problem_related.c 2007-01-30 19:04:35 UTC (rev 5928)
+++ mc/3D/CitcomS/branches/compressible/lib/Problem_related.c 2007-01-30 19:10:00 UTC (rev 5929)
@@ -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;
}
Modified: mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c 2007-01-30 19:04:35 UTC (rev 5928)
+++ mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c 2007-01-30 19:10:00 UTC (rev 5929)
@@ -408,18 +408,22 @@
- E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
if(E->control.mat_control==1) {
- visc1 = E->VIP[m][i];
+ /* visc1 = E->VIP[m][i];
visc2 = 2.0/(1./visc1 + 1.);
tempa_exp = tempa*
- exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
- - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
+ exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
+ - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
visc1 = tempa*E->viscosity.max_value;
if(tempa_exp > visc1) tempa_exp=visc1;
EEta[m][ (i-1)*vpts + jj ] = visc2*tempa_exp;
- /* if(E->parallel.me == 0 && visc1 < 1.0e-03)
- fprintf(stderr,"%f %f %e %e %e\n",zzz,temp,visc1,visc2,
- EEta[m][ (i-1)*vpts + jj ]); */
- }
+ */
+ visc1 = E->VIP[m][i];
+ visc2 = tempa*
+ exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
+ - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
+ if(visc1 <= 0.95) visc2=visc1;
+ EEta[m][ (i-1)*vpts + jj ] = visc2;
+ }
}
}
More information about the cig-commits
mailing list