https://github.com/geodynamics/citcoms
Tip revision: 74e1fe0aaa5811939f0483e545b9712885ca4722 authored by Eh Tan on 15 May 2003, 07:07:54 UTC
first created
first created
Tip revision: 74e1fe0
Initial_temperature.c
#include <math.h>
#include <sys/types.h>
#include "element_definitions.h"
#include "global_defs.h"
#include <stdlib.h> /* for "system" command */
#include <strings.h>
void convection_initial_temperature(E)
struct All_variables *E;
{
float find_age_in_MY();
int i,j,k,p,node,ii,jj,m,mm,ll;
int nox,noy,noz,noz2;
int gnox,gnoy,gnoz,nodeg;
char output_file[255],input_s[1000];
float lscale,ficenter,rcenter,dist2,lscalex,rscale,age;
double con,temp,t1,f1,r1;
float rad,beta,notusedhere;
float a,b,c,d,e,f,g;
FILE *fp,*fp1,*fp2;
float v1,v2,v3;
float e_4;
void temperatures_conform_bcs();
void thermal_buoyancy();
void parallel_process_termination();
void sphere_harmonics_layer();
void inv_sphere_harmonics();
const int dims=E->mesh.nsd;
rad = 180.0/M_PI;
e_4=1.e-4;
noy=E->lmesh.noy;
nox=E->lmesh.nox;
noz=E->lmesh.noz;
noz2=(E->mesh.noz-1)/2+1;
gnox=E->mesh.nox;
gnoy=E->mesh.noy;
gnoz=E->mesh.noz;
if ((E->control.restart || E->control.post_p)) {
/* used if restarting from a previous run. CPC 1/28/00 */
/*
E->monitor.solution_cycles=(E->control.restart)?E->control.restart:E->advection.max_timesteps;
*/
ii = E->monitor.solution_cycles_init;
sprintf(output_file,"%s.velo.%d.%d",E->control.old_P_file,E->parallel.me,ii);
fp=fopen(output_file,"r");
if (fp == NULL) {
fprintf(E->fp,"(Initial_temperature.c #1) Cannot open %s\n",output_file);
exit(8);
}
fgets(input_s,1000,fp);
sscanf(input_s,"%d %d %f",&ll,&mm,¬usedhere);
for(m=1;m<=E->sphere.caps_per_proc;m++) {
fgets(input_s,1000,fp);
sscanf(input_s,"%d %d",&ll,&mm);
for(i=1;i<=E->lmesh.nno;i++) {
fgets(input_s,1000,fp);
sscanf(input_s,"%g %g %g %f",&(v1),&(v2),&(v3),&(g));
/* E->sphere.cap[m].V[1][i] = d;
E->sphere.cap[m].V[1][i] = e;
E->sphere.cap[m].V[1][i] = f; */
E->T[m][i] = max(0.0,min(g,1.0));
}
}
fclose (fp);
/* use this section to include imposed ages at the surface CPC 4/27/00 */
if(E->control.lith_age==1) {
if(E->control.lith_age_time==1) {
/* if opening lithosphere age info every timestep - naming is different*/
age=find_age_in_MY(E);
sprintf(output_file,"%s%0.0f",E->control.lith_age_file,age);
if(E->parallel.me==0) {
fprintf(E->fp,"%s %s\n","Initial Lithosphere age info:",output_file);
}
}
else { /* just open lithosphere age info here*/
sprintf(output_file,"%s",E->control.lith_age_file);
}
fp1=fopen(output_file,"r");
if (fp1 == NULL) {
fprintf(E->fp,"(Convection.c #2) Cannot open %s\n",output_file);
exit(8);
}
for(i=1;i<=gnoy;i++)
for(j=1;j<=gnox;j++) {
node=j+(i-1)*gnox;
fscanf(fp1,"%f",&(E->age_t[node]));
E->age_t[node]=E->age_t[node]/E->data.scalet;
}
fclose(fp1);
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=noy;i++)
for(j=1;j<=nox;j++)
for(k=1;k<=noz;k++) {
nodeg=E->lmesh.nxs-1+j+(E->lmesh.nys+i-2)*gnox;
node=k+(j-1)*noz+(i-1)*nox*noz;
r1=E->sx[m][3][node];
if( r1 >= E->sphere.ro-E->control.lith_age_depth )
{ /* if closer than (lith_age_depth) from top */
temp = (E->sphere.ro-r1) *0.5 /sqrt(E->age_t[nodeg]);
E->T[m][node] = E->control.mantle_temp * erf(temp);
}
} /* end k */
} /* end lith age */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=noy;i++)
for(j=1;j<=nox;j++)
for(k=1;k<=noz;k++) {
node=k+(j-1)*noz+(i-1)*nox*noz;
r1=E->sx[m][3][node];
if(fabs(r1-E->sphere.ro)>=e_4 && fabs(r1-E->sphere.ri)>=e_4) {
E->sphere.cap[m].TB[1][node]=E->T[m][node];
E->sphere.cap[m].TB[2][node]=E->T[m][node];
E->sphere.cap[m].TB[3][node]=E->T[m][node];
}
}
}
else { /* if control.restart=0 => not using a restart file */
if(E->control.lith_age==1) {
/* used if the lithosphere age is given in an input file. CPC 1/28/00 */
gnox=E->mesh.nox;
gnoy=E->mesh.noy;
gnoz=E->mesh.noz;
if(E->control.lith_age_time==1) {
/* if opening lithosphere age info every timestep - naming is different*/
age=find_age_in_MY(E);
sprintf(output_file,"%s%0.0f",E->control.lith_age_file,age);
if(E->parallel.me==0) {
fprintf(E->fp,"%s %s\n","Initial Lithosphere age info:",output_file);
}
}
else { /* just open lithosphere age info here*/
sprintf(output_file,"%s",E->control.lith_age_file);
}
fp1=fopen(output_file,"r");
if (fp1 == NULL) {
fprintf(E->fp,"(Convection.c #3) Cannot open %s\n",output_file);
exit(8);
}
for(i=1;i<=gnoy;i++)
for(j=1;j<=gnox;j++) {
node=j+(i-1)*gnox;
fscanf(fp1,"%f",&(E->age_t[node]));
E->age_t[node]=E->age_t[node]/E->data.scalet;
}
fclose(fp1);
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=noy;i++)
for(j=1;j<=nox;j++)
for(k=1;k<=noz;k++) {
nodeg=E->lmesh.nxs-1+j+(E->lmesh.nys+i-2)*gnox;
node=k+(j-1)*noz+(i-1)*nox*noz;
E->T[m][node] = 1.0;
t1=E->sx[m][1][node];
f1=E->sx[m][2][node];
r1=E->sx[m][3][node];
temp = (E->sphere.ro-r1) *0.5 /sqrt(E->age_t[nodeg]);
E->T[m][node] = E->control.mantle_temp * erf(temp);
if(fabs(r1-E->sphere.ro)>=e_4 && fabs(r1-E->sphere.ri)>=e_4) {
E->sphere.cap[m].TB[1][node]=E->T[m][node];
E->sphere.cap[m].TB[2][node]=E->T[m][node];
E->sphere.cap[m].TB[3][node]=E->T[m][node];
}
} /* end k */
} /* end of lith-age */
else
/* used otherwise - Set your own initial temperature initial conditions! */
{
int number_of_perturbations;
int perturb_ll[32], perturb_mm[32], load_depth[32];
float perturb_mag[32];
float tlen, flen;
m = E->parallel.me;
tlen = M_PI / (E->control.theta_max - E->control.theta_min);
flen = M_PI / (E->control.fi_max - E->control.fi_min);
/* This part put a temperature anomaly at depth where the global
node number is equal to load_depth. The horizontal pattern of
the anomaly is given by spherical harmonic ll & mm. */
input_int("num_perturbations",&number_of_perturbations,"0,0,32",m);
if (number_of_perturbations > 0) {
if (! input_float_vector("perturbmag",number_of_perturbations,perturb_mag,m) ) {
fprintf(stderr,"Missing input parameter: 'perturbmag'\n");
parallel_process_termination();
}
if (! input_int_vector("perturbm",number_of_perturbations,perturb_mm,m) ) {
fprintf(stderr,"Missing input parameter: 'perturbm'\n");
parallel_process_termination();
}
if (! input_int_vector("perturbl",number_of_perturbations,perturb_ll,m) ) {;
fprintf(stderr,"Missing input parameter: 'perturbml'\n");
parallel_process_termination();
}
/* if (! input_int_vector("perturblayer",number_of_perturbations,load_depth,m) ) { */
/* fprintf(stderr,"Missing input parameter: 'perturblayer'\n"); */
/* parallel_process_termination(); */
/* } */
}
else {
number_of_perturbations = 1;
perturb_mag[0] = E->mesh.elz/(E->sphere.ro-E->sphere.ri);
perturb_mm[0] = 2;
perturb_ll[0] = 2;
/* load_depth[0] = noz/2; */
}
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=noy;i++)
for(j=1;j<=nox;j++)
for(k=1;k<=noz;k++) {
ii = k + E->lmesh.nzs - 1;
node=k+(j-1)*noz+(i-1)*nox*noz;
t1 = (E->sx[m][1][node] - E->control.theta_min) * tlen;
f1 = (E->sx[m][2][node] - E->control.fi_min) * flen;
r1 = E->sx[m][3][node] - E->sphere.ri;
E->T[m][node] = E->control.TBCbotval - (E->control.TBCtopval + E->control.TBCbotval)*r1/(E->sphere.ro - E->sphere.ri);
for (p=0; p<number_of_perturbations; p++) {
mm = perturb_mm[p];
ll = perturb_ll[p];
con = perturb_mag[p];
E->T[m][node] += con*cos(ll*f1)*cos(mm*t1)*sin(M_PI*r1/(E->sphere.ro - E->sphere.ri));
E->T[m][node] = max(min(E->T[m][node], 1.0), 0.0);
}
}
} /* end for else */
} /* end for else */
temperatures_conform_bcs(E);
if (E->control.verbose==1) {
fprintf(E->fp_out,"output_temperature\n");
for(m=1;m<=E->sphere.caps_per_proc;m++) {
fprintf(E->fp_out,"for cap %d\n",E->sphere.capid[m]);
for (j=1;j<=E->lmesh.nno;j++)
fprintf(E->fp_out,"X = %.6e Z = %.6e Y = %.6e T[%06d] = %.6e \n",E->sx[m][1][j],E->sx[m][2][j],E->sx[m][3][j],j,E->T[m][j]);
}
fflush(E->fp_out);
}
return;
}