[cig-commits] r11888 - in mc/3D/CitcomS/trunk: bin lib
gurnis at geodynamics.org
gurnis at geodynamics.org
Thu May 1 15:56:12 PDT 2008
Author: gurnis
Date: 2008-05-01 15:56:12 -0700 (Thu, 01 May 2008)
New Revision: 11888
Modified:
mc/3D/CitcomS/trunk/bin/Citcom.c
mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c
mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c
mc/3D/CitcomS/trunk/lib/Problem_related.c
mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c
mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
Log:
More additions to allow top surface temperatures to
be read in a each time step. Should read in every
time step for C only version. Additional changes
for Pyre version to follow; will test more once
these changes are made.
Small changes to viscosity structures when reading
in materials for each element. More documentation to
follow.
Modified: mc/3D/CitcomS/trunk/bin/Citcom.c
===================================================================
--- mc/3D/CitcomS/trunk/bin/Citcom.c 2008-05-01 22:15:52 UTC (rev 11887)
+++ mc/3D/CitcomS/trunk/bin/Citcom.c 2008-05-01 22:56:12 UTC (rev 11888)
@@ -224,6 +224,9 @@
#endif
if(E->control.vbcs_file==1)
read_velocity_boundary_from_file(E);
+
+ if(E->control.tbcs_file)
+ read_temperature_boundary_from_file(E);
/*
else
renew_top_velocity_boundary(E);
Modified: mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c 2008-05-01 22:15:52 UTC (rev 11887)
+++ mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c 2008-05-01 22:56:12 UTC (rev 11888)
@@ -121,6 +121,8 @@
if(E->mesh.toptbc == 1) {
horizontal_bc(E,E->sphere.cap[j].TB,noz,3,E->control.TBCtopval,TBZ,1,lev,j);
horizontal_bc(E,E->sphere.cap[j].TB,noz,3,E->control.TBCtopval,FBZ,0,lev,j);
+ if(E->control.tbcs_file)
+ read_temperature_boundary_from_file(E);
}
else {
horizontal_bc(E,E->sphere.cap[j].TB,noz,3,E->control.TBCtopval,TBZ,0,lev,j);
Modified: mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c 2008-05-01 22:15:52 UTC (rev 11887)
+++ mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c 2008-05-01 22:56:12 UTC (rev 11888)
@@ -47,7 +47,7 @@
FILE *fp1, *fp2;
float age, newage1, newage2;
char output_file1[255],output_file2[255];
- float *VB1[4],*VB2[4], inputage1, inputage2;
+ float *TB1, *TB2, *VB1[4],*VB2[4], inputage1, inputage2;
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;
@@ -189,6 +189,31 @@
#endif
break;
/* mode 4 is rayleigh control for GGRD, see below */
+
+ case 5: /* read temperature boundary conditions, top surface */
+ sprintf(output_file1,"%s%0.0f.%d",E->control.temperature_boundary_file,newage1,cap);
+ sprintf(output_file2,"%s%0.0f.%d",E->control.temperature_boundary_file,newage2,cap);
+ fp1=fopen(output_file1,"r");
+ if (fp1 == NULL) {
+ fprintf(E->fp,"(Problem_related #10) Cannot open %s\n",output_file1);
+ exit(8);
+ }
+ if (pos_age) {
+ fp2=fopen(output_file2,"r");
+ if (fp2 == NULL) {
+ fprintf(E->fp,"(Problem_related #11) Cannot open %s\n",output_file2);
+ exit(8);
+ }
+ }
+ if((E->parallel.me==0) && (output==1)) {
+ fprintf(E->fp,"Surface Temperature: Starting Age = %g, Elapsed time = %g, Current Age = %g\n",E->control.start_age,E->monitor.elapsed_time,age);
+ fprintf(E->fp,"Surface Temperature: File1 = %s\n",output_file1);
+ if (pos_age)
+ fprintf(E->fp,"Surface Temperature: File2 = %s\n",output_file2);
+ else
+ fprintf(E->fp,"Velocity: File2 = No file inputted (negative age)\n");
+ }
+ break;
} /* end switch */
@@ -285,13 +310,13 @@
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;
+ fprintf(stderr,"\nINSIDE llayer=&d",llayer);
}
}
for(i=1;i<=emax;i++) {
@@ -310,7 +335,7 @@
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];
+ /* E->mat[m][el] = LL1[elg]; */ /*get material numbers from radius internally */
} /* end for j */
} /* end for i */
@@ -334,6 +359,41 @@
#endif
break;
+ case 5: /* read temperature boundary conditions, top surface */
+ nnn=nox*noy;
+ TB1=(float*) malloc ((nnn+1)*sizeof(float));
+ TB2=(float*) malloc ((nnn+1)*sizeof(float));
+
+ for(i=1;i<=nnn;i++) {
+ fscanf(fp1,"%f",&(TB1[i]));
+ if (pos_age) {
+ fscanf(fp2,"%f",&(TB2[i]));
+ }
+ }
+ fclose(fp1);
+ if (pos_age) fclose(fp2);
+
+ if(E->parallel.me_loc[3]==E->parallel.nprocz-1 ) {
+ for(k=1;k<=noy1;k++)
+ for(i=1;i<=nox1;i++) {
+ nodeg = E->lmesh.nxs+i-1 + (E->lmesh.nys+k-2)*nox;
+ nodel = (k-1)*nox1*noz1 + (i-1)*noz1+noz1;
+ if (pos_age) { /* positive ages - we must interpolate */
+ E->sphere.cap[m].TB[1][nodel] = (TB1[nodeg] + (TB2[nodeg]-TB1[nodeg])/(newage2-newage1)*(age-newage1));
+ E->sphere.cap[m].TB[2][nodel] = (TB1[nodeg] + (TB2[nodeg]-TB1[nodeg])/(newage2-newage1)*(age-newage1));
+ E->sphere.cap[m].TB[3][nodel] = (TB1[nodeg] + (TB2[nodeg]-TB1[nodeg])/(newage2-newage1)*(age-newage1));
+ }
+ else { /* negative ages - don't do the interpolation */
+ E->sphere.cap[m].TB[1][nodel] = TB1[nodeg];
+ E->sphere.cap[m].TB[2][nodel] = TB1[nodeg];
+ E->sphere.cap[m].TB[3][nodel] = TB1[nodeg];
+ }
+ }
+ } /* end of E->parallel.me_loc[3]==E->parallel.nproczl-1 */
+ free ((void *) TB1);
+ free ((void *) TB2);
+ break;
+
} /* end switch */
} /* end for m */
Modified: mc/3D/CitcomS/trunk/lib/Problem_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Problem_related.c 2008-05-01 22:15:52 UTC (rev 11887)
+++ mc/3D/CitcomS/trunk/lib/Problem_related.c 2008-05-01 22:56:12 UTC (rev 11888)
@@ -63,8 +63,18 @@
return;
}
+/*=======================================================================
+ read temperature at the top surface from files
+=========================================================================*/
+void read_temperature_boundary_from_file(E)
+ struct All_variables *E;
+{
+ (E->solver.read_input_files_for_timesteps)(E,5,1); /* read temperature(5) and output(1) */
+ return;
+}
+
/*=======================================================================
Open restart file to get initial elapsed time, or calculate the right value
=========================================================================*/
Modified: mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c 2008-05-01 22:15:52 UTC (rev 11887)
+++ mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c 2008-05-01 22:56:12 UTC (rev 11888)
@@ -143,6 +143,9 @@
if(E->mesh.toptbc == 1) {
horizontal_bc(E,E->sphere.cap[j].TB,noz,3,E->control.TBCtopval,TBZ,1,lev,j);
horizontal_bc(E,E->sphere.cap[j].TB,noz,3,E->control.TBCtopval,FBZ,0,lev,j);
+ if(E->control.tbcs_file) {
+ read_temperature_boundary_from_file(E); /* read in the temperature boundary condition from file */
+ }
}
else {
horizontal_bc(E,E->sphere.cap[j].TB,noz,3,E->control.TBCtopval,TBZ,0,lev,j);
Modified: mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c 2008-05-01 22:15:52 UTC (rev 11887)
+++ mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c 2008-05-01 22:56:12 UTC (rev 11888)
@@ -47,7 +47,7 @@
FILE *fp1, *fp2;
float age, newage1, newage2;
char output_file1[255],output_file2[255];
- float *VB1[4],*VB2[4], inputage1, inputage2;
+ float *TB1, *TB2, *VB1[4],*VB2[4], inputage1, inputage2;
int nox,noz,noy,nnn,nox1,noz1,noy1,lev;
int i,ii,ll,mm,j,k,n,nodeg,nodel,node;
int intage, pos_age;
@@ -63,7 +63,8 @@
int llayer;
int layers();
- /* if( E->parallel.me == 0) fprintf(stderr, "\nINSIDE regional_read_input_files_for_timesteps action=%d\n",action); */
+ /*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;
@@ -188,7 +189,34 @@
}
#endif
break;
+
/* mode 4 is rayleigh control for GGRD, see below */
+
+ case 5: /* read temperature boundary conditions, top surface */
+ sprintf(output_file1,"%s%0.0f",E->control.temperature_boundary_file,newage1);
+ sprintf(output_file2,"%s%0.0f",E->control.temperature_boundary_file,newage2);
+ fp1=fopen(output_file1,"r");
+ if (fp1 == NULL) {
+ fprintf(E->fp,"(Problem_related #10) Cannot open %s\n",output_file1);
+ exit(8);
+ }
+ if (pos_age) {
+ fp2=fopen(output_file2,"r");
+ if (fp2 == NULL) {
+ fprintf(E->fp,"(Problem_related #11) Cannot open %s\n",output_file2);
+ exit(8);
+ }
+ }
+ if((E->parallel.me==0) && (output==1)) {
+ fprintf(E->fp,"Surface Temperature: Starting Age = %g, Elapsed time = %g, Current Age = %g\n",E->control.start_age,E->monitor.elapsed_time,age);
+ fprintf(E->fp,"Surface Temperature: File1 = %s\n",output_file1);
+ if (pos_age)
+ fprintf(E->fp,"Surface Temperature: File2 = %s\n",output_file2);
+ else
+ fprintf(E->fp,"Surface Temperature: File2 = No file inputted (negative age)\n");
+ }
+ break;
+
} /* end switch */
@@ -286,7 +314,6 @@
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);
@@ -309,7 +336,7 @@
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];
+ /* E->mat[1][el] = LL1[elg]; */ /*use the mat numbers base on radius*/
} /* end for j */
} /* end for i */
@@ -332,6 +359,43 @@
#endif
break;
+ case 5: /* read temperature boundary conditions, top surface */
+ nnn=nox*noy;
+ TB1=(float*) malloc ((nnn+1)*sizeof(float));
+ TB2=(float*) malloc ((nnn+1)*sizeof(float));
+
+ for(i=1;i<=nnn;i++) {
+ fscanf(fp1,"%f",&(TB1[i]));
+ if (pos_age) {
+ fscanf(fp2,"%f",&(TB2[i]));
+ }
+ }
+ fclose(fp1);
+ if (pos_age) fclose(fp2);
+
+ if(E->parallel.me_loc[3]==E->parallel.nprocz-1 ) {
+ for(k=1;k<=noy1;k++)
+ for(i=1;i<=nox1;i++) {
+ nodeg = E->lmesh.nxs+i-1 + (E->lmesh.nys+k-2)*nox;
+ nodel = (k-1)*nox1*noz1 + (i-1)*noz1+noz1;
+ if (pos_age) { /* positive ages - we must interpolate */
+ E->sphere.cap[1].TB[1][nodel] = (TB1[nodeg] + (TB2[nodeg]-TB1[nodeg])/(newage2-newage1)*(age-newage1));
+ E->sphere.cap[1].TB[2][nodel] = (TB1[nodeg] + (TB2[nodeg]-TB1[nodeg])/(newage2-newage1)*(age-newage1));
+ E->sphere.cap[1].TB[2][nodel] = (TB1[nodeg] + (TB2[nodeg]-TB1[nodeg])/(newage2-newage1)*(age-newage1));
+ }
+ else { /* negative ages - don't do the interpolation */
+ E->sphere.cap[1].TB[1][nodel] = TB1[nodeg];
+ E->sphere.cap[1].TB[2][nodel] = TB1[nodeg];
+ E->sphere.cap[1].TB[3][nodel] = TB1[nodeg];
+ }
+ }
+ } /* end of E->parallel.me_loc[3]==E->parallel.nprocz-1 */
+ free ((void *) TB1);
+ free ((void *) TB2);
+
+ break;
+
+
} /* end switch */
return;
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2008-05-01 22:15:52 UTC (rev 11887)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2008-05-01 22:56:12 UTC (rev 11888)
@@ -420,51 +420,48 @@
break;
- case 5: /* this still needs to be documented, who wrote this? */
+ case 5:
/* same as rheol 3, except alternative margin, VIP, formulation */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nel;i++) {
l = E->mat[m][i];
-
tempa = E->viscosity.N0[l-1];
+ /* fprintf(stderr,"\nINSIDE visc_from_T, l=%d, tempa=%g",l,tempa);*/
j = 0;
for(kk=1;kk<=ends;kk++) {
TT[kk] = E->T[m][E->ien[m][i].node[kk]];
- zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
+ /* zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]); */
}
for(jj=1;jj<=vpts;jj++) {
temp=0.0;
- zzz=0.0;
+ /* zzz=0.0; */
for(kk=1;kk<=ends;kk++) {
TT[kk]=max(TT[kk],zero);
temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
- zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+ /* zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)]; */
}
- if(E->control.mat_control==0){
+ if(E->control.mat_control==0)
EEta[m][ (i-1)*vpts + jj ] = tempa*
exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
- E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
- }else{
- /* 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]) );
- visc1 = tempa*E->viscosity.max_value;
- if(tempa_exp > visc1) tempa_exp=visc1;
- EEta[m][ (i-1)*vpts + jj ] = visc2*tempa_exp;
- */
- visc1 = E->VIP[m][i];
+ if(E->control.mat_control==1) {
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;
+ if(E->viscosity.MAX) {
+ if(visc2 > E->viscosity.max_value)
+ visc2 = E->viscosity.max_value;
+ }
+ if(E->viscosity.MIN) {
+ if(visc2 < E->viscosity.min_value)
+ visc2 = E->viscosity.min_value;
+ }
+ EEta[m][ (i-1)*vpts + jj ] = E->VIP[m][i]*visc2;
}
}
More information about the cig-commits
mailing list