[cig-commits] r6206 - in mc/3D/CitcomS/trunk: . CitcomS/Components
lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Thu Mar 8 16:25:04 PST 2007
Author: tan2
Date: 2007-03-08 16:25:03 -0800 (Thu, 08 Mar 2007)
New Revision: 6206
Added:
mc/3D/CitcomS/trunk/lib/Composition_related.c
mc/3D/CitcomS/trunk/lib/composition_related.h
Modified:
mc/3D/CitcomS/trunk/AUTHORS
mc/3D/CitcomS/trunk/CitcomS/Components/Output.py
mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py
mc/3D/CitcomS/trunk/lib/Element_calculations.c
mc/3D/CitcomS/trunk/lib/Full_obsolete.c
mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
mc/3D/CitcomS/trunk/lib/Global_operations.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Makefile.am
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/lib/tracer_defs.h
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
A big patch for tracers and composition:
* Seperating composition codes from tracer codes, several struct members in E->trace are moved to E->composition
* Add chemical buoyancy term
* Renamed thermal_buoyancy() to get_buoyancy() to reflect the fact that this function computes buoyancy according to temperature/composition/phase
* Fix and optimize the fix_*() functions
* Remove some static variables, or converte them to const
* Change the type of E->SinCos from float to double, and replace E->trace.DSinCos
* Provide default values for "required" input parameters
* Renamed tracer_restart to tracer_ic_method, also remapped the options
* Removed option for Cartesian tracer input and trace.iwrite_tracers_every
* Replaced trace.itracer_type by composition.ichemical_buoyancy
* Optionally output tracers and composition fields
Modified: mc/3D/CitcomS/trunk/AUTHORS
===================================================================
--- mc/3D/CitcomS/trunk/AUTHORS 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/AUTHORS 2007-03-09 00:25:03 UTC (rev 6206)
@@ -11,6 +11,7 @@
Clint Conrad
Michael Gurnis
Eun-seo Choi
+ Allen McNamara
Also see the file ChangeLog.
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Output.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Output.py 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Output.py 2007-03-09 00:25:03 UTC (rev 6206)
@@ -51,7 +51,8 @@
output_format = inv.str("output_format", default="ascii",
validator=inv.choice(["ascii",
"hdf5"]))
- output_optional = inv.str("output_optional", default="surf,botm")
+ output_optional = inv.str("output_optional",
+ default="surf,botm,tracer,comp_el")
# max. degree for spherical harmonics output
output_ll_max = inv.int("output_ll_max", default=20)
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py 2007-03-09 00:25:03 UTC (rev 6206)
@@ -60,24 +60,21 @@
tracer = inv.bool("tracer", default=False)
- # tracer_restart=-1 (read from file) TODO: rename to tracer_init_method
- # tracer_restart=0 (generate array)
- # tracer_restart=1 (read from scratch disks)
- tracer_restart = inv.int("tracer_restart", default=0)
+ # tracer_ic_method=0 (random generated array)
+ # tracer_ic_method=1 (all proc read the same file)
+ # tracer_ic_method=2 (each proc reads its restart file)
+ tracer_ic_method = inv.int("tracer_ic_method", default=0)
- # (tracer_restart == 0)
+ # (tracer_ic_method == 0)
tracers_per_element = inv.int("tracers_per_element", default=10)
- # TODO: remove
- z_interface = inv.float("z_interface", default=0.5)
- # (tracer_restart == -1 or 1)
+ # (tracer_ic_method == 1)
tracer_file = inv.str("tracer_file", default="tracer.dat")
- # icartesian_or_spherical=0 (cartesian coordinate input) */
- # icartesian_or_spherical=1 (spherical coordinate input) */
- cartesian_or_spherical_input = inv.int("cartesian_or_spherical_input",
- default=1)
+ # TODO: remove
+ z_interface = inv.float("z_interface", default=0.5)
+
# Advection Scheme
# itracer_advection_scheme=1
@@ -99,9 +96,7 @@
# Analytical Test Function
analytical_tracer_test = inv.int("analytical_tracer_test", default=0)
- # itracer_type=0 passive
- # itracer_type=1 active
- tracer_type = inv.int("tracer_type", default=1)
+ chemical_buoyancy = inv.int("chemical_buoyancy", default=1)
# ibuoy_type=0 (absolute method, not implemented)
# ibuoy_type=1 (ratio method)
@@ -116,11 +111,8 @@
compositional_prefactor = inv.float("compositional_prefactor",
default=1.0)
- # Output frequency TODO: remove
- write_tracers_every = inv.int("write_tracers_every", default=1000000)
-
# version
__id__ = "$Id$"
Added: mc/3D/CitcomS/trunk/lib/Composition_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Composition_related.c 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/Composition_related.c 2007-03-09 00:25:03 UTC (rev 6206)
@@ -0,0 +1,660 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *<LicenseText>
+ *
+ * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
+ * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
+ * Copyright (C) 1994-2005, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+#include <math.h>
+#include "global_defs.h"
+#include "parsing.h"
+#include "parallel_related.h"
+#include "composition_related.h"
+
+
+static void initialize_old_composition(struct All_variables *E);
+
+
+void composition_input(struct All_variables *E)
+{
+ int m = E->parallel.me;
+
+ input_int("chemical_buoyancy",&(E->composition.ichemical_buoyancy),
+ "1,0,nomax",m);
+
+ if (E->composition.ichemical_buoyancy==1) {
+
+ input_double("buoyancy_ratio",
+ &(E->composition.buoyancy_ratio),"1.0",m);
+
+ /* ibuoy_type=0 (absolute method) */
+ /* ibuoy_type=1 (ratio method) */
+
+ input_int("buoy_type",&(E->composition.ibuoy_type),"1,0,nomax",m);
+ if (E->composition.ibuoy_type!=1) {
+ fprintf(stderr,"Terror-Sorry, only ratio method allowed now\n");
+ fflush(stderr);
+ parallel_process_termination();
+ }
+
+ if (E->trace.ic_method==2) {
+ input_int("reset_initial_composition",
+ &(E->composition.ireset_initial_composition),"0",m);
+ }
+
+ input_double("z_interface",&(E->composition.z_interface),"0.5",m);
+
+ }
+
+
+ /* compositional rheology */
+
+ /* icompositional_rheology=0 (off) */
+ /* icompositional_rheology=1 (on) */
+
+ input_int("compositional_rheology",
+ &(E->composition.icompositional_rheology),"1,0,nomax",m);
+
+ if (E->composition.icompositional_rheology==1) {
+ input_double("compositional_prefactor",
+ &(E->composition.compositional_rheology_prefactor),
+ "1.0",m);
+ }
+
+}
+
+
+
+void composition_setup(struct All_variables *E)
+{
+
+ E->composition.on = 0;
+ if (E->composition.ichemical_buoyancy==1 ||
+ E->composition.icompositional_rheology)
+ E->composition.on = 1;
+
+ if (E->composition.ichemical_buoyancy==1 && E->composition.ibuoy_type==1) {
+ initialize_old_composition(E);
+ fill_composition(E);
+ }
+
+}
+
+
+
+void init_tracer_composition(struct All_variables *E)
+{
+ int j, kk, number_of_tracers;
+ double rad;
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ number_of_tracers = E->trace.itrac[j][0];
+ for (kk=1;kk<=number_of_tracers;kk++) {
+ rad = E->trace.rtrac[j][2][kk];
+
+ if (rad<=E->composition.z_interface) E->trace.etrac[j][0][kk]=1.0;
+ if (rad>E->composition.z_interface) E->trace.etrac[j][0][kk]=0.0;
+ }
+ }
+ return;
+}
+
+
+void write_composition_instructions(struct All_variables *E)
+{
+
+ if (E->composition.ichemical_buoyancy==0)
+ fprintf(E->trace.fpt,"Passive Tracers\n");
+
+ if (E->composition.ichemical_buoyancy==1)
+ fprintf(E->trace.fpt,"Active Tracers\n");
+
+
+ if (E->composition.ibuoy_type==1) fprintf(E->trace.fpt,"Ratio Method\n");
+ if (E->composition.ibuoy_type==0) fprintf(E->trace.fpt,"Absolute Method\n");
+
+ fprintf(E->trace.fpt,"Buoyancy Ratio: %f\n", E->composition.buoyancy_ratio);
+
+
+ if (E->composition.ireset_initial_composition==0)
+ fprintf(E->trace.fpt,"Using old initial composition from tracer files\n");
+ else
+ fprintf(E->trace.fpt,"Resetting initial composition\n");
+
+
+
+ if (E->composition.icompositional_rheology==0) {
+ fprintf(E->trace.fpt,"Compositional Rheology - OFF\n");
+ }
+ else if (E->composition.icompositional_rheology>0) {
+ fprintf(E->trace.fpt,"Compositional Rheology - ON\n");
+ fprintf(E->trace.fpt,"Compositional Prefactor: %f\n",
+ E->composition.compositional_rheology_prefactor);
+ }
+
+
+ fflush(E->trace.fpt);
+ fflush(stderr);
+
+}
+
+
+/************ FILL COMPOSITION ************************/
+void fill_composition(struct All_variables *E)
+{
+
+ void compute_elemental_composition_ratio_method();
+ void map_composition_to_nodes();
+
+ /* Currently, only the ratio method works here. */
+ /* Will have to come back here to include the absolute method. */
+
+ /* ratio method */
+
+ if (E->composition.ibuoy_type==1)
+ {
+ compute_elemental_composition_ratio_method(E);
+ }
+
+ /* absolute method */
+
+ if (E->composition.ibuoy_type!=1)
+ {
+ fprintf(E->trace.fpt,"Error(compute...)-only ratio method now\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ /* Map elemental composition to nodal points */
+
+ map_composition_to_nodes(E);
+
+ return;
+}
+
+
+
+/************ INITIALIZE OLD COMPOSITION ************************/
+static void initialize_old_composition(struct All_variables *E)
+{
+
+ char output_file[200];
+ char input_s[1000];
+
+ int ibottom_node;
+ int kk;
+ int j;
+
+ double zbottom;
+ double time;
+
+ FILE *fp;
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ if ((E->trace.oldel[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(fill old composition)-no memory 324c\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+
+ if ((E->trace.ic_method==0)||(E->trace.ic_method==1))
+ {
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+
+ ibottom_node=E->ien[j][kk].node[1];
+ zbottom=E->sx[j][3][ibottom_node];
+
+ if (zbottom<E->composition.z_interface) E->trace.oldel[j][kk]=1.0;
+ if (zbottom>=E->composition.z_interface) E->trace.oldel[j][kk]=0.0;
+
+ } /* end kk */
+ } /* end j */
+ }
+
+
+ /* Else read from file */
+
+
+ else if (E->trace.ic_method==2)
+ {
+
+ /* first look for backing file */
+
+ sprintf(output_file,"%s.comp_el.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
+ if ( (fp=fopen(output_file,"r"))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKMerror(Initialize Old Composition)-FILE NOT EXIST: %s\n", output_file);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ fgets(input_s,200,fp);
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ fgets(input_s,200,fp);
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ fgets(input_s,200,fp);
+ sscanf(input_s,"%lf",&E->trace.oldel[j][kk]);
+ }
+ }
+
+ fclose(fp);
+
+ } /* endif */
+
+
+
+ return;
+}
+
+
+
+/*********** COMPUTE ELEMENTAL COMPOSITION RATIO METHOD ***/
+/* */
+/* This function computes the composition per element. */
+/* This function computes the composition per element. */
+/* Integer array ieltrac stores tracers per element. */
+/* Double array celtrac stores the sum of tracer composition */
+
+void compute_elemental_composition_ratio_method(E)
+ struct All_variables *E;
+{
+
+ int kk;
+ int numtracers;
+ int nelem;
+ int j;
+ int iempty=0;
+
+ double comp;
+
+ static int been_here=0;
+
+ if (been_here==0)
+ {
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ if ((E->trace.ieltrac[j]=(int *)malloc((E->lmesh.nel+1)*sizeof(int)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 5u83a\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if ((E->trace.celtrac[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 58hy8\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if ((E->trace.comp_el[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 8989y\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+ been_here++;
+
+ }
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ /* first zero arrays */
+
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ E->trace.ieltrac[j][kk]=0;
+ E->trace.celtrac[j][kk]=0.0;
+ }
+
+ numtracers=E->trace.itrac[j][0];
+
+ /* Fill ieltrac and celtrac */
+
+
+ for (kk=1;kk<=numtracers;kk++)
+ {
+
+ nelem=E->trace.itrac[j][kk];
+ E->trace.ieltrac[j][nelem]++;
+
+ comp=E->trace.etrac[j][0][kk];
+
+ if (comp>1.0000001)
+ {
+ fprintf(E->trace.fpt,"ERROR(compute elemental)-not ready for comp>1 yet (%f)(tr. %d) \n",comp,kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ E->trace.celtrac[j][nelem]=E->trace.celtrac[j][nelem]+comp;
+
+ }
+
+ /* Check for empty entries and compute ratio. */
+ /* If no tracers are in an element, use previous composition */
+
+ iempty=0;
+
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+
+ if (E->trace.ieltrac[j][kk]==0)
+ {
+ iempty++;
+ E->trace.comp_el[j][kk]=E->trace.oldel[j][kk];
+ }
+ else if (E->trace.ieltrac[j][kk]>0)
+ {
+ E->trace.comp_el[j][kk]=E->trace.celtrac[j][kk]/(1.0*E->trace.ieltrac[j][kk]);
+ }
+
+ if (E->trace.comp_el[j][kk]>(1.000001) || E->trace.comp_el[j][kk]<(-0.000001))
+ {
+ fprintf(E->trace.fpt,"ERROR(compute elemental)-noway (3u5hd)\n");
+ fprintf(E->trace.fpt,"COMPEL: %f (%d)(%d)\n",E->trace.comp_el[j][kk],kk,E->trace.ieltrac[j][kk]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ if (iempty>0)
+ {
+
+ /*
+ fprintf(E->trace.fpt,"%d empty elements filled with old values (%f percent)\n",iempty, (100.0*iempty/E->lmesh.nel));
+ fflush(E->trace.fpt);
+ */
+
+ if ((1.0*iempty/E->lmesh.nel)>0.80)
+ {
+ fprintf(E->trace.fpt,"WARNING(compute_elemental...)-number of tracers is REALLY LOW\n");
+ fflush(E->trace.fpt);
+ if (E->trace.itracer_warnings==1) exit(10);
+ }
+ }
+
+ /* Fill oldel */
+
+
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ E->trace.oldel[j][kk]=E->trace.comp_el[j][kk];
+ }
+
+ } /* end j */
+
+ E->trace.istat_iempty=E->trace.istat_iempty+iempty;
+
+ return;
+}
+
+/********** MAP COMPOSITION TO NODES ****************/
+/* */
+
+
+void map_composition_to_nodes(E)
+ struct All_variables *E;
+{
+
+ int kk;
+ int nelem, nodenum;
+ int j;
+
+
+ static int been_here=0;
+
+ if (been_here==0)
+ {
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ if ((E->trace.comp_node[j]=(double *)malloc((E->lmesh.nno+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(map_compostion_to_nodes)-no memory 983rk\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+ been_here++;
+ }
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ /* first, initialize node array */
+
+ for (kk=1;kk<=E->lmesh.nno;kk++)
+ {
+ E->trace.comp_node[j][kk]=0.0;
+ }
+
+ /* Loop through all elements */
+
+ for (nelem=1;nelem<=E->lmesh.nel;nelem++)
+ {
+
+ /* for each element, loop through element nodes */
+
+ /* weight composition */
+
+ for (nodenum=1;nodenum<=8;nodenum++)
+ {
+
+ E->trace.comp_node[j][E->ien[j][nelem].node[nodenum]] +=
+ E->trace.comp_el[j][nelem]*
+ E->TWW[E->mesh.levmax][j][nelem].node[nodenum];
+
+ }
+
+ } /* end nelem */
+ } /* end j */
+
+ /* akm modified exchange node routine for doubles */
+
+ (E->exchange_node_d)(E,E->trace.comp_node,E->mesh.levmax);
+
+ /* Divide by nodal volume */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=E->lmesh.nno;kk++)
+ {
+ E->trace.comp_node[j][kk] *= E->MASS[E->mesh.levmax][j][kk];
+ }
+
+ /* testing */
+ /*
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ fprintf(E->trace.fpt,"%d %f\n",kk,E->trace.comp_el[j][kk]);
+ }
+
+ for (kk=1;kk<=E->lmesh.nno;kk++)
+ {
+ fprintf(E->trace.fpt,"%d %f %f\n",kk,E->sx[j][3][kk],E->trace.comp_node[j][kk]);
+ }
+ */
+
+
+ } /* end j */
+
+ return;
+}
+
+
+/*********** GET BULK COMPOSITION *******************************/
+
+void get_bulk_composition(E)
+ struct All_variables *E;
+
+{
+
+ char output_file[200];
+ char input_s[1000];
+
+ double return_bulk_value_d();
+ double volume;
+ double rdum1;
+ double rdum2;
+ double rdum3;
+
+ int ival=0;
+ int idum1;
+ int istep;
+
+
+ FILE *fp;
+
+ static int been_here=0;
+
+
+ /* ival=0 returns integral not average */
+
+ volume=return_bulk_value_d(E,E->composition.comp_node,ival);
+
+ E->composition.bulk_composition=volume;
+
+ /* Here we assume if restart = 1 or 0 tracers are reset */
+ /* if restart = 2 tracers may or may not be reset */
+ /* (read initial composition from file) */
+
+
+ if (been_here==0)
+ {
+ if (E->composition.ireset_initial_composition==1)
+ {
+ E->composition.initial_bulk_composition=volume;
+ }
+ else
+ {
+
+ if (E->trace.ic_method!=2)
+ {
+ fprintf(E->trace.fpt,"ERROR(bulk composition)-wrong reset,restart combo\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ sprintf(output_file,"%s.comp.%d.%d",E->control.old_P_file,
+ E->parallel.me,E->monitor.solution_cycles);
+
+ fp=fopen(output_file,"r");
+ fgets(input_s,200,fp);
+ sscanf(input_s,"%d %d %lf %lf %lf",
+ &istep,&idum1,&rdum1,&rdum2,&rdum3);
+
+ E->composition.initial_bulk_composition=rdum2;
+ fclose(fp);
+
+ if (istep!=E->monitor.solution_cycles)
+ {
+ fprintf(E->trace.fpt,"ERROR(get_bulk_composition) %d %d\n",
+ istep,E->monitor.solution_cycles);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
+
+ E->composition.error_fraction=((volume-E->composition.initial_bulk_composition)/
+ E->composition.initial_bulk_composition);
+
+ parallel_process_sync(E);
+
+ been_here++;
+ return;
+}
+
+
+
+/******************* READ COMP ***********************************************/
+/* */
+/* This function is similar to read_temp. It is used to read the composition */
+/* from file for post-proceesing. */
+
+void read_comp(E)
+ struct All_variables *E;
+{
+ int i,ii,m,mm,ll;
+ char output_file[255],input_s[1000];
+
+ double g;
+ FILE *fp;
+
+
+
+ ii = E->monitor.solution_cycles;
+ sprintf(output_file,"%s.comp.%d.%d",E->control.old_P_file,E->parallel.me,ii);
+
+ if ((fp=fopen(output_file,"r"))==NULL)
+ {
+ fprintf(stderr,"ERROR(read_temp) - %s not found\n",output_file);
+ fflush(stderr);
+ exit(10);
+ }
+
+ fgets(input_s,1000,fp);
+ sscanf(input_s,"%d %d %f",&ll,&mm,&E->monitor.elapsed_time);
+
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ {
+ E->trace.comp_node[m]=(double *)malloc((E->lmesh.nno+1)*sizeof(double));
+
+ fgets(input_s,1000,fp);
+ sscanf(input_s,"%d %d",&ll,&mm);
+ for(i=1;i<=E->lmesh.nno;i++)
+ {
+ if (fgets(input_s,1000,fp)==NULL)
+ {
+ fprintf(stderr,"ERROR(read_comp) -data for node %d not found\n",i);
+ fflush(stderr);
+ exit(10);
+ }
+ sscanf(input_s,"%lf",&g);
+ E->trace.comp_node[m][i] = g;
+
+ }
+ }
+
+ fclose (fp);
+ return;
+}
+
+
Modified: mc/3D/CitcomS/trunk/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Element_calculations.c 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/Element_calculations.c 2007-03-09 00:25:03 UTC (rev 6206)
@@ -68,16 +68,16 @@
double elt_f[24];
int m,a,e,i;
+ void get_buoyancy();
void get_elt_f();
void get_elt_tr();
void strip_bcs_from_residual();
- void thermal_buoyancy();
const int neq=E->lmesh.neq;
const int nel=E->lmesh.nel;
const int lev=E->mesh.levmax;
- thermal_buoyancy(E,E->buoyancy);
+ get_buoyancy(E,E->buoyancy);
for(m=1;m<=E->sphere.caps_per_proc;m++) {
@@ -114,16 +114,16 @@
double elt_f[24];
int m,a,e,i;
+ void get_buoyancy();
void get_elt_f();
void get_elt_tr_pseudo_surf();
void strip_bcs_from_residual();
- void thermal_buoyancy();
const int neq=E->lmesh.neq;
const int nel=E->lmesh.nel;
const int lev=E->mesh.levmax;
- thermal_buoyancy(E,E->buoyancy);
+ get_buoyancy(E,E->buoyancy);
for(m=1;m<=E->sphere.caps_per_proc;m++) {
Modified: mc/3D/CitcomS/trunk/lib/Full_obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_obsolete.c 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/Full_obsolete.c 2007-03-09 00:25:03 UTC (rev 6206)
@@ -796,12 +796,97 @@
/*************************************************************************/
-/* from */
+/* from Full_tracer_advection.c */
/*************************************************************************/
+/* */
+/* This function writes the radial distribution of tracers */
+/* (horizontally averaged) */
+void write_radial_horizontal_averages(E)
+ struct All_variables *E;
+{
+
+ char output_file[200];
+
+ int j;
+ int kk;
+ double halfpoint;
+ double *reltrac[13];
+
+ static int been_here=0;
+
+ void return_horiz_ave();
+ void return_elementwise_horiz_ave();
+
+ FILE *fp2;
+
+ if (been_here==0)
+ {
+ E->trace.Have_C=(double *)malloc((E->lmesh.noz+2)*sizeof(double));
+ E->trace.Havel_tracers=(double *)malloc((E->lmesh.elz+2)*sizeof(double));
+ }
+
+ /* Tracers */
+
+ /* first, change from int to double */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ reltrac[j]=(double *) malloc((E->lmesh.nel+1)*sizeof(double));
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ reltrac[j][kk]=(1.0*E->trace.ieltrac[j][kk]);
+ }
+ }
+
+ return_elementwise_horiz_ave(E,reltrac,E->trace.Havel_tracers);
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ free(reltrac[j]);
+ }
+
+ if (E->parallel.me<E->parallel.nprocz)
+ {
+ sprintf(output_file,"%s.ave_tracers.%d.%d",E->control.data_file,E->parallel.me,E->monitor.solution_cycles);
+ fp2=fopen(output_file,"w");
+ for(kk=1;kk<=E->lmesh.elz;kk++)
+ {
+ halfpoint=0.5*(E->sx[1][3][kk+1]+E->sx[1][3][kk]);
+ fprintf(fp2,"%.4e %.4e\n",halfpoint,E->trace.Havel_tracers[kk]);
+ }
+ fclose(fp2);
+ }
+
+ /* Composition */
+
+ if (E->composition.chemical_buoyancy==1)
+ {
+ return_horiz_ave(E,E->trace.comp_node,E->trace.Have_C);
+
+
+ if (E->parallel.me<E->parallel.nprocz)
+ {
+ sprintf(output_file,"%s.ave_c.%d.%d",E->control.data_file,E->parallel.me,E->monitor.solution_cycles);
+ fp2=fopen(output_file,"w");
+ for(kk=1;kk<=E->lmesh.noz;kk++)
+ {
+ fprintf(fp2,"%.4e %.4e\n",E->sx[1][3][kk],E->trace.Have_C[kk]);
+ }
+ fclose(fp2);
+
+ }
+ }
+
+ been_here++;
+
+ return;
+}
+
+
/*************************************************************************/
/* from */
/*************************************************************************/
Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-03-09 00:25:03 UTC (rev 6206)
@@ -31,10 +31,22 @@
#include "global_defs.h"
#include "parsing.h"
#include "parallel_related.h"
+#include "composition_related.h"
static void get_neighboring_caps(struct All_variables *E);
static void pdebug(struct All_variables *E, int i);
+static void fix_radius(struct All_variables *E,
+ double *radius, double *theta, double *phi,
+ double *x, double *y, double *z);
+static void fix_theta_phi(double *theta, double *phi);
+static void fix_phi(double *phi);
+static void predict_tracers(struct All_variables *E);
+static void correct_tracers(struct All_variables *E);
+static int isum_tracers(struct All_variables *E);
+
+
+
/******* FULL TRACER INPUT *********************/
void full_tracer_input(E)
@@ -43,95 +55,29 @@
{
int m = E->parallel.me;
- /* itracer_type=0 passive */
- /* itracer_type=1 active */
+ /* Initial condition, this option is ignored if E->control.restart is 1,
+ * ie. restarted from a previous run */
+ /* tracer_ic_method=0 (random generated array) */
+ /* tracer_ic_method=1 (all proc read the same file) */
+ /* tracer_ic_method=2 (each proc reads its restart file) */
+ if(E->control.restart)
+ E->trace.ic_method = 2;
+ else {
+ input_int("tracer_ic_method",&(E->trace.ic_method),"0,0,nomax",m);
- E->trace.itracer_type=1;
- input_int("tracer_type",&(E->trace.itracer_type),"1,0,nomax",m);
-
- /* Restart options */
-
- /* itracer_restart=-1 (read from file) */
- /* itracer_restart=0 (generate array) */
- /* itracer_restart=1 (read from scratch disks) */
-
- E->trace.itracer_restart=0;
- input_int("tracer_restart",&(E->trace.itracer_restart),"0,0,nomax",m);
-
- /* Active tracer inputs */
-
- if (E->trace.itracer_type==1)
- {
-
- E->trace.buoyancy_ratio=1.0;
- input_double("buoyancy_ratio",&(E->trace.buoyancy_ratio),"1.0",m);
-
- E->trace.ireset_initial_composition=1;
- if (E->trace.itracer_restart==1)
- {
- E->trace.ireset_initial_composition=0;
- input_int("reset_initial_composition",&(E->trace.ireset_initial_composition),"0",m);
- }
-
- /* compositional rheology */
-
- /* ibuoy_type=1 (ratio method) */
-
- E->trace.ibuoy_type=1;
- input_int("buoy_type",&(E->trace.ibuoy_type),"1,0,nomax",m);
- if (E->trace.ibuoy_type!=1)
- {
- fprintf(stderr,"Terror-Sorry, only ratio method allowed now\n");
- fflush(stderr);
- parallel_process_termination();
- }
-
- }
-
- /* icompositional_rheology=0 (off) */
- /* icompositional_rheology=1 (on) */
-
- E->trace.icompositional_rheology=0;
- input_int("compositional_rheology",
- &(E->trace.icompositional_rheology),"1,0,nomax",m);
-
- E->trace.compositional_rheology_prefactor=1.0;
-
- if (E->trace.icompositional_rheology>0)
- {
- input_double("compositional_prefactor",
- &(E->trace.compositional_rheology_prefactor),"1.0",m);
- }
-
- /* ibuoy_type=0 (absolute method) */
-
- /* icartesian_or_spherical=0 (cartesian coordinate input) */
- /* icartesian_or_spherical=1 (spherical coordinate input) */
-
- E->trace.icartesian_or_spherical_input=0;
- if (E->trace.itracer_restart==-1)
- {
- sprintf(E->trace.tracer_file,"tracer.dat");
+ if (E->trace.ic_method==0)
+ input_int("tracers_per_element",&(E->trace.itperel),"10,0,nomax",m);
+ else if (E->trace.ic_method==1)
input_string("tracer_file",E->trace.tracer_file,"tracer.dat",m);
+ else if (E->trace.ic_method==2) {
+ }
+ else {
+ fprintf(stderr,"Sorry, tracer_ic_method only 0, 1 and 2 available\n");
+ fflush(stderr);
+ parallel_process_termination();
+ }
+ }
- input_int("cartesian_or_spherical_input",&(E->trace.icartesian_or_spherical_input),"1,0,nomax",m);
-
- }
-
-
- if (E->trace.itracer_restart==0)
- {
- E->trace.itperel=10;
- input_int("tracers_per_element",&(E->trace.itperel),"10,0,nomax",m);
-
- E->trace.z_interface=0.5;
- input_double("z_interface",&(E->trace.z_interface),"0.5",m);
- }
-
- if (E->trace.itracer_restart==1)
- {
- }
-
/* Advection Scheme */
/* itracer_advection_scheme=1 (simple predictor corrector -uses only V(to)) */
@@ -155,11 +101,6 @@
}
- /* Output Options */
-
- E->trace.iwrite_tracers_every=1000000;
- input_int("write_tracers_every",&(E->trace.iwrite_tracers_every),"1000000,0,nomax",m);
-
/* Interpolation Scheme */
/* itracer_interpolation_scheme=1 (gnometric projection) */
/* itracer_interpolation_scheme=2 (simple average) */
@@ -190,6 +131,8 @@
input_int("analytical_tracer_test",&(E->trace.ianalytical_tracer_test),
"0,0,nomax",m);
+
+ composition_input(E);
return;
}
@@ -205,9 +148,7 @@
void viscosity_checks();
void make_tracer_array();
void initialize_old_composition();
- void fill_composition();
void find_tracers();
- void setup_dsincos();
void make_regular_grid();
void initialize_tracer_elements();
void define_uv_space();
@@ -222,7 +163,7 @@
E->trace.itracing=1;
/* some obscure initial parameters */
-
+ /* This parameter specifies how close a tracer can get to the boundary */
E->trace.box_cushion=0.00001;
/* AKMA turn this back on after debugging */
@@ -231,14 +172,15 @@
/* Determine number of tracer quantities */
/* advection_quantites - those needed for advection */
-
+ // TODO: generalize it
if (E->trace.itracer_advection_scheme==1) E->trace.number_of_advection_quantities=12;
if (E->trace.itracer_advection_scheme==2) E->trace.number_of_advection_quantities=12;
/* extra_quantities - used for composition, etc. */
/* (can be increased for additional science i.e. tracing chemistry */
- if (E->trace.itracer_type==1) E->trace.number_of_extra_quantities=1;
+ // TODO: how to move this part to Composition_related.c?
+ if (E->composition.ichemical_buoyancy==1) E->trace.number_of_extra_quantities=1;
else E->trace.number_of_extra_quantities=0;
E->trace.number_of_tracer_quantities=E->trace.number_of_advection_quantities +
@@ -256,11 +198,7 @@
sprintf(output_file,"%s.tracer_log.%d",E->control.data_file,m);
E->trace.fpt=fopen(output_file,"w");
- /* reset time counters */
- E->trace.stat_trace_time=0.0;
- E->trace.stat_stokes_time=0.0;
-
/* reset statistical counters */
E->trace.istat_isend=0;
@@ -295,17 +233,26 @@
if (E->sphere.caps_per_proc>1)
{
- fprintf(E->trace.fpt,"This code has not been tested or tried for multiple caps per processor!\n");
- fprintf(E->trace.fpt,"I would certainly do some testing first \n");
+ fprintf(E->trace.fpt,"This code does not work for multiple caps per processor!\n");
fflush(E->trace.fpt);
parallel_process_termination();
}
get_neighboring_caps(E);
- if (E->trace.itracer_restart==0) make_tracer_array(E);
- else if (E->trace.itracer_restart==-1) read_tracer_file(E);
- else if (E->trace.itracer_restart==1) restart_tracers(E);
+ if (E->trace.ic_method==0) {
+ make_tracer_array(E);
+
+ if (E->composition.ichemical_buoyancy==1)
+ init_tracer_composition(E);
+ }
+ else if (E->trace.ic_method==1) {
+ read_tracer_file(E);
+
+ if (E->composition.ichemical_buoyancy==1)
+ init_tracer_composition(E);
+ }
+ else if (E->trace.ic_method==2) restart_tracers(E);
else
{
fprintf(E->trace.fpt,"Not ready for other inputs yet\n");
@@ -318,8 +265,6 @@
parallel_process_sync(E);
- setup_dsincos(E);
-
if (E->trace.itracer_interpolation_scheme==1)
{
define_uv_space(E);
@@ -342,20 +287,20 @@
find_tracers(E);
- check_sum(E);
+ /* total number of tracers */
- if (E->trace.itracer_type==1 && E->trace.ibuoy_type==1)
- {
- initialize_old_composition(E);
- fill_composition(E);
- }
+ E->trace.ilast_tracer_count = isum_tracers(E);
+ fprintf(E->trace.fpt, "Sum of Tracers: %d\n", E->trace.ilast_tracer_count);
if (E->trace.ianalytical_tracer_test==1)
{
+ //TODO: walk into this code...
analytical_test(E);
parallel_process_termination();
}
+ composition_setup(E);
+
tracer_post_processing(E);
return;
@@ -373,30 +318,19 @@
struct All_variables *E;
{
- static double start_time;
- static double end_time;
-
- void fill_composition();
- void thermal_buoyancy();
void check_sum();
void tracer_post_processing();
void write_tracers();
- void predict_tracers();
- void correct_tracers();
- double CPU_time0();
- start_time=CPU_time0();
-
- E->trace.stat_trace_time=0.0;
-
fprintf(E->trace.fpt,"STEP %d\n",E->monitor.solution_cycles);
fflush(E->trace.fpt);
/* presave last timesteps tracers as preback */
-
+ //TODO: use system("mv oldfile oldfile.bak\n") instead
if ( (E->monitor.solution_cycles % E->control.record_every)==0) {
- write_tracers(E,3);
+ //TODO: migrate to output_tracer
+ //write_tracers(E,3);
}
/* advect tracers */
@@ -404,106 +338,20 @@
predict_tracers(E);
correct_tracers(E);
- if (E->trace.itracer_type==1) {
+ check_sum(E);
+
+ //TODO: move
+ if (E->composition.ichemical_buoyancy==1) {
fill_composition(E);
}
- end_time=CPU_time0();
- E->trace.stat_trace_time=E->trace.stat_trace_time+(end_time-start_time);
-
- check_sum(E);
tracer_post_processing(E);
return;
}
-void tracing(E,itimes_here_this_step)
- struct All_variables *E;
- int itimes_here_this_step;
-{
- static double start_time;
- static double end_time;
-
- void fill_composition();
- void thermal_buoyancy();
- void check_sum();
- void tracer_post_processing();
- void write_tracers();
- void predict_tracers();
- void correct_tracers();
- double CPU_time0();
-
-
- parallel_process_sync(E);
- start_time=CPU_time0();
-
- if (itimes_here_this_step==1)
- {
-
- E->trace.stat_trace_time=0.0;
-
- fprintf(E->trace.fpt,"STEP %d\n",E->monitor.solution_cycles);
- fflush(E->trace.fpt);
-
- /* presave last timesteps tracers as preback */
-
- if ( (E->monitor.solution_cycles % E->control.record_every)==0)
- {
- write_tracers(E,3);
- }
- }
- else if (itimes_here_this_step==2)
- {
- E->trace.stat_stokes_time=start_time-end_time;
- }
-
- /* advect tracers */
-
- if (itimes_here_this_step==1)
- {
- if (E->trace.itracer_advection_scheme==1)
- {
- predict_tracers(E);
- correct_tracers(E);
- }
- else if (E->trace.itracer_advection_scheme==2)
- {
- predict_tracers(E);
- }
- }
- else if (itimes_here_this_step==2)
- {
- if (E->trace.itracer_advection_scheme==1)
- {
- }
- else if (E->trace.itracer_advection_scheme==2)
- {
- correct_tracers(E);
- }
- }
-
-
-
- if (E->trace.itracer_type==1)
- {
- fill_composition(E);
- }
-
- end_time=CPU_time0();
- E->trace.stat_trace_time=E->trace.stat_trace_time+(end_time-start_time);
-
-
- if (itimes_here_this_step==2)
- {
- check_sum(E);
- tracer_post_processing(E);
- }
-
- return;
-}
-
/********* TRACER POST PROCESSING ****************************************/
void tracer_post_processing(E)
@@ -517,25 +365,24 @@
double trace_fraction,total_time;
void get_bulk_composition();
- void write_compositional_field();
void write_tracers();
- void write_radial_horizontal_averages();
static int been_here=0;
- //if (E->trace.itracer_type==1) get_bulk_composition(E);
+ //TODO: fix this function
+ //if (E->composition.ichemical_buoyancy==1) get_bulk_composition(E);
if ( ((E->monitor.solution_cycles % E->control.record_every)==0) || (been_here==0))
{
- //TODO: calling this function will crash CitcomS
- //write_radial_horizontal_averages(E);
- write_tracers(E,1);
- //if (E->trace.itracer_type==1) write_compositional_field(E);
+ //TODO: migrate to output_tracer
+ //write_tracers(E,1);
}
- if ( ((E->monitor.solution_cycles % E->trace.iwrite_tracers_every)==0) &&
- ((E->monitor.solution_cycles!=E->control.restart)&&(E->monitor.solution_cycles!=0) ))
+
+ //TODO: move to Output.c
+ if ( ((E->monitor.solution_cycles!=E->control.restart)&&(E->monitor.solution_cycles!=0) ))
{
- write_tracers(E,2);
+ //TODO: migrate to output_tracer
+ //write_tracers(E,2);
}
@@ -544,7 +391,7 @@
{
fprintf(E->trace.fpt,"Number of tracers sent to other processors: %d\n",E->trace.istat_isend);
fprintf(E->trace.fpt,"Number of times element columns are checked: %d \n",E->trace.istat_elements_checked);
- if (E->trace.itracer_type==1)
+ if (E->composition.ichemical_buoyancy==1)
{
fprintf(E->trace.fpt,"Empty elements filled with old compositional values: %d (%f percent)\n",
E->trace.istat_iempty,(100.0*E->trace.istat_iempty/E->lmesh.nel));
@@ -563,12 +410,12 @@
E->trace.istat1=0;
/* compositional and error fraction data files */
-
+ //TODO: move
if (E->parallel.me==0)
{
if (been_here==0)
{
- if (E->trace.itracer_type==1)
+ if (E->composition.ichemical_buoyancy==1)
{
sprintf(output_file,"%s.error_fraction.data",E->control.data_file);
E->trace.fp_error_fraction=fopen(output_file,"w");
@@ -578,8 +425,9 @@
}
}
- if (E->trace.itracer_type==1)
+ if (E->composition.ichemical_buoyancy==1)
{
+ //TODO: to be init'd
fprintf(E->trace.fp_error_fraction,"%e %e\n",E->monitor.elapsed_time,E->trace.error_fraction);
fprintf(E->trace.fp_composition,"%e %e\n",E->monitor.elapsed_time,E->trace.bulk_composition);
@@ -590,19 +438,7 @@
}
- convection_time=E->trace.stat_stokes_time;
- tracer_time=E->trace.stat_trace_time;
- total_time=convection_time+tracer_time;
- trace_fraction=tracer_time/total_time;
-
- if (been_here!=0)
- {
- fprintf(E->trace.fpt,"CPU time: Convection: %f Tracing: %f (%f percent) Total: %f\n",
- convection_time,tracer_time,trace_fraction*100.0,total_time);
-
- }
-
fflush(E->trace.fpt);
if (been_here==0) been_here++;
@@ -610,92 +446,6 @@
return;
}
-/*********** WRITE RADIAL HORIZONTAL AVERAGES ***************************/
-/* */
-/* This function writes the radial distribution of tracers */
-/* (horizontally averaged) */
-
-void write_radial_horizontal_averages(E)
- struct All_variables *E;
-{
-
- char output_file[200];
-
- int j;
- int kk;
- double halfpoint;
- double *reltrac[13];
-
- static int been_here=0;
-
- void return_horiz_ave();
- void return_elementwise_horiz_ave();
-
- FILE *fp2;
-
- if (been_here==0)
- {
- E->trace.Have_C=(double *)malloc((E->lmesh.noz+2)*sizeof(double));
- E->trace.Havel_tracers=(double *)malloc((E->lmesh.elz+2)*sizeof(double));
- }
-
- /* Tracers */
-
- /* first, change from int to double */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- reltrac[j]=(double *) malloc((E->lmesh.nel+1)*sizeof(double));
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
- reltrac[j][kk]=(1.0*E->trace.ieltrac[j][kk]);
- }
- }
-
- return_elementwise_horiz_ave(E,reltrac,E->trace.Havel_tracers);
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- free(reltrac[j]);
- }
-
- if (E->parallel.me<E->parallel.nprocz)
- {
- sprintf(output_file,"%s.ave_tracers.%d.%d",E->control.data_file,E->parallel.me,E->monitor.solution_cycles);
- fp2=fopen(output_file,"w");
- for(kk=1;kk<=E->lmesh.elz;kk++)
- {
- halfpoint=0.5*(E->sx[1][3][kk+1]+E->sx[1][3][kk]);
- fprintf(fp2,"%.4e %.4e\n",halfpoint,E->trace.Havel_tracers[kk]);
- }
- fclose(fp2);
- }
-
- /* Composition */
-
- if (E->trace.itracer_type==1)
- {
- return_horiz_ave(E,E->trace.comp_node,E->trace.Have_C);
-
-
- if (E->parallel.me<E->parallel.nprocz)
- {
- sprintf(output_file,"%s.ave_c.%d.%d",E->control.data_file,E->parallel.me,E->monitor.solution_cycles);
- fp2=fopen(output_file,"w");
- for(kk=1;kk<=E->lmesh.noz;kk++)
- {
- fprintf(fp2,"%.4e %.4e\n",E->sx[1][3][kk],E->trace.Have_C[kk]);
- }
- fclose(fp2);
-
- }
- }
-
- been_here++;
-
- return;
-}
-
/*********** WRITE TRACERS **********************************************/
/* */
/* if iflag=1, rewrite over last save array */
@@ -738,7 +488,7 @@
for(j=1;j<=E->sphere.caps_per_proc;j++)
{
-
+ //TODO: move
/* fprintf(fpcomp,"%6d %6d %.5e %.5e %.5e %d\n",E->trace.itrac[j][0], */
/* E->monitor.solution_cycles,E->monitor.elapsed_time, */
/* E->trace.bulk_composition, */
@@ -751,13 +501,13 @@
phi=E->trace.rtrac[j][1][kk];
rad=E->trace.rtrac[j][2][kk];
- if (E->trace.itracer_type==1)
+ if (E->composition.ichemical_buoyancy==1)
{
comp=E->trace.etrac[j][0][kk];
fprintf(fpcomp,"%.12e %.12e %.12e %.12e \n",theta,phi,rad,comp);
fflush(fpcomp);
}
- else if (E->trace.itracer_type==0)
+ else if (E->composition.ichemical_buoyancy==0)
{
fprintf(fpcomp,"%.12e %.12e %.12e \n",theta,phi,rad);
fflush(fpcomp);
@@ -780,169 +530,6 @@
}
-
-
-
-/********** WRITE COMPOSITIONAL FIELD ***********************************/
-
-void write_compositional_field(E)
- struct All_variables *E;
-{
-
- char output_file1[1000];
- char output_file2[1000];
-
- int j,i;
-
- FILE *fp1;
- FILE *fp2;
-
-
- /* comp */
-
- sprintf(output_file1,"%s.comp.%d.%d",E->control.data_file,E->parallel.me,E->monitor.solution_cycles);
-
- fp1=fopen(output_file1,"w");
-
- fprintf(fp1,"%d %d %.5e %.5e %.5e\n",E->monitor.solution_cycles,E->lmesh.nno,E->monitor.elapsed_time,
- E->trace.initial_bulk_composition,E->trace.bulk_composition);
-
-
- for(j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- fprintf(fp1,"%3d %7d\n",j,E->lmesh.nno);
-
- for(i=1;i<=E->lmesh.nno;i++)
- {
- fprintf(fp1,"%.6e\n",E->trace.comp_node[j][i]);
- }
- }
-
- fflush(fp1);
-
- /* OLDCOMP */
-
- sprintf(output_file2,"%s.OLDCOMP.%d",E->control.data_file,E->parallel.me);
- fp2=fopen(output_file2,"w");
- fprintf(fp2,"%d %d %.5e \n",E->lmesh.nel,E->monitor.solution_cycles,E->monitor.elapsed_time);
- for(j=1;j<=E->sphere.caps_per_proc;j++)
- {
- fprintf(fp2,"%d %d\n",j,E->lmesh.nel);
- for(i=1;i<=E->lmesh.nel;i++)
- {
- fprintf(fp2,"%.6e\n",E->trace.oldel[j][i]);
- }
- }
- fflush(fp2);
-
- if ( ((E->monitor.solution_cycles % E->trace.iwrite_tracers_every)==0) &&
- ((E->monitor.solution_cycles!=E->control.restart)&&(E->monitor.solution_cycles!=0) ))
- {
- sprintf(output_file2,"%s.OLDCOMP.%d.%d",E->control.data_file,E->parallel.me,E->monitor.solution_cycles);
- fp2=fopen(output_file2,"w");
- fprintf(fp2,"%d %d %.5e \n",E->lmesh.nel,E->monitor.solution_cycles,E->monitor.elapsed_time);
- for(j=1;j<=E->sphere.caps_per_proc;j++)
- {
- fprintf(fp2,"%d %d\n",j,E->lmesh.nel);
- for(i=1;i<=E->lmesh.nel;i++)
- {
- fprintf(fp2,"%.6e\n",E->trace.oldel[j][i]);
- }
- }
-
- }
- fflush(fp2);
-
-
-
- return;
-}
-
-/*********** GET BULK COMPOSITION *******************************/
-
-void get_bulk_composition(E)
- struct All_variables *E;
-
-{
-
- char output_file[200];
- char input_s[1000];
-
- double return_bulk_value_d();
- double volume;
- double rdum1;
- double rdum2;
- double rdum3;
-
- int ival=0;
- int idum1;
- int istep;
-
-
- FILE *fp;
-
- static int been_here=0;
-
-
- /* ival=0 returns integral not average */
-
- volume=return_bulk_value_d(E,E->trace.comp_node,ival);
-
- E->trace.bulk_composition=volume;
-
- /* Here we assume if restart = -1 or 0 tracers are reset */
- /* if restart = 1 tracers may or may not be reset */
- /* (read initial composition from file) */
-
-
- if (been_here==0)
- {
- if (E->trace.ireset_initial_composition==1)
- {
- E->trace.initial_bulk_composition=volume;
- }
- else
- {
-
- if (E->trace.itracer_restart!=1)
- {
- fprintf(E->trace.fpt,"ERROR(bulk composition)-wrong reset,restart combo\n");
- fflush(E->trace.fpt);
- exit(10);
- }
-
- sprintf(output_file,"%s.comp.%d.%d",E->control.old_P_file,
- E->parallel.me,E->monitor.solution_cycles);
-
- fp=fopen(output_file,"r");
- fgets(input_s,200,fp);
- sscanf(input_s,"%d %d %lf %lf %lf",
- &istep,&idum1,&rdum1,&rdum2,&rdum3);
-
- E->trace.initial_bulk_composition=rdum2;
- fclose(fp);
-
- if (istep!=E->monitor.solution_cycles)
- {
- fprintf(E->trace.fpt,"ERROR(get_bulk_composition) %d %d\n",
- istep,E->monitor.solution_cycles);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- }
-
- E->trace.error_fraction=((volume-E->trace.initial_bulk_composition)/
- E->trace.initial_bulk_composition);
-
- parallel_process_sync(E);
-
- been_here++;
- return;
-}
-
-
/*********** PREDICT TRACERS **********************************************/
/* */
/* This function predicts tracers performing an euler step */
@@ -956,8 +543,7 @@
/* */
-void predict_tracers(E)
- struct All_variables *E;
+static void predict_tracers(struct All_variables *E)
{
int numtracers;
@@ -1069,8 +655,7 @@
/* */
-void correct_tracers(E)
- struct All_variables *E;
+static void correct_tracers(struct All_variables *E)
{
int j;
@@ -1255,7 +840,7 @@
int maxlevel=E->mesh.levmax;
- static double eps=-1e-4;
+ const double eps=-1e-4;
void get_radial_shape();
void sphere_to_cart();
@@ -1469,7 +1054,7 @@
int node1,node5;
double rad1,rad5,f1,f2,delrad;
- static double eps=1e-6;
+ const double eps=1e-6;
double top_bound=1.0+eps;
double bottom_bound=0.0-eps;
@@ -2819,15 +2404,21 @@
int icheck_element_column();
int icheck_column_neighbors();
void sphere_to_cart();
- void fix_theta();
- void fix_phi();
- static double two_pi=2.0*M_PI;
+ const double two_pi=2.0*M_PI;
elz=E->lmesh.elz;
nelsurf=E->lmesh.elx*E->lmesh.ely;
+ //TODO: find the bounding box of the mesh, if the box is too close to
+ // to core, set a flag (rotated_reggrid) to true and rotate the
+ // bounding box to the equator. Generate the regular grid with the new
+ // bounding box. The rotation should be a simple one, e.g.
+ // (theta, phi) -> (??)
+ // Whenever the regular grid is used, check the flat (rotated_reggrid),
+ // if true, rotate the checkpoint as well.
+
/* note, mesh is rotated along theta 22.5 degrees divided by elx. */
/* We at least want that much expansion here! Otherwise, theta min */
/* will not be valid near poles. We do a little more (x2) to be safe */
@@ -3429,12 +3020,12 @@
int icheck_element_column();
/*
- static int number_of_neighbors=24;
+ const int number_of_neighbors=24;
*/
/* maybe faster to only check inner ring */
- static int number_of_neighbors=8;
+ const int number_of_neighbors=8;
elx=E->lmesh.elx;
ely=E->lmesh.ely;
@@ -3522,444 +3113,30 @@
}
-/************ INITIALIZE OLD COMPOSITION ************************/
-void initialize_old_composition(E)
- struct All_variables *E;
-{
- char output_file[200];
- char input_s[1000];
-
- int ibottom_node;
- int kk;
- int istep;
- int idum;
- int j;
-
- double zbottom;
- double time;
-
- FILE *fp;
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- if ((E->trace.oldel[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(fill old composition)-no memory 324c\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
-
-
- if ((E->trace.itracer_restart==0)||(E->trace.itracer_restart==-1))
- {
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
-
- ibottom_node=E->ien[j][kk].node[1];
- zbottom=E->sx[j][3][ibottom_node];
-
- if (zbottom<E->trace.z_interface) E->trace.oldel[j][kk]=1.0;
- if (zbottom>=E->trace.z_interface) E->trace.oldel[j][kk]=0.0;
-
- } /* end kk */
- } /* end j */
- }
-
-
- /* Else read from file */
-
-
- else if (E->trace.itracer_restart==1)
- {
-
- /* first look for backing file */
-
-
- sprintf(output_file,"%s.OLDCOMP.%d.%d",E->control.old_P_file,E->parallel.me,E->control.restart);
- if ( (fp=fopen(output_file,"r"))==NULL)
- {
- sprintf(output_file,"%s.OLDCOMP.%d",E->control.old_P_file,E->parallel.me);
-
- if ( (fp=fopen(output_file,"r"))==NULL)
- {
- fprintf(E->trace.fpt,"AKMerror(Initialize Old Composition)-OLDCOMP NOT EXIST\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
-
- fgets(input_s,200,fp);
- sscanf(input_s,"%d %d %lf",&idum,&istep,&time);
-
- if (istep!=E->control.restart)
- {
- fprintf(E->trace.fpt,"Error(initialize old comp) %d %d\n",
- istep,E->monitor.solution_cycles);
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
- for(j=1;j<=E->sphere.caps_per_proc;j++)
- {
- fgets(input_s,200,fp);
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
- fgets(input_s,200,fp);
- sscanf(input_s,"%lf",&E->trace.oldel[j][kk]);
- }
- }
-
- fclose(fp);
-
- } /* endif */
-
-
-
- return;
-}
-
-/************ FILL COMPOSITION ************************/
-void fill_composition(E)
- struct All_variables *E;
-{
-
- void compute_elemental_composition_ratio_method();
- void map_composition_to_nodes();
-
- /* Currently, only the ratio method works here. */
- /* Will have to come back here to include the absolute method. */
-
- /* ratio method */
-
- if (E->trace.ibuoy_type==1)
- {
- compute_elemental_composition_ratio_method(E);
- }
-
- /* absolute method */
-
- if (E->trace.ibuoy_type!=1)
- {
- fprintf(E->trace.fpt,"Error(compute...)-only ratio method now\n");
- fflush(E->trace.fpt);
- exit(10);
- }
-
- /* Map elemental composition to nodal points */
-
- map_composition_to_nodes(E);
-
- return;
-}
-
-/*********** COMPUTE ELEMENTAL COMPOSITION RATIO METHOD ***/
-/* */
-/* This function computes the composition per element. */
-/* This function computes the composition per element. */
-/* Integer array ieltrac stores tracers per element. */
-/* Double array celtrac stores the sum of tracer composition */
-
-void compute_elemental_composition_ratio_method(E)
- struct All_variables *E;
-{
-
- int kk;
- int numtracers;
- int nelem;
- int j;
- int iempty=0;
-
- double comp;
-
- static int been_here=0;
-
- if (been_here==0)
- {
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- if ((E->trace.ieltrac[j]=(int *)malloc((E->lmesh.nel+1)*sizeof(int)))==NULL)
- {
- fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 5u83a\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- if ((E->trace.celtrac[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 58hy8\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- if ((E->trace.comp_el[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 8989y\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
-
- been_here++;
-
- }
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- /* first zero arrays */
-
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
- E->trace.ieltrac[j][kk]=0;
- E->trace.celtrac[j][kk]=0.0;
- }
-
- numtracers=E->trace.itrac[j][0];
-
- /* Fill ieltrac and celtrac */
-
-
- for (kk=1;kk<=numtracers;kk++)
- {
-
- nelem=E->trace.itrac[j][kk];
- E->trace.ieltrac[j][nelem]++;
-
- comp=E->trace.etrac[j][0][kk];
-
- if (comp>1.0000001)
- {
- fprintf(E->trace.fpt,"ERROR(compute elemental)-not ready for comp>1 yet (%f)(tr. %d) \n",comp,kk);
- fflush(E->trace.fpt);
- exit(10);
- }
-
- E->trace.celtrac[j][nelem]=E->trace.celtrac[j][nelem]+comp;
-
- }
-
- /* Check for empty entries and compute ratio. */
- /* If no tracers are in an element, use previous composition */
-
- iempty=0;
-
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
-
- if (E->trace.ieltrac[j][kk]==0)
- {
- iempty++;
- E->trace.comp_el[j][kk]=E->trace.oldel[j][kk];
- }
- else if (E->trace.ieltrac[j][kk]>0)
- {
- E->trace.comp_el[j][kk]=E->trace.celtrac[j][kk]/(1.0*E->trace.ieltrac[j][kk]);
- }
-
- if (E->trace.comp_el[j][kk]>(1.000001) || E->trace.comp_el[j][kk]<(-0.000001))
- {
- fprintf(E->trace.fpt,"ERROR(compute elemental)-noway (3u5hd)\n");
- fprintf(E->trace.fpt,"COMPEL: %f (%d)(%d)\n",E->trace.comp_el[j][kk],kk,E->trace.ieltrac[j][kk]);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- if (iempty>0)
- {
-
- /*
- fprintf(E->trace.fpt,"%d empty elements filled with old values (%f percent)\n",iempty, (100.0*iempty/E->lmesh.nel));
- fflush(E->trace.fpt);
- */
-
- if ((1.0*iempty/E->lmesh.nel)>0.80)
- {
- fprintf(E->trace.fpt,"WARNING(compute_elemental...)-number of tracers is REALLY LOW\n");
- fflush(E->trace.fpt);
- if (E->trace.itracer_warnings==1) exit(10);
- }
- }
-
- /* Fill oldel */
-
-
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
- E->trace.oldel[j][kk]=E->trace.comp_el[j][kk];
- }
-
- } /* end j */
-
- E->trace.istat_iempty=E->trace.istat_iempty+iempty;
-
- return;
-}
-
-/********** MAP COMPOSITION TO NODES ****************/
-/* */
-
-
-void map_composition_to_nodes(E)
- struct All_variables *E;
-{
-
- int kk;
- int nelem, nodenum;
- int j;
-
-
- static int been_here=0;
-
- if (been_here==0)
- {
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- if ((E->trace.comp_node[j]=(double *)malloc((E->lmesh.nno+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"AKM(map_compostion_to_nodes)-no memory 983rk\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
-
- been_here++;
- }
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- /* first, initialize node array */
-
- for (kk=1;kk<=E->lmesh.nno;kk++)
- {
- E->trace.comp_node[j][kk]=0.0;
- }
-
- /* Loop through all elements */
-
- for (nelem=1;nelem<=E->lmesh.nel;nelem++)
- {
-
- /* for each element, loop through element nodes */
-
- /* weight composition */
-
- for (nodenum=1;nodenum<=8;nodenum++)
- {
-
- E->trace.comp_node[j][E->ien[j][nelem].node[nodenum]] +=
- E->trace.comp_el[j][nelem]*
- E->TWW[E->mesh.levmax][j][nelem].node[nodenum];
-
- }
-
- } /* end nelem */
- } /* end j */
-
- /* akm modified exchange node routine for doubles */
-
- (E->exchange_node_d)(E,E->trace.comp_node,E->mesh.levmax);
-
- /* Divide by nodal volume */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->lmesh.nno;kk++)
- {
- E->trace.comp_node[j][kk] *= E->MASS[E->mesh.levmax][j][kk];
- }
-
- /* testing */
- /*
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
- fprintf(E->trace.fpt,"%d %f\n",kk,E->trace.comp_el[j][kk]);
- }
-
- for (kk=1;kk<=E->lmesh.nno;kk++)
- {
- fprintf(E->trace.fpt,"%d %f %f\n",kk,E->sx[j][3][kk],E->trace.comp_node[j][kk]);
- }
- */
-
-
- } /* end j */
-
- return;
-}
-
/**** WRITE TRACE INSTRUCTIONS ***************/
void write_trace_instructions(E)
struct All_variables *E;
{
-
-
-
-
fprintf(E->trace.fpt,"\nTracing Activated! (proc: %d)\n",E->parallel.me);
fprintf(E->trace.fpt," Allen K. McNamara 12-2003\n\n");
- if (E->trace.itracer_type==0)
+ if (E->trace.ic_method==0)
{
- fprintf(E->trace.fpt,"Passive Tracers\n");
- }
-
- if (E->trace.itracer_type==1)
- {
- fprintf(E->trace.fpt,"Active Tracers\n");
- if (E->trace.ibuoy_type==1) fprintf(E->trace.fpt,"Ratio Method\n");
- if (E->trace.ibuoy_type==0) fprintf(E->trace.fpt,"Absolute Method\n");
- fprintf(E->trace.fpt,"Buoyancy Ratio: %f\n", E->trace.buoyancy_ratio);
-
- if (E->trace.ireset_initial_composition==0)
- {
- fprintf(E->trace.fpt,"Using old initial composition from tracer files\n");
- }
- else
- {
- fprintf(E->trace.fpt,"Resetting initial composition\n");
- }
-
- fflush(E->trace.fpt);
- fflush(stderr);
-
- }
-
- if (E->trace.itracer_restart==0)
- {
fprintf(E->trace.fpt,"Generating New Tracer Array\n");
fprintf(E->trace.fpt,"Tracers per element: %d\n",E->trace.itperel);
- if (E->trace.itracer_type==1)
+ /* TODO: move
+ if (E->composition.ichemical_buoyancy==1)
{
- fprintf(E->trace.fpt,"Interface Height: %f\n",E->trace.z_interface);
+ fprintf(E->trace.fpt,"Interface Height: %f\n",E->composition.z_interface);
}
+ */
}
- if (E->trace.itracer_restart==-1)
+ if (E->trace.ic_method==1)
{
fprintf(E->trace.fpt,"Reading tracer file %s\n",E->trace.tracer_file);
- if (E->trace.icartesian_or_spherical_input==0)
- {
- fprintf(E->trace.fpt,"Coordinates are read as Cartesian\n");
- fprintf(E->trace.fpt," x,y,z format\n");
- }
- else if (E->trace.icartesian_or_spherical_input==1)
- {
- fprintf(E->trace.fpt,"Coordinates are read as Spherical\n");
- fprintf(E->trace.fpt," theta,phi,rad format\n");
- }
- else
- {
- fprintf(E->trace.fpt,"ERROR-Cartesian or Spherical input ?\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
}
- if (E->trace.itracer_restart==1)
+ if (E->trace.ic_method==2)
{
fprintf(E->trace.fpt,"Restarting Tracers\n");
}
@@ -4004,21 +3181,7 @@
parallel_process_termination();
}
- fprintf(E->trace.fpt,"Writing Tracers to File every %d steps\n",
- E->trace.iwrite_tracers_every);
-
- if (E->trace.icompositional_rheology==0)
- {
- fprintf(E->trace.fpt,"Compositional Rheology - OFF\n");
- }
- else if (E->trace.icompositional_rheology>0)
- {
- fprintf(E->trace.fpt,"Compositional Rheology - ON\n");
- fprintf(E->trace.fpt,"Compositional Prefactor: %f\n",
- E->trace.compositional_rheology_prefactor);
- }
-
/* regular grid stuff */
fprintf(E->trace.fpt,"Regular Grid-> deltheta: %f delphi: %f\n",
@@ -4056,6 +3219,7 @@
fflush(stderr);
}
+ write_composition_instructions(E);
return;
}
@@ -4073,7 +3237,7 @@
char input_s[1000];
int j,kk;
- int idum1,istep;
+ int idum1,ncolumns;
int numtracers;
double rdum1;
@@ -4088,22 +3252,14 @@
FILE *fp1;
- /* first try old storage file */
+ sprintf(output_file,"%s.tracer.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
- sprintf(output_file,"%s.tracers.%d.%d",E->control.old_P_file,E->parallel.me,E->control.restart);
-
- /* if not, try last backing file */
-
if ( (fp1=fopen(output_file,"r"))==NULL)
- {
- sprintf(output_file,"%s.tracers.%d",E->control.old_P_file,E->parallel.me);
- if ( (fp1=fopen(output_file,"r"))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(restart tracers)-file not found %s\n",output_file);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-file not found %s\n",output_file);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
fprintf(stderr,"Restarting Tracers from %s\n",output_file);
fflush(stderr);
@@ -4112,20 +3268,21 @@
for(j=1;j<=E->sphere.caps_per_proc;j++)
{
fgets(input_s,200,fp1);
- sscanf(input_s,"%d %d %lf %lf %lf %d",
- &numtracers,&istep,&rdum1,&E->trace.bulk_composition,&E->trace.initial_bulk_composition,&idum1);
+ sscanf(input_s,"%d %d %d %lf",
+ &idum1, &numtracers, &ncolumns, &rdum1);
+ /* some error control */
+ if (E->composition.ichemical_buoyancy==0 && ncolumns!=3) {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-wrong # of columns\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if (E->composition.ichemical_buoyancy==1 && ncolumns!=4) {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-wrong # of columns\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
-
- /* some error control */
-
- if (istep!=E->control.restart)
- {
- fprintf(E->trace.fpt,"ERROR(restart tracers)- step? %d %d\n",istep,E->control.restart);
- fflush(E->trace.fpt);
- exit(10);
- }
-
/* allocate memory for tracer arrays */
initialize_tracer_arrays(E,j,numtracers);
@@ -4135,11 +3292,12 @@
{
comp=0.0;
fgets(input_s,200,fp1);
- if (E->trace.itracer_type==0)
+ if (E->composition.ichemical_buoyancy==0)
{
sscanf(input_s,"%lf %lf %lf\n",&theta,&phi,&rad);
}
- else if (E->trace.itracer_type==1)
+ //TODO: move
+ else if (E->composition.ichemical_buoyancy==1)
{
sscanf(input_s,"%lf %lf %lf %lf\n",&theta,&phi,&rad,&comp);
}
@@ -4162,7 +3320,8 @@
E->trace.rtrac[j][3][kk]=x;
E->trace.rtrac[j][4][kk]=y;
E->trace.rtrac[j][5][kk]=z;
- if (E->trace.itracer_type==1) E->trace.etrac[j][0][kk]=comp;
+ //TODO: move
+ if (E->composition.ichemical_buoyancy==1) E->trace.etrac[j][0][kk]=comp;
}
@@ -4186,7 +3345,6 @@
int kk;
int tracers_cap;
int j;
- int ithiscap;
int ival;
int number_of_tries=0;
int max_tries;
@@ -4210,14 +3368,9 @@
for (j=1;j<=E->sphere.caps_per_proc;j++)
{
- ithiscap=E->sphere.capid[j];
-
-
processor_fraction=( ( pow(E->sx[j][3][E->lmesh.noz],3.0)-pow(E->sx[j][3][1],3.0))/
(pow(E->sphere.ro,3.0)-pow(E->sphere.ri,3.0)));
-
- tracers_cap=((E->lmesh.nel*E->parallel.nprocz)*E->trace.itperel
- *processor_fraction);
+ tracers_cap=E->mesh.nel*E->trace.itperel*processor_fraction;
/*
fprintf(stderr,"AA: proc frac: %f (%d) %d %d %f %f\n",processor_fraction,tracers_cap,E->lmesh.nel,E->parallel.nprocz, E->sx[j][3][E->lmesh.noz],E->sx[j][3][1]);
*/
@@ -4286,13 +3439,6 @@
E->trace.rtrac[j][4][kk]=y;
E->trace.rtrac[j][5][kk]=z;
-
- if (E->trace.itracer_type==1)
- {
- if (rad<=E->trace.z_interface) E->trace.etrac[j][0][kk]=1.0;
- if (rad>E->trace.z_interface) E->trace.etrac[j][0][kk]=0.0;
- }
-
next_try:
;
} /* end while */
@@ -4323,7 +3469,6 @@
int icheck;
int iestimate;
int icushion;
- int ithiscap;
int j;
int icheck_cap();
@@ -4357,10 +3502,6 @@
for (j=1;j<=E->sphere.caps_per_proc;j++)
{
- ithiscap=E->sphere.capid[j];
-
- E->trace.itracsize[j]=iestimate;
-
initialize_tracer_arrays(E,j,iestimate);
for (kk=1;kk<=number_of_tracers;kk++)
@@ -4368,20 +3509,10 @@
fgets(input_s,200,fptracer);
sscanf(input_s,"%lf %lf %lf",&rdum1,&rdum2,&rdum3);
- if (E->trace.icartesian_or_spherical_input==0)
- {
- x=rdum1;
- y=rdum2;
- z=rdum3;
- cart_to_sphere(E,x,y,z,&theta,&phi,&rad);
- }
- if (E->trace.icartesian_or_spherical_input==1)
- {
- theta=rdum1;
- phi=rdum2;
- rad=rdum3;
- sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
- }
+ theta=rdum1;
+ phi=rdum2;
+ rad=rdum3;
+ sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
/* make sure theta, phi is in range, and radius is within bounds */
@@ -4411,15 +3542,6 @@
E->trace.rtrac[j][4][E->trace.itrac[j][0]]=y;
E->trace.rtrac[j][5][E->trace.itrac[j][0]]=z;
- /* Assign composition for active tracers */
-
- if (E->trace.itracer_type==1)
- {
- if (rad<=E->trace.z_interface) E->trace.etrac[j][0][kk]=1.0;
- if (rad>E->trace.z_interface) E->trace.etrac[j][0][kk]=0.0;
- }
-
-
next_tracer:
;
} /* end kk, number of tracers */
@@ -4654,23 +3776,13 @@
double test_point[4];
double rnode[5][10];
+ int lev = E->mesh.levmax;
int ival;
int kk;
int node;
int icheck_bounds();
- /* AKMA REMOVE FOR SPEED */
-
- if ((nel<1) || (nel>E->lmesh.nel) )
- {
- fprintf(E->trace.fpt,"ERROR(icheck element column)-Weird nel: %d\n",nel);
- fprintf(E->trace.fpt,"x: %f y: %f z: %f rad: %f\n",x,y,z,rad);
- fflush(E->trace.fpt);
- exit(10);
- }
-
-
E->trace.istat_elements_checked++;
/* surface coords of element nodes */
@@ -4687,10 +3799,10 @@
rnode[kk][4]=E->sx[j][1][node];
rnode[kk][5]=E->sx[j][2][node];
- rnode[kk][6]=E->trace.DSinCos[j][2][node]; /* cos(theta) */
- rnode[kk][7]=E->trace.DSinCos[j][0][node]; /* sin(theta) */
- rnode[kk][8]=E->trace.DSinCos[j][3][node]; /* cos(phi) */
- rnode[kk][9]=E->trace.DSinCos[j][1][node]; /* sin(phi) */
+ rnode[kk][6]=E->SinCos[lev][j][2][node]; /* cos(theta) */
+ rnode[kk][7]=E->SinCos[lev][j][0][node]; /* sin(theta) */
+ rnode[kk][8]=E->SinCos[lev][j][3][node]; /* cos(phi) */
+ rnode[kk][9]=E->SinCos[lev][j][1][node]; /* sin(phi) */
}
@@ -4958,101 +4070,44 @@
return;
}
-/*******************************************************************/
-/* Function SETUP DSINCOS */
-/* */
-/* open memory for and fill DSinCos */
-/* similar as E.SinCos but double instead of float */
-/* DSinCos[cap][0][node]=sint */
-/* DSinCos[cap][1][node]=sinf */
-/* DSinCos[cap][2][node]=cost */
-/* DSinCos[cap][3][node]=cosf */
-/* */
-void setup_dsincos(E)
- struct All_variables *E;
-{
-
- int j,kk;
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=0;kk<=3;kk++)
- {
- if ((E->trace.DSinCos[j][kk]=(double *)malloc((E->lmesh.nno+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(setup_dsincos)-not enough memory(ib)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- for (kk=1;kk<=(E->lmesh.nno);kk++)
- {
- E->trace.DSinCos[j][0][kk]=sin(E->sx[j][1][kk]);
- E->trace.DSinCos[j][1][kk]=sin(E->sx[j][2][kk]);
- E->trace.DSinCos[j][2][kk]=cos(E->sx[j][1][kk]);
- E->trace.DSinCos[j][3][kk]=cos(E->sx[j][2][kk]);
- }
-
- } /* end each cap */
-
- return;
-}
-
/************ FIX RADIUS ********************************************/
/* This function moves particles back in bounds if they left */
/* during advection */
-void fix_radius(E,radius,theta,phi,x,y,z)
- struct All_variables *E;
- double *radius;
- double *theta;
- double *phi;
- double *x;
- double *y;
- double *z;
-
+static void fix_radius(struct All_variables *E,
+ double *radius, double *theta, double *phi,
+ double *x, double *y, double *z)
{
-
double sint,cost,sinf,cosf,rad;
+ double max_radius, min_radius;
- static int been_here=0;
- static double max_radius;
- static double min_radius;
+ max_radius = E->sphere.ro - E->trace.box_cushion;
+ min_radius = E->sphere.ri + E->trace.box_cushion;
- if (been_here==0)
- {
- max_radius=E->sphere.ro-E->trace.box_cushion;
- min_radius=E->sphere.ri+E->trace.box_cushion;
- been_here++;
- }
+ if (*radius > max_radius) {
+ *radius=max_radius;
+ rad=max_radius;
+ cost=cos(*theta);
+ sint=sqrt(1.0-cost*cost);
+ cosf=cos(*phi);
+ sinf=sin(*phi);
+ *x=rad*sint*cosf;
+ *y=rad*sint*sinf;
+ *z=rad*cost;
+ }
+ if (*radius < min_radius) {
+ *radius=min_radius;
+ rad=min_radius;
+ cost=cos(*theta);
+ sint=sqrt(1.0-cost*cost);
+ cosf=cos(*phi);
+ sinf=sin(*phi);
+ *x=rad*sint*cosf;
+ *y=rad*sint*sinf;
+ *z=rad*cost;
+ }
- if (*radius > max_radius)
- {
- *radius=max_radius;
- rad=max_radius;
- cost=cos(*theta);
- sint=sqrt(1.0-cost*cost);
- cosf=cos(*phi);
- sinf=sin(*phi);
- *x=rad*sint*cosf;
- *y=rad*sint*sinf;
- *z=rad*cost;
- }
- if (*radius < min_radius)
- {
- *radius=min_radius;
- rad=min_radius;
- cost=cos(*theta);
- sint=sqrt(1.0-cost*cost);
- cosf=cos(*phi);
- sinf=sin(*phi);
- *x=rad*sint*cosf;
- *y=rad*sint*sinf;
- *z=rad*cost;
- }
-
-
return;
}
@@ -5064,46 +4119,39 @@
/* between 0 and 2 PI */
/* */
-void fix_phi(phi)
- double *phi;
-
+static void fix_phi(double *phi)
{
+ const double two_pi=2.0*M_PI;
- static double two_pi=2.0*M_PI;
+ double d2 = floor(*phi / two_pi);
+ *phi -= two_pi * d2;
- do_again:
-
- if (*phi<0.0) *phi=*phi+two_pi;
- else if (*phi>=two_pi) *phi=*phi-two_pi;
-
- if ((*phi<0.0) || (*phi >=two_pi)) goto do_again;
-
-
return;
}
/******************************************************************/
-/* FIX THETA */
+/* FIX THETA PHI */
/* */
/* This function constrains the value of theta to be */
-/* between 0 and PI */
+/* between 0 and PI, and */
+/* this function constrains the value of phi to be */
+/* between 0 and 2 PI */
/* */
-void fix_theta(theta)
- double *theta;
+static void fix_theta_phi(double *theta, double *phi)
{
+ const double two_pi=2.0*M_PI;
+ double d, d2;
- static double two_pi=2.0*M_PI;
+ d = floor(*theta/M_PI);
- do_again2:
+ *theta -= M_PI * d;
+ *phi += M_PI * d;
- if (*theta<0.0) *theta=fabs(*theta);
- else if (*theta>M_PI) *theta=two_pi-*theta;
+ d2 = floor(*phi / two_pi);
- if ((*theta<0.0) || (*theta >M_PI)) goto do_again2;
+ *phi -= two_pi * d2;
-
-
return;
}
@@ -5697,6 +4745,7 @@
for (kk=1;kk<=2;kk++)
{
+ //TODO: allocate for surface nodes only to save memory
if ((E->trace.UV[j][kk]=(double *)malloc((numnodes+1)*sizeof(double)))==NULL)
{
fprintf(E->trace.fpt,"Error(define uv)-not enough memory(a)\n");
@@ -5786,6 +4835,7 @@
{
for (kk=1;kk<=9;kk++)
{
+ //TODO: allocate for surface elements only to save memory
if ((E->trace.shape_coefs[j][iwedge][kk]=
(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
{
@@ -5884,6 +4934,7 @@
int j,m,i;
int kk;
+ int lev = E->mesh.levmax;
double sint,sinf,cost,cosf;
double v_theta,v_phi,v_rad;
@@ -5915,10 +4966,10 @@
for (i=1;i<=E->lmesh.nno;i++)
{
- sint=E->trace.DSinCos[m][0][i];
- sinf=E->trace.DSinCos[m][1][i];
- cost=E->trace.DSinCos[m][2][i];
- cosf=E->trace.DSinCos[m][3][i];
+ sint=E->SinCos[lev][m][0][i];
+ sinf=E->SinCos[lev][m][1][i];
+ cost=E->SinCos[lev][m][2][i];
+ cosf=E->SinCos[lev][m][3][i];
v_theta=E->sphere.cap[m].V[1][i];
v_phi=E->sphere.cap[m].V[2][i];
@@ -5953,14 +5004,8 @@
double *phi;
double *rad;
{
-
- void fix_radius();
- void fix_theta();
- void fix_phi();
-
+ fix_theta_phi(theta, phi);
fix_radius(E,rad,theta,phi,x,y,z);
- fix_theta(theta);
- fix_phi(phi);
return;
}
@@ -5975,20 +5020,8 @@
int number,iold_number;
- static int been_here=0;
-
- int isum_tracers();
-
number=isum_tracers(E);
- if (been_here==0)
- {
- E->trace.ilast_tracer_count=number;
- fprintf(E->trace.fpt,"Sum of Tracers: %d\n",number);
- been_here++;
- return;
- }
-
iold_number=E->trace.ilast_tracer_count;
if (number!=iold_number)
@@ -6008,8 +5041,7 @@
/* */
/* This function uses MPI to sum all tracers and returns number of them. */
-int isum_tracers(E)
- struct All_variables *E;
+static int isum_tracers(struct All_variables *E)
{
@@ -6464,60 +5496,6 @@
}
-/******************* READ COMP ***********************************************/
-/* */
-/* This function is similar to read_temp. It is used to read the composition */
-/* from file for post-proceesing. */
-
-void read_comp(E)
- struct All_variables *E;
-{
- int i,ii,m,mm,ll;
- char output_file[255],input_s[1000];
-
- double g;
- FILE *fp;
-
-
-
- ii = E->monitor.solution_cycles;
- sprintf(output_file,"%s.comp.%d.%d",E->control.old_P_file,E->parallel.me,ii);
-
- if ((fp=fopen(output_file,"r"))==NULL)
- {
- fprintf(stderr,"ERROR(read_temp) - %s not found\n",output_file);
- fflush(stderr);
- exit(10);
- }
-
- fgets(input_s,1000,fp);
- sscanf(input_s,"%d %d %f",&ll,&mm,&E->monitor.elapsed_time);
-
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- {
- E->trace.comp_node[m]=(double *)malloc((E->lmesh.nno+1)*sizeof(double));
-
- fgets(input_s,1000,fp);
- sscanf(input_s,"%d %d",&ll,&mm);
- for(i=1;i<=E->lmesh.nno;i++)
- {
- if (fgets(input_s,1000,fp)==NULL)
- {
- fprintf(stderr,"ERROR(read_comp) -data for node %d not found\n",i);
- fflush(stderr);
- exit(10);
- }
- sscanf(input_s,"%lf",&g);
- E->trace.comp_node[m][i] = g;
-
- }
- }
-
- fclose (fp);
- return;
-}
-
-
/******************* get_neighboring_caps ************************************/
/* */
/* Communicate with neighboring processors to get their cap boundaries, */
Modified: mc/3D/CitcomS/trunk/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-03-09 00:25:03 UTC (rev 6206)
@@ -359,7 +359,6 @@
{
void get_global_shape_fn();
- void float_global_operation();
double rtf[4][9];
int n,i,j,el,m;
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-03-09 00:25:03 UTC (rev 6206)
@@ -73,7 +73,6 @@
void set_sphere_harmonics (struct All_variables*);
void set_starting_age(struct All_variables*);
void tracer_initial_settings(struct All_variables*);
-double CPU_time0();
@@ -208,7 +207,7 @@
/* first the problem type (defines subsequent behaviour) */
- input_string("Problem",E->control.PROBLEM_TYPE,NULL,m);
+ input_string("Problem",E->control.PROBLEM_TYPE,"convection",m);
if ( strcmp(E->control.PROBLEM_TYPE,"convection") == 0) {
E->control.CONVECTION = 1;
set_convection_defaults(E);
@@ -226,7 +225,7 @@
set_convection_defaults(E);
}
- input_string("Geometry",E->control.GEOMETRY,NULL,m);
+ input_string("Geometry",E->control.GEOMETRY,"sphere",m);
if ( strcmp(E->control.GEOMETRY,"cart2d") == 0)
{ E->control.CART2D = 1;
(E->solver.set_2dc_defaults)(E);}
@@ -247,7 +246,7 @@
E->control.CART2D = 1;
(E->solver.set_2dc_defaults)(E); }
- input_string("Solver",E->control.SOLVER_TYPE,NULL,m);
+ input_string("Solver",E->control.SOLVER_TYPE,"cgrad",m);
if ( strcmp(E->control.SOLVER_TYPE,"cgrad") == 0)
{ E->control.CONJ_GRAD = 1;
set_cg_defaults(E);}
@@ -496,7 +495,7 @@
}
for(d=0;d<=3;d++)
- E->SinCos[i][j][d] = (float *) malloc((nno+1)*sizeof(float));
+ E->SinCos[i][j][d] = (double *) malloc((nno+1)*sizeof(double));
E->IEN[i][j] = (struct IEN *) malloc((nel+2)*sizeof(struct IEN));
E->EL[i][j] = (struct SUBEL *) malloc((nel+2)*sizeof(struct SUBEL));
@@ -999,6 +998,9 @@
E->output.botm = 0;
E->output.geoid = 0;
E->output.horiz_avg = 0;
+ E->output.tracer = 0;
+ E->output.comp_el = 0;
+ E->output.comp_nd = 0;
while(1) {
/* get next field */
@@ -1035,6 +1037,12 @@
}
else if(strcmp(prev, "horiz_avg")==0)
E->output.horiz_avg = 1;
+ else if(strcmp(prev, "tracer")==0)
+ E->output.tracer = 1;
+ else if(strcmp(prev, "comp_el")==0)
+ E->output.comp_el = 1;
+ else if(strcmp(prev, "comp_nd")==0)
+ E->output.comp_nd = 1;
else
if(E->parallel.me == 0)
fprintf(stderr, "Warning: unknown field for output_optional: %s\n", prev);
Modified: mc/3D/CitcomS/trunk/lib/Makefile.am
===================================================================
--- mc/3D/CitcomS/trunk/lib/Makefile.am 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/Makefile.am 2007-03-09 00:25:03 UTC (rev 6206)
@@ -62,6 +62,7 @@
BC_util.c \
Citcom_init.c \
citcom_init.h \
+ Composition_related.c \
Construct_arrays.c \
Convection.c \
convection_variables.h \
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2007-03-09 00:25:03 UTC (rev 6206)
@@ -36,6 +36,8 @@
#include "parsing.h"
#include "output.h"
+void output_comp_nd(struct All_variables *, int);
+void output_comp_el(struct All_variables *, int);
void output_coord(struct All_variables *);
void output_mat(struct All_variables *);
void output_velo(struct All_variables *, int);
@@ -60,7 +62,7 @@
int m = E->parallel.me;
input_string("output_format", E->output.format, "ascii",m);
- input_string("output_optional", E->output.optional, "surf,botm",m);
+ input_string("output_optional", E->output.optional, "surf,botm,tracer,comp_el",m);
}
@@ -79,11 +81,6 @@
output_surf_botm(E, cycles);
- /* output tracer location if using tracer */
- //TODO: output_tracer() is not working for full version of tracers
- //if(E->control.tracer==1)
- //output_tracer(E, cycles);
-
/* optiotnal output below */
/* compute and output geoid (in spherical harmonics coeff) */
if (E->output.geoid == 1)
@@ -98,6 +95,15 @@
if (E->output.horiz_avg == 1)
output_horiz_avg(E, cycles);
+ if(E->output.tracer == 1 && E->control.tracer == 1)
+ output_tracer(E, cycles);
+
+ if (E->output.comp_nd == 1 && E->composition.on)
+ output_comp_nd(E, cycles);
+
+ if (E->output.comp_el == 1 && E->composition.on)
+ output_comp_el(E, cycles);
+
return;
}
@@ -387,7 +393,7 @@
void output_tracer(struct All_variables *E, int cycles)
{
- int n;
+ int j, n, ncolumns;
char output_file[255];
FILE *fp1;
@@ -395,18 +401,124 @@
E->parallel.me, cycles);
fp1 = output_open(output_file);
- fprintf(fp1,"%.5e\n",E->monitor.elapsed_time);
+ //TODO: unify
- for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++) {
- fprintf(fp1,"%.4e %.4e %.4e %.4e\n", E->Tracer.itcolor[n], E->Tracer.tracer_x[n],E->Tracer.tracer_y[n],E->Tracer.tracer_z[n]);
+ if(E->parallel.nprocxy==1) {
+ /* regional model */
+ fprintf(fp1,"%d %d %.5e\n", cycles, E->Tracer.LOCAL_NUM_TRACERS,
+ E->monitor.elapsed_time);
+
+ for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++) {
+ fprintf(fp1,"%.4e %.4e %.4e %.4e\n", E->Tracer.itcolor[n], E->Tracer.tracer_x[n],E->Tracer.tracer_y[n],E->Tracer.tracer_z[n]);
+ }
}
+ else {
+ /* global model */
+ if (E->composition.ichemical_buoyancy==1)
+ ncolumns = 4;
+ else if (E->composition.ichemical_buoyancy==0)
+ ncolumns = 3;
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ fprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.itrac[j][0],
+ ncolumns, E->monitor.elapsed_time);
+
+ if (E->composition.ichemical_buoyancy==1) {
+ for(n=1;n<=E->trace.itrac[j][0];n++) {
+ fprintf(fp1,"%.12e %.12e %.12e %.12e\n",
+ E->trace.rtrac[j][0][n],
+ E->trace.rtrac[j][1][n],
+ E->trace.rtrac[j][2][n],
+ E->trace.etrac[j][0][n]);
+ }
+ }
+ else if (E->composition.ichemical_buoyancy==0) {
+ for(n=1;n<=E->trace.itrac[j][0];n++) {
+ fprintf(fp1,"%.12e %.12e %.12e\n",
+ E->trace.rtrac[j][0][n],
+ E->trace.rtrac[j][1][n],
+ E->trace.rtrac[j][2][n]);
+ }
+ }
+
+ }
+ }
+
fclose(fp1);
-
return;
}
+void output_comp_nd(struct All_variables *E, int cycles)
+{
+ int i, j;
+ char output_file[255];
+ FILE *fp1;
+
+ sprintf(output_file,"%s.comp_nd.%d.%d", E->control.data_file,
+ E->parallel.me, cycles);
+ fp1 = output_open(output_file);
+
+ //TODO: unify
+
+ if(E->parallel.nprocxy==1) {
+ }
+ else {
+ fprintf(fp1,"%d %d %.5e %.5e %.5e\n",
+ cycles, E->lmesh.nno,
+ E->monitor.elapsed_time,
+ E->trace.initial_bulk_composition,
+ E->trace.bulk_composition);
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ fprintf(fp1,"%3d %7d\n", j, E->lmesh.nno);
+ for(i=1;i<=E->lmesh.nno;i++) {
+ fprintf(fp1,"%.6e\n",E->trace.comp_node[j][i]);
+ }
+ }
+
+ }
+
+ fclose(fp1);
+ return;
+}
+
+
+void output_comp_el(struct All_variables *E, int cycles)
+{
+ int i, j;
+ char output_file[255];
+ FILE *fp1;
+
+ sprintf(output_file,"%s.comp_el.%d.%d", E->control.data_file,
+ E->parallel.me, cycles);
+ fp1 = output_open(output_file);
+
+ //TODO: unify
+
+ if(E->parallel.nprocxy==1) {
+ }
+ else {
+ fprintf(fp1,"%d %d %.5e %.5e %.5e\n",
+ cycles, E->lmesh.nel,
+ E->monitor.elapsed_time,
+ E->trace.initial_bulk_composition,
+ E->trace.bulk_composition);
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ fprintf(fp1,"%3d %7d\n", j, E->lmesh.nel);
+ for(i=1;i<=E->lmesh.nel;i++) {
+ fprintf(fp1,"%.6e\n",E->trace.oldel[j][i]);
+ }
+ }
+
+ }
+
+ fclose(fp1);
+ return;
+}
+
+
void output_time(struct All_variables *E, int cycles)
{
double CPU_time0();
Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-03-09 00:25:03 UTC (rev 6206)
@@ -126,10 +126,10 @@
}
-void thermal_buoyancy(struct All_variables *E, double **buoy)
+void get_buoyancy(struct All_variables *E, double **buoy)
{
int i,m;
- double temp,*H;
+ double temp,temp2,*H;
void remove_horiz_ave();
H = (double *)malloc( (E->lmesh.noz+1)*sizeof(double));
@@ -140,6 +140,14 @@
for(i=1;i<=E->lmesh.nno;i++)
buoy[m][i] = temp * E->T[m][i];
+ /* chemical buoyancy */
+ if(E->composition.ichemical_buoyancy==1) {
+ temp2 = E->composition.buoyancy_ratio * temp;
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=E->lmesh.nno;i++)
+ buoy[m][i] -= temp2 * E->trace.comp_node[m][i];
+ }
+
phase_change_apply_410(E, buoy);
phase_change_apply_670(E, buoy);
phase_change_apply_cmb(E, buoy);
Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-03-09 00:25:03 UTC (rev 6206)
@@ -420,7 +420,7 @@
{
int a, node;
- float sint, cost, sinf, cosf;
+ double sint, cost, sinf, cosf;
const int dims=E->mesh.nsd;
const int ends=enodes[E->mesh.nsd];
const int nno=E->lmesh.nno;
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-03-09 00:25:03 UTC (rev 6206)
@@ -99,7 +99,7 @@
{
int m = E->parallel.me;
- input_string("Viscosity",E->viscosity.STRUCTURE,NULL,m);
+ input_string("Viscosity",E->viscosity.STRUCTURE,"system",m);
input_int ("visc_smooth_method",&(E->viscosity.smooth_cycles),"0",m);
if ( strcmp(E->viscosity.STRUCTURE,"system") == 0)
Added: mc/3D/CitcomS/trunk/lib/composition_related.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/composition_related.h 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/composition_related.h 2007-03-09 00:25:03 UTC (rev 6206)
@@ -0,0 +1,33 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ *
+ *<LicenseText>
+ *
+ * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
+ * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
+ * Copyright (C) 1994-2005, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+
+
+void composition_input(struct All_variables *E);
+void composition_setup(struct All_variables *E);
+void write_composition_instructions(struct All_variables *E);
+void init_tracer_composition(struct All_variables *E);
+void fill_composition(struct All_variables *E);
+
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-03-09 00:25:03 UTC (rev 6206)
@@ -493,6 +493,7 @@
int read_slabgeoid;
int pseudo_free_surf;
+
int tracer;
double theta_min, theta_max, fi_min, fi_max;
@@ -654,9 +655,38 @@
int botm; /* whether to output bottom data */
int geoid; /* whether to output geoid/topo spherial harmonics */
int horiz_avg; /* whether to output horizontal averaged profile */
+ int tracer; /* whether to output tracer coordinate */
+ int comp_el; /* whether to output composition at elements */
+ int comp_nd; /* whether to output composition at nodes */
+
};
+struct COMPOSITION {
+
+ int ichemical_buoyancy;
+ int icompositional_rheology;
+
+ /* if any of the above flags is true, turn this flag on */
+ int on;
+
+ double buoyancy_ratio;
+ int ireset_initial_composition;
+ int ibuoy_type;
+ int *ieltrac[13];
+ double *celtrac[13];
+ double *comp_el[13];
+ double *comp_node[13];
+ double initial_bulk_composition;
+ double bulk_composition;
+ double error_fraction;
+
+ double z_interface;
+
+ double compositional_rheology_prefactor;
+};
+
+
#ifdef USE_HDF5
#include "hdf5_related.h"
#endif
@@ -685,15 +715,18 @@
struct SLICE slice;
struct Parallel parallel;
struct SPHERE sphere;
+ struct Bdry boundary;
+ struct SBC sbc;
+ struct Output output;
+
/* for regional tracer*/
struct Tracer Tracer;
+
/* for global tracer*/
struct TRACE trace;
- struct Bdry boundary;
- struct SBC sbc;
- struct Output output;
- int filed[20];
+ /* for chemical convection & composition rheology */
+ struct COMPOSITION composition;
struct COORD *eco[NCS];
struct IEN *ien[NCS]; /* global */
@@ -733,7 +766,7 @@
double *SX[MAX_LEVELS][NCS][4],*X[MAX_LEVELS][NCS][4];
double *sx[NCS][4],*x[NCS][4];
double *surf_det[NCS][5];
- float *SinCos[MAX_LEVELS][NCS][4];
+ double *SinCos[MAX_LEVELS][NCS][4];
float *TT;
float *V[NCS][4],*GV[NCS][4],*GV1[NCS][4];
Modified: mc/3D/CitcomS/trunk/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-03-09 00:25:03 UTC (rev 6206)
@@ -45,7 +45,7 @@
float *FI_LOC_ELEM_T;
float *R_LOC_ELEM_T;
-} r;
+};
struct TRACE{
@@ -63,24 +63,17 @@
int itracer_warnings;
int itracing;
-int itracer_type;
int ianalytical_tracer_test;
-int icompositional_rheology;
-int ibuoy_type;
-int itracer_restart;
-int icartesian_or_spherical_input;
+int ic_method;
int itperel;
int itracer_advection_scheme;
int itracer_interpolation_scheme;
-int iwrite_tracers_every;
-int ireset_initial_composition;
+int icompositional_rheology;
+
int itracsize[13];
int *itrac[13];
-double buoyancy_ratio;
-double z_interface;
-double compositional_rheology_prefactor;
double box_cushion;
double *rtrac[13][100];
double *etrac[13][100];
@@ -122,10 +115,6 @@
int istat1;
int istat_elements_checked;
int ilast_tracer_count;
-double stat_time_start_trace;
-double stat_time_end_trace;
-double stat_trace_time;
-double stat_stokes_time;
/* Mesh information */
double xcap[13][5];
@@ -153,8 +142,6 @@
double sin_theta_f;
double *shape_coefs[13][3][10];
-double *DSinCos[13][4];
-
double *V0_cart[13][4];
double initial_bulk_composition;
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2007-03-08 22:48:41 UTC (rev 6205)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2007-03-09 00:25:03 UTC (rev 6206)
@@ -577,14 +577,32 @@
getStringProperty(properties, "tracer_file", E->control.tracer_file, fp);
}
else if(E->parallel.nprocxy == 12) {
- getIntProperty(properties, "tracer_restart", E->trace.itracer_restart, fp);
- getIntProperty(properties, "tracers_per_element", E->trace.itperel, fp);
- getDoubleProperty(properties, "z_interface", E->trace.z_interface, fp);
+ if(E->control.restart)
+ E->trace.ic_method = 2;
+ else {
- getStringProperty(properties, "tracer_file", E->trace.tracer_file, fp);
- getIntProperty(properties, "cartesian_or_spherical_input", E->trace.icartesian_or_spherical_input, fp);
+ getIntProperty(properties, "tracer_ic_method",
+ E->trace.ic_method, fp);
+ if (E->trace.ic_method==0)
+ getIntProperty(properties, "tracers_per_element",
+ E->trace.itperel, fp);
+ else if (E->trace.ic_method==1)
+ getStringProperty(properties, "tracer_file",
+ E->trace.tracer_file, fp);
+ else if (E->trace.ic_method==2) {
+ }
+ else {
+ fprintf(stderr,"Sorry, tracer_ic_method only 0, 1 and 2 available\n");
+ fflush(stderr);
+ parallel_process_termination();
+ }
+ }
+
+
+ getDoubleProperty(properties, "z_interface", E->composition.z_interface, fp);
+
getIntProperty(properties, "tracer_advection_scheme", E->trace.itracer_advection_scheme, fp);
getIntProperty(properties, "tracer_interpolation_scheme", E->trace.itracer_interpolation_scheme, fp);
@@ -594,15 +612,14 @@
E->trace.delphi[0] = tmp;
getIntProperty(properties, "analytical_tracer_test", E->trace.ianalytical_tracer_test, fp);
- getIntProperty(properties, "tracer_type", E->trace.itracer_type, fp);
- getIntProperty(properties, "buoy_type", E->trace.ibuoy_type, fp);
- getDoubleProperty(properties, "buoyancy_ratio", E->trace.buoyancy_ratio, fp);
- getIntProperty(properties, "reset_initial_composition", E->trace.ireset_initial_composition, fp);
+ getIntProperty(properties, "chemical_buoyancy",
+ E->composition.ichemical_buoyancy, fp);
+ getIntProperty(properties, "buoy_type", E->composition.ibuoy_type, fp);
+ getDoubleProperty(properties, "buoyancy_ratio", E->composition.buoyancy_ratio, fp);
+ getIntProperty(properties, "reset_initial_composition", E->composition.ireset_initial_composition, fp);
getIntProperty(properties, "compositional_rheology", E->trace.icompositional_rheology, fp);
- getDoubleProperty(properties, "compositional_prefactor", E->trace.compositional_rheology_prefactor, fp);
+ getDoubleProperty(properties, "compositional_prefactor", E->composition.compositional_rheology_prefactor, fp);
- getIntProperty(properties, "write_tracers_every", E->trace.iwrite_tracers_every, fp);
-
}
PUTS(("\n"));
More information about the cig-commits
mailing list