[cig-commits] r6229 - mc/3D/CitcomS/trunk/lib
gurnis at geodynamics.org
gurnis at geodynamics.org
Mon Mar 12 13:08:14 PDT 2007
Author: gurnis
Date: 2007-03-12 13:08:12 -0700 (Mon, 12 Mar 2007)
New Revision: 6229
Modified:
mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
mc/3D/CitcomS/trunk/lib/Lith_age.c
mc/3D/CitcomS/trunk/lib/Regional_lith_age_read_files.c
mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
Log:
Modified Regional_read_input_from_files.c so that it now opens and reads in material files. Material files are organized by elements. Previsouly only worked for Full
Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-03-12 16:35:30 UTC (rev 6228)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-03-12 20:08:12 UTC (rev 6229)
@@ -170,6 +170,7 @@
}
if(E->control.lith_age) {
+ if(E->parallel.me==0) fprintf(stderr,"PG_timestep_solve\n");
lith_age_conform_tbc(E);
assimilate_lith_conform_bcs(E);
}
Modified: mc/3D/CitcomS/trunk/lib/Lith_age.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Lith_age.c 2007-03-12 16:35:30 UTC (rev 6228)
+++ mc/3D/CitcomS/trunk/lib/Lith_age.c 2007-03-12 20:08:12 UTC (rev 6229)
@@ -75,6 +75,7 @@
gnox=E->mesh.nox;
gnoy=E->mesh.noy;
+ if (E->parallel.me == 0 ) fprintf(stderr,"INSIDE lith_age_init\n");
E->age_t=(float*) malloc((gnox*gnoy+1)*sizeof(float));
if(E->control.lith_age_time==1) {
@@ -287,6 +288,7 @@
output = 1;
E->control.lith_age_old_cycles = E->monitor.solution_cycles;
}
+ if (E->parallel.me == 0) fprintf(stderr,"INSIDE lith_age_conform_tbc\n");
(E->solver.lith_age_read_files)(E,output);
}
Modified: mc/3D/CitcomS/trunk/lib/Regional_lith_age_read_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_lith_age_read_files.c 2007-03-12 16:35:30 UTC (rev 6228)
+++ mc/3D/CitcomS/trunk/lib/Regional_lith_age_read_files.c 2007-03-12 20:08:12 UTC (rev 6229)
@@ -42,6 +42,7 @@
int i,j,node;
int intage, pos_age;
+
nox=E->mesh.nox;
noy=E->mesh.noy;
noz=E->mesh.noz;
@@ -51,6 +52,9 @@
lev=E->mesh.levmax;
age=find_age_in_MY(E);
+ if(E->parallel.me==0) {
+ fprintf(stderr,"INSIDE regional_lith_age_read_files, age= %g\n",age);
+ }
if (age < 0.0) { /* age is negative -> use age=0 for input files */
intage = 0;
Modified: mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c 2007-03-12 16:35:30 UTC (rev 6228)
+++ mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c 2007-03-12 20:08:12 UTC (rev 6229)
@@ -48,9 +48,20 @@
int nox,noz,noy,nnn,nox1,noz1,noy1,lev;
int i,ii,ll,mm,j,k,n,nodeg,nodel,node;
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();
+
+ if( E->parallel.me == 0) fprintf(stderr, "\nINSIDE regional_read_input_files_for_timesteps action=%d\n",action);
+
nox=E->mesh.nox;
noy=E->mesh.noy;
noz=E->mesh.noz;
@@ -59,6 +70,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 */
@@ -100,6 +117,31 @@
}
break;
+ case 3: /* read element materials */
+
+ sprintf(output_file1,"%s%0.0f.0",E->control.mat_file,newage1);
+ sprintf(output_file2,"%s%0.0f.0",E->control.mat_file,newage2);
+ 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");
+ }
+
} /* end switch */
@@ -147,6 +189,48 @@
free ((void *) VB2[i]);
}
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 (el=1; el<=elx*ely*elz; el++) {
+ nodea = E->ien[1][el].node[2];
+ llayer = layers(E,1,nodea);
+ if (llayer) { /* for layers:1-lithosphere,2-upper, 3-trans, and 4-lower mantle */
+ E->mat[1][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 (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[1][el] = VIP1[elg]+(VIP2[elg]-VIP1[elg])/(newage2-newage1)*(age-newage1);
+ E->mat[1][el] = LL1[elg];
+
+ } /* end for j */
+ } /* end for i */
+ } /* end for k */
+
+ free ((void *) VIP1);
+ free ((void *) VIP2);
+ free ((void *) LL1);
+ free ((void *) LL2);
+
} /* end switch */
return;
More information about the cig-commits
mailing list