[cig-commits] r7901 - in mc/3D/CitcomS/trunk: bin lib
becker at geodynamics.org
becker at geodynamics.org
Tue Aug 28 17:51:22 PDT 2007
Author: becker
Date: 2007-08-28 17:51:20 -0700 (Tue, 28 Aug 2007)
New Revision: 7901
Modified:
mc/3D/CitcomS/trunk/bin/Citcom.c
mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
mc/3D/CitcomS/trunk/lib/Initial_temperature.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Output_gzdir.c
mc/3D/CitcomS/trunk/lib/Output_h5.c
mc/3D/CitcomS/trunk/lib/Process_buoyancy.c
mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/convection_variables.h
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
When checking in my changes, I had to resolve a conflict for
lib/Instructions.c by hand, which I hope I did properly. Here are my changes:
- renamed CONTOL structure members ORTHO and ORTHOZ to CITCOM_ORTHO and CITCOM_ORTHOZ
Those were never used and conflicted with definitions in GMT gmt.h
- Added a higher frequency heat flow output option like so:
write_q_files=1 # option to write heat flux to files qt.dat and qb.dat
# at intervals smaller than storage_spacing (0)
- Added the capability to read in initial temperatures from netcdf grd files, if -USE_GGRD is used.
Lot of options, like so:
#
# read initial temperature conditions from grd files (default values in parentheses)
#
tic_method=4 # read initial temperature from netcdf GRD files (off)
ggrd_tinit_scale_with_prem=off # scale the temperature with PREM densities (off)
ggrd_tinit_scale=1.0 # scaling factor to apply to read in scalars f (1.0)
ggrd_tinit_offset=-0.5 # offset, T = f * scale + offset + tm (0.0)
# where tm is the mean between top and bottom TBC values
# if the bottom is flux, will use 1 for bottom TBC value
ggrd_tinit_gfile="../../data/tomography/s20a_smean_new_age/t" # prefix of grd files, will
# try to read t.1.grd, t.2.grd ... t.n.grd
# where n is the number of layers in the depth file
ggrd_tinit_dfile="../../data/tomography/s20a_smean_new_age/tdepth.dat" # file with layer depths in km from bottom up
ggrd_tinit_prem_file="../../progs/src/hc-svn/prem/prem.dat" # PREM data file
ggrd_tinit_override_tbc=on # override temperature boundary conditions (off)
Modified: mc/3D/CitcomS/trunk/bin/Citcom.c
===================================================================
--- mc/3D/CitcomS/trunk/bin/Citcom.c 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/bin/Citcom.c 2007-08-29 00:51:20 UTC (rev 7901)
@@ -169,6 +169,9 @@
tracer_advection(E);
general_stokes_solver(E);
+ if(E->output.write_q_files)
+ if ((E->monitor.solution_cycles % E->output.write_q_files)==0)
+ heat_flux(E);
if ((E->monitor.solution_cycles % E->control.record_every)==0) {
(E->problem_output)(E, E->monitor.solution_cycles);
Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-08-29 00:51:20 UTC (rev 7901)
@@ -29,6 +29,9 @@
#include "global_defs.h"
#include "parallel_related.h"
+#ifdef USE_GGRD
+void ggrd_full_temp_init(struct All_variables *);
+#endif
void get_r_spacing_fine(double *, double , double ,int ,double , double ,double , double,struct All_variables *);
@@ -270,8 +273,9 @@
gnoz=E->mesh.noz;
- if (E->convection.tic_method == 0) {
+ switch (E->convection.tic_method){
+ case 0:
/* set up a linear temperature profile first */
for(m=1;m<=E->sphere.caps_per_proc;m++)
@@ -311,8 +315,10 @@
E->T[m][node] += con*modified_plgndr_a(ll,mm,t1)*cos(mm*f1);
}
}
- }
- else if (E->convection.tic_method == 3) {
+ break;
+
+ case 3:
+
/* set up a linear temperature profile first */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=noy;i++)
@@ -354,11 +360,24 @@
*sin(M_PI*(r1-E->sphere.ri)/(E->sphere.ro-E->sphere.ri));
}
}
+ break;
+ case 1:
+ case 2:
+ break;
+ case 4:
+#ifdef USE_GGRD
+ ggrd_full_temp_init(E);
+#else
+ fprintf(stderr,"tic_method 4 only works for USE_GGRD compiled code\n");
+ parallel_process_termination();
+#endif
+ break;
+ default: /* unknown option */
+ fprintf(stderr,"Invalid value of 'tic_method'\n");
+ parallel_process_termination();
+ break;
}
- else if (E->convection.tic_method == 1) {
- }
-
temperatures_conform_bcs(E);
return;
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2007-08-29 00:51:20 UTC (rev 7901)
@@ -45,6 +45,8 @@
void ggrd_init_tracer_flavors(struct All_variables *);
int layers_r(struct All_variables *,float );
void myerror(struct All_variables *,char *);
+void ggrd_full_temp_init(struct All_variables *);
+void ggrd_reg_temp_init(struct All_variables *);
/*
@@ -135,21 +137,224 @@
}
+/*
+initialize temperatures from grd files for spherical geometry
+*/
+void ggrd_full_temp_init(struct All_variables *E)
+{
+
+ MPI_Status mpi_stat;
+ int mpi_rc;
+ int mpi_inmsg, mpi_success_message = 1;
+ double temp1,tbot,tgrad,tmean,tadd,rho_prem;
+ int i,j,k,m,node,noxnoz,nox,noy,noz;
+ noy=E->lmesh.noy;
+ nox=E->lmesh.nox;
+ noz=E->lmesh.noz;
+ noxnoz = nox * noz;
+ if(E->parallel.me == 0)
+ fprintf(stderr,"ggrd_full_temp_init: using GMT grd files for temperatures\n");
+ /*
+
+
+ read in tempeatures/density from GMT grd files
+
+
+ */
+ /*
+
+ begin MPI synchronization part
+
+ */
+ if(E->parallel.me > 0){
+ /*
+ wait for the previous processor
+ */
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1),
+ 0, MPI_COMM_WORLD, &mpi_stat);
+ }
+
+ if(E->convection.ggrd_tinit_scale_with_prem){/* initialize PREM */
+ if(prem_read_model(E->convection.prem.model_filename,
+ &E->convection.prem, (E->parallel.me==0)))
+ myerror(E,"PREM init error");
+ }
+ /*
+ initialize the GMT grid files
+ */
+ E->convection.ggrd_tinit_d[0].init = FALSE;
+ if(ggrd_grdtrack_init_general(TRUE,E->convection.ggrd_tinit_gfile,
+ E->convection.ggrd_tinit_dfile,"-L", /* this could also be -Lg for global */
+ E->convection.ggrd_tinit_d,(E->parallel.me == 0),
+ FALSE))
+ myerror(E,"grd init error");
+ if(E->parallel.me < E->parallel.nproc-1){
+ /* tell the next processor to go ahead with the init step */
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 0, MPI_COMM_WORLD);
+ }else{
+ fprintf(stderr,"ggrd_full_temp_init: last processor (%i) done with grd init\n",
+ E->parallel.me);
+ }
+ /*
+
+ interpolate densities to temperature given PREM variations
+
+ */
+ if(E->mesh.bottbc == 1){
+ /* bottom has specified temperature */
+ tbot = E->control.TBCbotval;
+ }else{
+ /*
+ bottom has specified heat flux start with unity bottom temperature
+ */
+ tbot = 1.0;
+ }
+ /*
+ mean temp is (top+bot)/2 + offset
+ */
+ tmean = (tbot + E->control.TBCtopval)/2.0 + E->convection.ggrd_tinit_offset;
+ 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 numbers */
+ node=k+(j-1)*noz+(i-1)*noxnoz;
+ /*
+ get interpolated velocity anomaly
+ */
+ if(!ggrd_grdtrack_interpolate_rtp((double)E->sx[m][3][node],
+ (double)E->sx[m][1][node],
+ (double)E->sx[m][2][node],
+ E->convection.ggrd_tinit_d,&tadd,
+ FALSE))
+ myerror(E,"ggrd_full_temp_init");
+ if(E->convection.ggrd_tinit_scale_with_prem){
+ /*
+ get the PREM density at r for additional scaling
+ */
+ prem_get_rho(&rho_prem,(double)E->sx[m][3][node],&E->convection.prem);
+ if(rho_prem < 3200.0)
+ rho_prem = 3200.0; /* we don't want the density of water */
+ /*
+ assign temperature
+ */
+ E->T[m][node] = tmean + tadd * E->convection.ggrd_tinit_scale *
+ rho_prem / E->data.density;
+ }else{
+ /* no PREM scaling */
+ E->T[m][node] = tmean + tadd * E->convection.ggrd_tinit_scale;
+ }
+ if(E->convection.ggrd_tinit_override_tbc){
+ if((k == 1) && (E->mesh.bottbc == 1)){ /* bottom TBC */
+ 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];
+ }
+ if((k == noz) && (E->mesh.toptbc == 1)){ /* top TBC */
+ 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];
+ }
+ }
+
+ //if(E->convection.ggrd_tinit_limit_unity)
+ /* limit to >0 */
+ //E->T[m][node] = min(max(E->T[m][node], 0.0),1.0);
+ }
+ /*
+ free the structure, not needed anymore since T should now
+ change internally
+ */
+ ggrd_grdtrack_free_gstruc(E->convection.ggrd_tinit_d);
+ /*
+ end temperature/density from GMT grd init
+ */
+}
+/*
+initialize temperatures from grd files for spherical geometry, regional version
+(this is identical to global version but for flags for GMT interpolation)
+*/
+void ggrd_reg_temp_init(struct All_variables *E)
+{
+ MPI_Status mpi_stat;
+ int mpi_rc, mpi_inmsg, mpi_success_message = 1;
+ double temp1,tbot,tgrad,tmean,tadd,rho_prem;
+ int i,j,k,m,ii,node,noxnoz,nox,noz,noy;
+ noy=E->lmesh.noy;
+ nox=E->lmesh.nox;
+ noz=E->lmesh.noz;
+ noxnoz = nox * noz;
+ if(E->parallel.me==0)
+ fprintf(stderr,"ggrd_reg_temp_init: using GMT grd files for temperatures\n");
+ if(E->parallel.me > 0)
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), 0, MPI_COMM_WORLD, &mpi_stat);
+ if(E->convection.ggrd_tinit_scale_with_prem){
+ if(prem_read_model(E->convection.prem.model_filename,&E->convection.prem, (E->parallel.me==0)))
+ myerror(E,"prem error");
+ }
+ if(ggrd_grdtrack_init_general(TRUE,E->convection.ggrd_tinit_gfile,
+ E->convection.ggrd_tinit_dfile,"", /* this part differnent from global */
+ E->convection.ggrd_tinit_d,(E->parallel.me==0),
+ FALSE))
+ myerror(E,"grdtrack init error");
+ if(E->parallel.me < E->parallel.nproc-1){
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 0, MPI_COMM_WORLD);
+ }else{
+ fprintf(stderr,"ggrd_reg_temp_init: last processor (%i) done with grd init\n",
+ E->parallel.me);
+ }
+ tmean = (E->control.TBCbotval + E->control.TBCtopval)/2.0 + E->convection.ggrd_tinit_offset;
+ 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)*noxnoz;
+ if(!ggrd_grdtrack_interpolate_rtp((double)E->sx[m][3][node],(double)E->sx[m][1][node],
+ (double)E->sx[m][2][node],E->convection.ggrd_tinit_d,&tadd,FALSE))
+ myerror(E,"ggrd_reg_temp_init");
+ if(E->convection.ggrd_tinit_scale_with_prem){
+ prem_get_rho(&rho_prem,(double)E->sx[m][3][node],&E->convection.prem);
+ if(rho_prem < 3200.0)
+ rho_prem = 3200.0;
+ E->T[m][node] = tmean + tadd * E->convection.ggrd_tinit_scale * rho_prem / E->data.density;
+ }else{
+ E->T[m][node] = tmean + tadd * E->convection.ggrd_tinit_scale;
+ }
+ if(E->convection.ggrd_tinit_override_tbc){
+ if((k == 1) && (E->mesh.bottbc == 1)){ /* bottom TBC */
+ 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];
+ }
+ if((k == noz) && (E->mesh.toptbc == 1)){ /* top TBC */
+ 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];
+ }
+ }
+ }
+ ggrd_grdtrack_free_gstruc(E->convection.ggrd_tinit_d);
+}
+
+
+
+
+
#endif
Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2007-08-29 00:51:20 UTC (rev 7901)
@@ -74,10 +74,15 @@
When tic_method is 3, the temperature is a linear profile + perturbation
for whole mantle.
+
+ tic_method is 4: read in initial temperature distribution from a set of netcdf grd
+ files. this required the GGRD extension to be compiled in
+
*/
-
- if (E->convection.tic_method == 0 || E->convection.tic_method == 3 ) {
+ switch(E->convection.tic_method){
+ case 0:
+ case 3:
/* 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. */
@@ -111,14 +116,14 @@
E->convection.perturb_ll[0] = 2;
E->convection.load_depth[0] = (noz+1)/2;
}
+
+ break;
+ case 1: /* case 1 */
- }
- else if (E->convection.tic_method == 1) {
-
input_float("half_space_age", &(E->convection.half_space_age), "40.0,1e-3,nomax", m);
+ break;
- }
- else if (E->convection.tic_method == 2) {
+ case 2: /* case 2 */
input_float("half_space_age", &(E->convection.half_space_age), "40.0,1e-3,nomax", m);
if( ! input_float_vector("blob_center", 3, E->convection.blob_center, m)) {
assert( E->sphere.caps == 12 || E->sphere.caps == 1 );
@@ -135,10 +140,38 @@
}
input_float("blob_radius", &(E->convection.blob_radius), "0.063,0.0,1.0", m);
input_float("blob_dT", &(E->convection.blob_dT), "0.18,nomin,nomax", m);
- }
- else {
+ break;
+ case 4:
+ /*
+ case 4: initial temp from grd files
+ */
+#ifdef USE_GGRD
+ /* read in some more parameters */
+ /* scale the anomalies with PREM densities */
+ input_boolean("ggrd_tinit_scale_with_prem",&(E->convection.ggrd_tinit_scale_with_prem),"off",E->parallel.me);
+ /* scaling factor for the grids */
+ input_double("ggrd_tinit_scale",&(E->convection.ggrd_tinit_scale),"1.0",E->parallel.me); /* scale */
+ /* temperature offset factor */
+ input_double("ggrd_tinit_offset",&(E->convection.ggrd_tinit_offset),"0.0",E->parallel.me); /* offset */
+ /* grid name, without the .i.grd suffix */
+ input_string("ggrd_tinit_gfile",E->convection.ggrd_tinit_gfile,"",E->parallel.me); /* grids */
+ input_string("ggrd_tinit_dfile",E->convection.ggrd_tinit_dfile,"",E->parallel.me); /* depth.dat
+ layers
+ of
+ grids*/
+ /* override temperature boundary condition? */
+ input_boolean("ggrd_tinit_override_tbc",&(E->convection.ggrd_tinit_override_tbc),"off",E->parallel.me);
+ input_string("ggrd_tinit_prem_file",E->convection.prem.model_filename,"", E->parallel.me); /* PREM model filename */
+#else
+ fprintf(stderr,"tic_method 4 only works for USE_GGRD compiled code\n");
+ parallel_process_termination();
+#endif
+
+ break;
+ default: /* unknown option */
fprintf(stderr,"Invalid value of 'tic_method'\n");
parallel_process_termination();
+ break;
}
return;
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-29 00:51:20 UTC (rev 7901)
@@ -78,7 +78,9 @@
void viscosity_input(struct All_variables*);
void get_vtk_filename(char *,int,struct All_variables *,int);
void myerror(struct All_variables *,char *);
+void open_qfiles(struct All_variables *) ;
+
void initial_mesh_solver_setup(struct All_variables *E)
{
@@ -404,7 +406,14 @@
input_int("storage_spacing",&(E->control.record_every),"10",m);
input_int("checkpointFrequency",&(E->control.checkpoint_frequency),"100",m);
input_int("cpu_limits_in_seconds",&(E->control.record_all_until),"5",m);
+ input_int("write_q_files",&(E->output.write_q_files),"0",m);/* write additional
+ heat flux files? */
+ if(E->output.write_q_files){ /* make sure those get written at
+ least as often as velocities */
+ E->output.write_q_files = min(E->output.write_q_files,E->control.record_every);
+ }
+
input_boolean("precond",&(E->control.precondition),"off",m);
input_int("mg_cycle",&(E->control.mg_cycle),"2,0,nomax",m);
input_int("down_heavy",&(E->control.down_heavy),"1,0,nomax",m);
@@ -776,8 +785,8 @@
/* SECOND: values for which an obvious default setting is useful */
- E->control.ORTHO = 1; /* for orthogonal meshes by default */
- E->control.ORTHOZ = 1; /* for orthogonal meshes by default */
+ E->control.CITCOM_ORTHO = 1; /* for orthogonal meshes by default */
+ E->control.CITCOM_ORTHOZ = 1; /* for orthogonal meshes by default */
E->control.stokes=0;
@@ -1061,7 +1070,32 @@
return;
}
+void open_qfiles(struct All_variables *E) /* additional heat
+ flux output */
+{
+ char output_file[255];
+ /* only one CPU will write to those */
+
+ /* top heat flux and other stat quantities */
+ if (strcmp(E->output.format, "ascii-gz") == 0)
+ sprintf(output_file,"%s/qt.dat", E->control.data_dir);
+ else
+ sprintf(output_file,"%s.qt.dat", E->control.data_file);
+ E->output.fpqt = output_open(output_file);
+ /* bottom heat flux and other stat quantities */
+ if (strcmp(E->output.format, "ascii-gz") == 0)
+ sprintf(output_file,"%s/qb.dat", E->control.data_dir);
+ else
+ sprintf(output_file,"%s.qb.dat", E->control.data_file);
+ E->output.fpqb = output_open(output_file);
+
+
+
+ return;
+}
+
+
static void output_parse_optional(struct All_variables *E)
{
char* strip(char*);
@@ -1261,6 +1295,11 @@
open_log(E);
open_time(E);
open_info(E);
+ if(E->output.write_q_files)
+ open_qfiles(E);
+ else{
+ E->output.fpqt = E->output.fpqb = NULL;
+ }
if (strcmp(E->output.format, "ascii") == 0) {
E->problem_output = output;
@@ -1305,7 +1344,13 @@
if (E->fp_out)
fclose(E->fp_out);
+
+ if(E->output.fpqt)
+ fclose(E->output.fpqt);
+ if(E->output.fpqb)
+ fclose(E->output.fpqb);
+
#ifdef USE_GZDIR
/*
remove VTK geo file in case we used that for IO
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2007-08-29 00:51:20 UTC (rev 7901)
@@ -230,7 +230,9 @@
FILE* fp2;
float *topo;
- heat_flux(E);
+ if(E->output.write_q_files == 0) /* else, the heat flux will have
+ been computed already */
+ heat_flux(E);
get_STD_topo(E,E->slice.tpg,E->slice.tpgb,E->slice.divg,E->slice.vort,cycles);
if (E->output.surf && (E->parallel.me_loc[3]==E->parallel.nprocz-1)) {
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-08-29 00:51:20 UTC (rev 7901)
@@ -739,8 +739,9 @@
char output_file[255];
gzFile *fp2;
float *topo;
-
- heat_flux(E);
+ if(E->output.write_q_files == 0) /* else, the heat flux will have
+ been computed already */
+ heat_flux(E);
get_STD_topo(E,E->slice.tpg,E->slice.tpgb,E->slice.divg,E->slice.vort,cycles);
if (E->output.surf && (E->parallel.me_loc[3]==E->parallel.nprocz-1)) {
@@ -1144,6 +1145,7 @@
/* E->sphere.cap[m].V[1][i] = v1;
E->sphere.cap[m].V[1][i] = v2;
E->sphere.cap[m].V[1][i] = v3; */
+ /* I don't like that */
//E->T[m][i] = max(0.0,min(g,1.0));
E->T[m][i] = g;
}
Modified: mc/3D/CitcomS/trunk/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_h5.c 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/lib/Output_h5.c 2007-08-29 00:51:20 UTC (rev 7901)
@@ -796,8 +796,11 @@
mx = scalar->block[1];
my = scalar->block[2];
+ if(E->output.write_q_files == 0) /* else, the heat flux will have
+ been computed already */
+ heat_flux(E);
- heat_flux(E);
+
get_STD_topo(E, E->slice.tpg, E->slice.tpgb, E->slice.divg, E->slice.vort, cycles);
Modified: mc/3D/CitcomS/trunk/lib/Process_buoyancy.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Process_buoyancy.c 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/lib/Process_buoyancy.c 2007-08-29 00:51:20 UTC (rev 7901)
@@ -167,6 +167,8 @@
if (E->parallel.me==E->parallel.nprocz-1) {
fprintf(stderr,"surface heat flux= %f\n",sum_h[0]);
/*fprintf(E->fp,"surface heat flux= %f\n",sum_h[0]);*/
+ if(E->output.write_q_files)
+ fprintf(E->output.fpqt,"%.5e %.5e\n",E->monitor.elapsed_time,sum_h[0]);
}
}
@@ -176,6 +178,9 @@
if (E->parallel.me==0) {
fprintf(stderr,"bottom heat flux= %f\n",sum_h[2]);
fprintf(E->fp,"bottom heat flux= %f\n",sum_h[2]);
+ if(E->output.write_q_files)
+ fprintf(E->output.fpqb,"%.5e %.5e\n",E->monitor.elapsed_time,sum_h[2]);
+
}
}
Modified: mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2007-08-29 00:51:20 UTC (rev 7901)
@@ -31,6 +31,10 @@
#include "parallel_related.h"
void get_r_spacing_fine(double *, double , double ,int ,double , double ,double , double,struct All_variables *);
+#ifdef USE_GGRD
+void ggrd_reg_temp_init(struct All_variables *);
+#endif
+
/* Setup global mesh parameters */
void regional_global_derived_values(E)
struct All_variables *E;
@@ -333,6 +337,14 @@
int mm, ll;
double con, temp;
+ double theta_center;
+ double fi_center;
+ double r_center;
+ double radius;
+ double amp;
+ double x_center,y_center,z_center;
+ double theta,fi,r,x,y,z,distance;
+
double tlen = M_PI / (E->control.theta_max - E->control.theta_min);
double flen = M_PI / (E->control.fi_max - E->control.fi_min);
@@ -341,9 +353,9 @@
noz=E->lmesh.noz;
gnoz=E->mesh.noz;
- if (E->convection.tic_method == 0) {
+ switch (E->convection.tic_method){
+ case 0:
-
/* set up a linear temperature profile first */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=noy;i++)
@@ -387,8 +399,9 @@
}
}
- }
- else if (E->convection.tic_method == 1) {
+ break;
+ case 1:
+
/* set up a top thermal boundary layer */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=noy;i++)
@@ -400,19 +413,19 @@
E->T[m][node] = E->control.TBCbotval*erf(temp);
}
- }
- else if (E->convection.tic_method == 2) {
- double temp;
- double theta_center = E->convection.blob_center[0];
- double fi_center = E->convection.blob_center[1];
- double r_center = E->convection.blob_center[2];
- double radius = E->convection.blob_radius;
- double amp = E->convection.blob_dT;
- double x_center,y_center,z_center;
- double theta,fi,r,x,y,z,distance;
+ break;
- fprintf(stderr,"center=%e %e %e radius=%e dT=%e\n",theta_center,fi_center,r_center,radius,amp);
- /* set up a thermal boundary layer first */
+ case 2:
+
+
+ theta_center = E->convection.blob_center[0];
+ fi_center = E->convection.blob_center[1];
+ r_center = E->convection.blob_center[2];
+ radius = E->convection.blob_radius;
+ amp = E->convection.blob_dT;
+
+ fprintf(stderr,"center=%e %e %e radius=%e dT=%e\n",theta_center,fi_center,r_center,radius,amp);
+ /* set up a thermal boundary layer first */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=noy;i++)
for(j=1;j<=nox;j++)
@@ -446,8 +459,9 @@
if (distance < radius)
E->T[m][node] += amp * exp(-1.0*distance/radius);
}
- }
- else if (E->convection.tic_method == 3) {
+ break;
+ case 3:
+
/* set up a linear temperature profile first */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=noy;i++)
@@ -484,7 +498,23 @@
*sin(M_PI*(r1-E->sphere.ri)/(E->sphere.ro-E->sphere.ri));
}
}
+ break;
+ case 4: /* from grd files */
+#ifdef USE_GGRD
+ ggrd_reg_temp_init(E);
+#else
+ fprintf(stderr,"tic_method 4 only works for USE_GGRD compiled code\n");
+ parallel_process_termination();
+#endif
+ break;
+
+ default: /* unknown option */
+ fprintf(stderr,"Invalid value of 'tic_method'\n");
+ parallel_process_termination();
+ break;
}
+
+
temperatures_conform_bcs(E);
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-29 00:51:20 UTC (rev 7901)
@@ -652,6 +652,7 @@
double vmean,cc_loc;
int m,l,z,jj,kk,i;
+
const int vpts = vpoints[E->mesh.nsd];
const int nel = E->lmesh.nel;
const int ends = enodes[E->mesh.nsd];
Modified: mc/3D/CitcomS/trunk/lib/convection_variables.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/convection_variables.h 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/lib/convection_variables.h 2007-08-29 00:51:20 UTC (rev 7901)
@@ -25,6 +25,8 @@
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
+
+
struct CONVECTION { /* information controlling convection problems */
int tic_method;
@@ -48,6 +50,19 @@
float lambda[10];
} heat_sources;
+
+#ifdef USE_GGRD
+ /* for temperature init from grd files */
+ int ggrd_tinit,ggrd_tinit_scale_with_prem;
+ int ggrd_tinit_override_tbc;
+ double ggrd_tinit_scale,ggrd_tinit_offset,ggrd_vstage_transition;
+ char ggrd_tinit_gfile[1000];
+ char ggrd_tinit_dfile[1000];
+ struct ggrd_gt ggrd_tinit_d[1];
+ struct ggrd_t ggrd_time_hist;
+ struct prem_model prem;
+#endif
+
} convection;
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-28 22:52:45 UTC (rev 7900)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-29 00:51:20 UTC (rev 7901)
@@ -39,6 +39,11 @@
#include <stdlib.h>
#include "mpi.h"
+#ifdef USE_GGRD
+#include "hc.h"
+#endif
+
+
#ifdef USE_HDF5
#include "hdf5.h"
#endif
@@ -49,6 +54,9 @@
#else
+
+
+
/* Macros */
#define max(A,B) (((A) > (B)) ? (A) : (B))
#define min(A,B) (((A) < (B)) ? (A) : (B))
@@ -382,11 +390,12 @@
struct CONTROL {
int PID;
-
+
char output_written_external_command[500]; /* a unix command to run when output files have been created */
+
+ /* this clashed with a GMT definition of ORTHO TWB */
+ int CITCOM_ORTHO,CITCOM_ORTHOZ; /* indicates levels of mesh symmetry */
- int ORTHO,ORTHOZ; /* indicates levels of mesh symmetry */
-
char data_prefix[50];
char data_prefix_old[50];
@@ -616,6 +625,10 @@
/* flags used by GZDIR */
struct gzd_struc gzdir;
+
+
+ int write_q_files;
+ FILE *fpqt,*fpqb; /* additional heat flux output */
};
More information about the cig-commits
mailing list