[cig-commits] r7837 - in mc/3D/CitcomS/trunk: CitcomS CitcomS/Components CitcomS/Solver bin examples/Cookbook5 lib module tests tests/checkpoint

tan2 at geodynamics.org tan2 at geodynamics.org
Thu Aug 16 14:36:20 PDT 2007


Author: tan2
Date: 2007-08-16 14:36:18 -0700 (Thu, 16 Aug 2007)
New Revision: 7837

Added:
   mc/3D/CitcomS/trunk/lib/Checkpoints.c
   mc/3D/CitcomS/trunk/lib/checkpoints.h
   mc/3D/CitcomS/trunk/tests/checkpoint/
   mc/3D/CitcomS/trunk/tests/checkpoint/restart.cfg
   mc/3D/CitcomS/trunk/tests/checkpoint/run.sh
   mc/3D/CitcomS/trunk/tests/checkpoint/tracer.cfg
Modified:
   mc/3D/CitcomS/trunk/CitcomS/Components/IC.py
   mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py
   mc/3D/CitcomS/trunk/CitcomS/Controller.py
   mc/3D/CitcomS/trunk/CitcomS/Solver/CoupledSolver.py
   mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py
   mc/3D/CitcomS/trunk/bin/Citcom.c
   mc/3D/CitcomS/trunk/examples/Cookbook5/cookbook5.cfg
   mc/3D/CitcomS/trunk/lib/Composition_related.c
   mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
   mc/3D/CitcomS/trunk/lib/Initial_temperature.c
   mc/3D/CitcomS/trunk/lib/Instructions.c
   mc/3D/CitcomS/trunk/lib/Lith_age.c
   mc/3D/CitcomS/trunk/lib/Makefile.am
   mc/3D/CitcomS/trunk/lib/Output.c
   mc/3D/CitcomS/trunk/lib/Problem_related.c
   mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c
   mc/3D/CitcomS/trunk/lib/Tracer_setup.c
   mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
   mc/3D/CitcomS/trunk/lib/global_defs.h
   mc/3D/CitcomS/trunk/lib/initial_temperature.h
   mc/3D/CitcomS/trunk/lib/lith_age.h
   mc/3D/CitcomS/trunk/module/bindings.c
   mc/3D/CitcomS/trunk/module/initial_conditions.c
   mc/3D/CitcomS/trunk/module/initial_conditions.h
   mc/3D/CitcomS/trunk/module/outputs.c
   mc/3D/CitcomS/trunk/module/outputs.h
   mc/3D/CitcomS/trunk/module/setProperties.c
Log:
Merging the checkpoint branch to trunk.

The checkpoint file is called [datafile].chkpt.[rank].[step] and is outputed with 
an interval of controller.checkpointFrequency.

To read back the checkpoint file, set solver.ic.restart=on. The old checkpoint
file, [datafile_old].chkpt.[rank].[solution_cycles_init] will be read in. 

The old way of restarting (read in temperature from the velo files) will still
be available through solver.ic.tic_method=-1.

Setting solver.tracer.ic_method=2 will read the old *.tracer.* files as before,
but the composition field will be calculated according to the tracers, not read
from the *.comp_el.* files. The parameter solver.tracer.reset_initial_composition 
becomes obsolete.


Modified: mc/3D/CitcomS/trunk/CitcomS/Components/IC.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/IC.py	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/IC.py	2007-08-16 21:36:18 UTC (rev 7837)
@@ -54,6 +54,8 @@
 
 
     def launch(self):
+        self.initMaterial()
+        self.initTracer()
         self.initTemperature()
         self.initPressure()
         self.initVelocity()
@@ -61,13 +63,21 @@
         return
 
 
+    def initMaterial(self):
+        from CitcomSLib import initialize_material
+        initialize_material(self.all_variables)
+        return
 
+
+    def initTracer(self):
+        from CitcomSLib import init_tracer_composition
+        init_tracer_composition(self.all_variables)
+        return
+
+
     def initTemperature(self):
-        from CitcomSLib import constructTemperature, restartTemperature
-        if self.inventory.restart:
-            restartTemperature(self.all_variables)
-        else:
-            constructTemperature(self.all_variables)
+        from CitcomSLib import constructTemperature
+        constructTemperature(self.all_variables)
         return
 
 
@@ -106,7 +116,7 @@
         zero_elapsed_time = pyre.inventory.bool("zero_elapsed_time", default=True)
 
         tic_method = pyre.inventory.int("tic_method", default=0,
-                            validator=pyre.inventory.choice([0, 1, 2, 3]))
+                            validator=pyre.inventory.choice([-1, 0, 1, 2, 3]))
 
         # for tic_method=0 or 3
         num_perturbations = pyre.inventory.int("num_perturbations", default=1,

Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py	2007-08-16 21:36:18 UTC (rev 7837)
@@ -62,7 +62,7 @@
 
         # 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=2 (each proc reads its own file)
         tracer_ic_method = inv.int("tracer_ic_method", default=0)
 
         # (tracer_ic_method == 0)
@@ -94,6 +94,8 @@
         # ibuoy_type=1 (ratio method)
         buoy_type = inv.int("buoy_type", default=1)
         buoyancy_ratio = inv.float("buoyancy_ratio", default=1.0)
+
+        # This is not used anymore and is left here for backward compatibility
         reset_initial_composition = inv.bool("reset_initial_composition",
                                              default=False)
 

Modified: mc/3D/CitcomS/trunk/CitcomS/Controller.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Controller.py	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/CitcomS/Controller.py	2007-08-16 21:36:18 UTC (rev 7837)
@@ -67,6 +67,8 @@
 
         # do io for 0th step
         self.save()
+
+        ### XXX: if stokes: advection tracers and terminate
         return
 
 
@@ -145,6 +147,7 @@
 
     def save(self):
         self.solver.save(self.inventory.monitoringFrequency)
+        self.solver.checkpoint(self.inventory.checkpointFrequency)
         return
 
 
@@ -154,3 +157,4 @@
         import pyre.inventory
 
         monitoringFrequency = pyre.inventory.int("monitoringFrequency", default=100)
+        checkpointFrequency = pyre.inventory.int("checkpointFrequency", default=100)

Modified: mc/3D/CitcomS/trunk/CitcomS/Solver/CoupledSolver.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Solver/CoupledSolver.py	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/CitcomS/Solver/CoupledSolver.py	2007-08-16 21:36:18 UTC (rev 7837)
@@ -58,8 +58,11 @@
 
 
     def launch(self, application):
-        self._setup()
 
+        #TODO: checkpoint doesn't contain coupler information yet
+        if self.restart:
+            pass
+
         self.coupler.launch(self)
 
         ic = self.inventory.ic
@@ -86,7 +89,7 @@
 
 
     def advectTracers(self):
-        # override Solver.advectTracers to do nothing, since tracer module
+        # override Solver.advectTracers, since tracer module
         # doesn't work in coupled run
         return
 
@@ -139,9 +142,16 @@
         return
 
 
+    def checkpoint(self, checkpointFrequency):
+        Solver.checkpoint(self, checkpointFrequency)
 
+        if not (self.step % checkpointFrequency):
+            #TODO: checkpoint for coupler
+            pass
+        return
 
 
+
 # version
 __id__ = "$Id$"
 

Modified: mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py	2007-08-16 21:36:18 UTC (rev 7837)
@@ -26,7 +26,7 @@
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 #
 
-from CitcomSLib import CPU_time, output, output_time, return_dt, return_t, return_step
+from CitcomSLib import CPU_time, output, output_time, output_checkpoint, return_dt, return_t, return_step
 from pyre.components.Component import Component
 import journal
 
@@ -110,24 +110,30 @@
             print >> stream, "[CitcomS.controller]"
             print >> stream, ("monitoringFrequency=%d" %
                 application.controller.inventory.monitoringFrequency)
+            print >> stream, ("checkpointFrequency=%d" %
+                application.controller.inventory.checkpointFrequency)
             print >> stream
 
         self.setProperties(stream)
 
-        self.restart = self.inventory.ic.inventory.restart
+        self._setup()
 
-        self.ic_initTemperature = self.inventory.ic.initTemperature
         return
 
 
     def launch(self, application):
-        self._setup()
+        if self.inventory.ic.inventory.restart:
+            from CitcomSLib import readCheckpoint
+            readCheckpoint(self.all_variables)
 
-        # initial conditions
-        ic = self.inventory.ic
-        ic.launch()
+            # XXX: if post_processing
+            # calling post_processing() and terminate
+        else:
+            # initial conditions
+            ic = self.inventory.ic
+            ic.launch()
 
-        self.solveVelocities()
+            self.solveVelocities()
         return
 
 
@@ -224,6 +230,14 @@
         return
 
 
+    def checkpoint(self, checkpointFrequency):
+        step = self.step
+
+        if not (step % checkpointFrequency):
+            output_checkpoint(self.all_variables)
+        return
+
+
     def setProperties(self, stream):
 
         from CitcomSLib import Solver_set_properties

Modified: mc/3D/CitcomS/trunk/bin/Citcom.c
===================================================================
--- mc/3D/CitcomS/trunk/bin/Citcom.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/bin/Citcom.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -35,6 +35,7 @@
 #include "citcom_init.h"
 #include "output.h"
 #include "parallel_related.h"
+#include "checkpoints.h"
 
 extern int Emergency_stop;
 
@@ -48,6 +49,8 @@
   void general_stokes_solver();
   void general_stokes_solver_pseudo_surf();
   void read_instructions();
+  void initial_setup();
+  void initial_conditions();
   void solve_constrained_flow();
   void solve_derived_velocities();
   void process_temp_field();
@@ -84,6 +87,7 @@
 
   start_time = time = CPU_time0();
   read_instructions(E, argv[1]);
+  initial_setup(E);
 
   cpu_time_on_vp_it = CPU_time0();
   initial_time = cpu_time_on_vp_it - time;
@@ -94,19 +98,27 @@
     fflush(E->fp);
   }
 
-  if (E->control.post_p)   {
-    post_processing(E);
-    parallel_process_termination();
+  /* This if-block is replaced by CitcomS.Solver.launch()*/
+  if (E->control.restart || E->control.post_p) {
+      read_checkpoint(E);
+
+      if (E->control.post_p) {
+          post_processing(E);
+          parallel_process_termination();
+      }
   }
+  else {
+      initial_conditions(E);
 
-  if(E->control.pseudo_free_surf) {
-    if(E->mesh.topvbc == 2)
-	    general_stokes_solver_pseudo_surf(E);
-    else
-	    assert(0);
+      if(E->control.pseudo_free_surf) {
+          if(E->mesh.topvbc == 2)
+              general_stokes_solver_pseudo_surf(E);
+          else
+              assert(0);
+      }
+      else
+          general_stokes_solver(E);
   }
-  else
-    general_stokes_solver(E);
 
   (E->problem_output)(E, E->monitor.solution_cycles);
 
@@ -122,8 +134,6 @@
     parallel_process_termination();
   }
 
-  (E->next_buoyancy_field_init)(E);
-
   while ( E->control.keep_going   &&  (Emergency_stop == 0) )   {
 
     /* The next few lines of code were replaced by
@@ -167,6 +177,10 @@
     /* information about simulation time and wall clock time */
     output_time(E, E->monitor.solution_cycles);
 
+    if ((E->monitor.solution_cycles % E->control.checkpoint_frequency)==0) {
+	output_checkpoint(E);
+    }
+
     if(E->control.mat_control==1)
       read_mat_from_file(E);
     /*

Modified: mc/3D/CitcomS/trunk/examples/Cookbook5/cookbook5.cfg
===================================================================
--- mc/3D/CitcomS/trunk/examples/Cookbook5/cookbook5.cfg	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/examples/Cookbook5/cookbook5.cfg	2007-08-16 21:36:18 UTC (rev 7837)
@@ -24,8 +24,8 @@
 vel_bound_file = ./velocity/bvel.dat
 
 [CitcomS.solver.ic]
-restart = on
 solution_cycles_init = 0
+tic_method = -1
 
 [CitcomS.solver.mesher]
 coor = on

Added: mc/3D/CitcomS/trunk/lib/Checkpoints.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Checkpoints.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/Checkpoints.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -0,0 +1,483 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *<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 <sys/file.h>
+#include "element_definitions.h"
+#include "global_defs.h"
+
+/* Private function prototypes */
+static void backup_file(const char *output_file);
+static void write_sentinel(FILE *fp);
+static void read_sentinel(FILE *fp, int me);
+static void general_checkpoint(struct All_variables *E, FILE *fp);
+
+static void general_checkpoint(struct All_variables *E, FILE *fp);
+static void tracer_checkpoint(struct All_variables *E, FILE *fp);
+static void composition_checkpoint(struct All_variables *E, FILE *fp);
+static void energy_checkpoint(struct All_variables *E, FILE *fp);
+static void momentum_checkpoint(struct All_variables *E, FILE *fp);
+
+static void read_general_checkpoint(struct All_variables *E, FILE *fp);
+static void read_tracer_checkpoint(struct All_variables *E, FILE *fp);
+static void read_composition_checkpoint(struct All_variables *E, FILE *fp);
+static void read_energy_checkpoint(struct All_variables *E, FILE *fp);
+static void read_momentum_checkpoint(struct All_variables *E, FILE *fp);
+
+
+void output_checkpoint(struct All_variables *E)
+{
+    char output_file[255];
+    FILE *fp1;
+
+    sprintf(output_file, "%s.chkpt.%d.%d", E->control.data_file,
+            E->parallel.me, E->monitor.solution_cycles);
+
+    /* Disable the backup since the filename is unique. */
+    /* backup_file(output_file); */
+
+    fp1 = fopen(output_file, "wb");
+
+    /* checkpoint for general information */
+    /* this must be the first to be checkpointed */
+    general_checkpoint(E, fp1);
+
+    /* checkpoint for tracer/composition */
+    if(E->control.tracer) {
+        tracer_checkpoint(E, fp1);
+
+        if(E->composition.on)
+            composition_checkpoint(E, fp1);
+    }
+
+    /* checkpoint for energy equation */
+    energy_checkpoint(E, fp1);
+
+    /* checkpoint for momentum equation */
+    momentum_checkpoint(E, fp1);
+
+    fclose(fp1);
+    return;
+}
+
+
+void read_checkpoint(struct All_variables *E)
+{
+    void initialize_material(struct All_variables *E);
+
+    char output_file[255];
+    FILE *fp;
+
+    /* open the checkpoint file */
+    snprintf(output_file, 254, "%s.chkpt.%d.%d", E->control.old_P_file,
+             E->parallel.me, E->monitor.solution_cycles_init);
+    fp = fopen(output_file, "rb");
+    if(fp == NULL) {
+        fprintf(stderr, "Cannot open file: %s\n", output_file);
+        exit(-1);
+    }
+
+    /* check mesh information in the checkpoint file */
+    read_general_checkpoint(E, fp);
+
+    /* init E->mat */
+    initialize_material(E);
+
+    /* read tracer/composition information in the checkpoint file */
+    if(E->control.tracer) {
+        read_tracer_checkpoint(E, fp);
+
+        if(E->composition.on)
+            read_composition_checkpoint(E, fp);
+    }
+
+    /* read energy information in the checkpoint file */
+    read_energy_checkpoint(E, fp);
+
+    /* read momentum information in the checkpoint file */
+    read_momentum_checkpoint(E, fp);
+
+    fclose(fp);
+    return;
+}
+
+
+static void backup_file(const char *output_file)
+{
+    char bak_file[255];
+    int ierr;
+
+    /* check the existence of output_file */
+    if(access(output_file, F_OK) == 0) {
+        /* if exist, renamed it to back up */
+        sprintf(bak_file, "%s.bak", output_file);
+        ierr = rename(output_file, bak_file);
+        if(ierr != 0) {
+            fprintf(stderr, "Warning, cannot backup checkpoint files\n");
+        }
+    }
+
+    return;
+}
+
+
+static void write_sentinel(FILE *fp)
+{
+    int a[4] = {0, 0, 0, 0};
+
+    fwrite(a, sizeof(int), 4, fp);
+}
+
+
+static void read_sentinel(FILE *fp, int me)
+{
+    int i, a[4];
+    char nonzero = 0;
+
+    fread(a, sizeof(int), 4, fp);
+
+    /* check whether a[i] are all zero */
+    for(i=0; i<4; i++)
+        nonzero |= a[i];
+
+    if(nonzero) {
+        fprintf(stderr, "Error in reading checkpoint file: wrong sentinel, "
+                "me=%d\n", me);
+        exit(-1);
+    }
+
+    return;
+}
+
+
+static void general_checkpoint(struct All_variables *E, FILE *fp)
+{
+    /* write mesh information */
+    fwrite(&(E->lmesh.nox), sizeof(int), 1, fp);
+    fwrite(&(E->lmesh.noy), sizeof(int), 1, fp);
+    fwrite(&(E->lmesh.noz), sizeof(int), 1, fp);
+    fwrite(&(E->parallel.nprocx), sizeof(int), 1, fp);
+    fwrite(&(E->parallel.nprocy), sizeof(int), 1, fp);
+    fwrite(&(E->parallel.nprocz), sizeof(int), 1, fp);
+    fwrite(&(E->sphere.caps_per_proc), sizeof(int), 1, fp);
+
+    /* write timing information */
+    fwrite(&(E->monitor.solution_cycles), sizeof(int), 1, fp);
+    fwrite(&(E->monitor.elapsed_time), sizeof(float), 1, fp);
+    fwrite(&(E->advection.timestep), sizeof(float), 1, fp);
+    fwrite(&(E->control.start_age), sizeof(float), 1, fp);
+
+    return;
+}
+
+
+static void read_general_checkpoint(struct All_variables *E, FILE *fp)
+{
+    int tmp[7];
+    double dtmp;
+
+    /* read mesh information */
+    fread(tmp, sizeof(int), 7, fp);
+
+    if((tmp[0] != E->lmesh.nox) ||
+       (tmp[1] != E->lmesh.noy) ||
+       (tmp[2] != E->lmesh.noz) ||
+       (tmp[3] != E->parallel.nprocx) ||
+       (tmp[4] != E->parallel.nprocy) ||
+       (tmp[5] != E->parallel.nprocz) ||
+       (tmp[6] != E->sphere.caps_per_proc)) {
+
+        fprintf(stderr, "Error in reading checkpoint file: mesh, me=%d\n",
+                E->parallel.me);
+        fprintf(stderr, "%d %d %d %d %d %d %d\n",
+                tmp[0], tmp[1], tmp[2], tmp[3],
+                tmp[4], tmp[5], tmp[6]);
+        exit(-1);
+    }
+
+    /* read timing information */
+    fread(&(E->monitor.solution_cycles), sizeof(int), 1, fp);
+    fread(&(E->monitor.elapsed_time), sizeof(float), 1, fp);
+    fread(&(E->advection.timestep), sizeof(float), 1, fp);
+    fread(&(E->control.start_age), sizeof(float), 1, fp);
+
+    E->advection.timesteps = E->monitor.solution_cycles;
+
+    return;
+}
+
+
+static void tracer_checkpoint(struct All_variables *E, FILE *fp)
+{
+    int m, i;
+
+    write_sentinel(fp);
+
+    fwrite(&(E->trace.number_of_basic_quantities), sizeof(int), 1, fp);
+    fwrite(&(E->trace.number_of_extra_quantities), sizeof(int), 1, fp);
+    fwrite(&(E->trace.nflavors), sizeof(int), 1, fp);
+    fwrite(&(E->trace.ilast_tracer_count), sizeof(int), 1, fp);
+
+    for(m=1; m<=E->sphere.caps_per_proc; m++)
+        fwrite(&(E->trace.ntracers[m]), sizeof(int), 1, fp);
+
+    for(m=1; m<=E->sphere.caps_per_proc; m++) {
+        for(i=0; i<E->trace.number_of_basic_quantities; i++) {
+            fwrite(E->trace.basicq[m][i], sizeof(double),
+                   E->trace.ntracers[m]+1, fp);
+        }
+        for(i=0; i<E->trace.number_of_extra_quantities; i++) {
+            fwrite(E->trace.extraq[m][i], sizeof(double),
+                   E->trace.ntracers[m]+1, fp);
+        }
+        fwrite(E->trace.ielement[m], sizeof(int),
+               E->trace.ntracers[m]+1, fp);
+        for(i=0; i<E->trace.nflavors; i++) {
+            fwrite(E->trace.ntracer_flavor[m][i], sizeof(double),
+                   E->lmesh.nel+1, fp);
+        }
+    }
+
+    return;
+}
+
+
+static void read_tracer_checkpoint(struct All_variables *E, FILE *fp)
+{
+    void allocate_tracer_arrays();
+
+    int m, i, itmp;
+
+    read_sentinel(fp, E->parallel.me);
+
+    fread(&itmp, sizeof(int), 1, fp);
+    if (itmp != E->trace.number_of_basic_quantities) {
+        fprintf(stderr, "Error in reading checkpoint file: tracer basicq, me=%d\n",
+                E->parallel.me);
+        fprintf(stderr, "%d\n", itmp);
+        exit(-1);
+
+    }
+
+    fread(&itmp, sizeof(int), 1, fp);
+    if (itmp != E->trace.number_of_extra_quantities) {
+        fprintf(stderr, "Error in reading checkpoint file: tracer extraq, me=%d\n",
+                E->parallel.me);
+        fprintf(stderr, "%d\n", itmp);
+        exit(-1);
+
+    }
+
+    fread(&itmp, sizeof(int), 1, fp);
+    if (itmp != E->trace.nflavors) {
+        fprintf(stderr, "Error in reading checkpoint file: tracer nflavors, me=%d\n",
+                E->parallel.me);
+        fprintf(stderr, "%d\n", itmp);
+        exit(-1);
+
+    }
+
+    fread(&itmp, sizeof(int), 1, fp);
+    E->trace.ilast_tracer_count = itmp;
+
+    /* # of tracers, allocate memory */
+    for(m=1; m<=E->sphere.caps_per_proc; m++) {
+        fread(&itmp, sizeof(int), 1, fp);
+        allocate_tracer_arrays(E, m, itmp);
+        E->trace.ntracers[m] = itmp;
+    }
+
+    /* read tracer data */
+    for(m=1; m<=E->sphere.caps_per_proc; m++) {
+        for(i=0; i<E->trace.number_of_basic_quantities; i++) {
+            fread(E->trace.basicq[m][i], sizeof(double),
+                  E->trace.ntracers[m]+1, fp);
+        }
+        for(i=0; i<E->trace.number_of_extra_quantities; i++) {
+            fread(E->trace.extraq[m][i], sizeof(double),
+                  E->trace.ntracers[m]+1, fp);
+        }
+        fread(E->trace.ielement[m], sizeof(int),
+              E->trace.ntracers[m]+1, fp);
+        for(i=0; i<E->trace.nflavors; i++) {
+            fread(E->trace.ntracer_flavor[m][i], sizeof(double),
+                  E->lmesh.nel+1, fp);
+        }
+    }
+
+    return;
+}
+
+
+static void composition_checkpoint(struct All_variables *E, FILE *fp)
+{
+    int m;
+
+    write_sentinel(fp);
+
+    fwrite(&(E->composition.bulk_composition), sizeof(double), 1, fp);
+    fwrite(&(E->composition.initial_bulk_composition), sizeof(double), 1, fp);
+
+    for(m=1; m<=E->sphere.caps_per_proc; m++) {
+        fwrite(E->composition.comp_el[m], sizeof(double),
+               E->lmesh.nel+1, fp);
+    }
+    for(m=1; m<=E->sphere.caps_per_proc; m++) {
+        fwrite(E->composition.comp_node[m], sizeof(double),
+               E->lmesh.nno+1, fp);
+    }
+
+    return;
+}
+
+
+static void read_composition_checkpoint(struct All_variables *E, FILE *fp)
+{
+    double tmp;
+    int m;
+
+    read_sentinel(fp, E->parallel.me);
+
+    fread(&tmp, sizeof(double), 1, fp);
+    E->composition.bulk_composition = tmp;
+
+    fread(&tmp, sizeof(double), 1, fp);
+    E->composition.initial_bulk_composition = tmp;
+
+    for(m=1; m<=E->sphere.caps_per_proc; m++) {
+        fread(E->composition.comp_el[m], sizeof(double),
+              E->lmesh.nel+1, fp);
+    }
+
+    for(m=1; m<=E->sphere.caps_per_proc; m++) {
+        fread(E->composition.comp_node[m], sizeof(double),
+              E->lmesh.nno+1, fp);
+    }
+
+    /* preventing uninitialized access */
+    E->trace.istat_iempty = 0;
+    E->composition.error_fraction = E->composition.bulk_composition
+        / E->composition.initial_bulk_composition - 1.0;
+
+    return;
+}
+
+
+static void energy_checkpoint(struct All_variables *E, FILE *fp)
+{
+    int m;
+
+    write_sentinel(fp);
+
+    for(m=1; m<=E->sphere.caps_per_proc; m++) {
+        fwrite(E->T[m], sizeof(double), E->lmesh.nno+1, fp);
+        fwrite(E->Tdot[m], sizeof(double), E->lmesh.nno+1, fp);
+    }
+
+    return;
+}
+
+
+static void read_energy_checkpoint(struct All_variables *E, FILE *fp)
+{
+    int m;
+
+    read_sentinel(fp, E->parallel.me);
+
+    for(m=1; m<=E->sphere.caps_per_proc; m++) {
+        fread(E->T[m], sizeof(double), E->lmesh.nno+1, fp);
+        fread(E->Tdot[m], sizeof(double), E->lmesh.nno+1, fp);
+    }
+
+    return;
+}
+
+
+static void momentum_checkpoint(struct All_variables *E, FILE *fp)
+{
+    int m, i;
+    int lev = E->mesh.levmax;
+
+    write_sentinel(fp);
+
+    fwrite(&(E->monitor.vdotv), sizeof(float), 1, fp);
+    fwrite(&(E->monitor.incompressibility), sizeof(float), 1, fp);
+
+    for(m=1; m<=E->sphere.caps_per_proc; m++) {
+        /* Pressure at equation points and nodes */
+        /* Writing E->NP instead of calling p_to_nodes() because p_to_nodes()
+         * contains parallel communication */
+        fwrite(E->P[m], sizeof(double), E->lmesh.npno+1, fp);
+        fwrite(E->NP[m], sizeof(float), E->lmesh.nno+1, fp);
+
+        /* velocity at equation points */
+        fwrite(E->U[m], sizeof(double), E->lmesh.neq+2, fp);
+
+        /* viscosity at quadrature points and node points */
+        fwrite(E->EVI[lev][m], sizeof(float),
+               (E->lmesh.nel+2)*vpoints[E->mesh.nsd], fp);
+        fwrite(E->VI[lev][m], sizeof(float), E->lmesh.nno+2, fp);
+    }
+
+    return;
+}
+
+
+static void read_momentum_checkpoint(struct All_variables *E, FILE *fp)
+{
+    void v_from_vector();
+
+    int m, i;
+    int lev = E->mesh.levmax;
+
+    read_sentinel(fp, E->parallel.me);
+
+    fread(&(E->monitor.vdotv), sizeof(float), 1, fp);
+    fread(&(E->monitor.incompressibility), sizeof(float), 1, fp);
+
+    for(m=1; m<=E->sphere.caps_per_proc; m++) {
+        /* Pressure at equation points and nodes */
+        fread(E->P[m], sizeof(double), E->lmesh.npno+1, fp);
+        fread(E->NP[m], sizeof(float), E->lmesh.nno+1, fp);
+
+        /* velocity at equation points */
+        fread(E->U[m], sizeof(double), E->lmesh.neq+2, fp);
+
+        /* viscosity at quadrature points and node points */
+        fread(E->EVI[lev][m], sizeof(float),
+              (E->lmesh.nel+2)*vpoints[E->mesh.nsd], fp);
+        fread(E->VI[lev][m], sizeof(float), E->lmesh.nno+2, fp);
+    }
+
+    /* update velocity array */
+    v_from_vector(E);
+
+    return;
+}
+
+
+

Modified: mc/3D/CitcomS/trunk/lib/Composition_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Composition_related.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/Composition_related.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -35,7 +35,6 @@
 
 static void allocate_composition_memory(struct All_variables *E);
 static void compute_elemental_composition_ratio_method(struct All_variables *E);
-static void init_composition(struct All_variables *E);
 static void init_bulk_composition(struct All_variables *E);
 static void check_initial_composition(struct All_variables *E);
 static void map_composition_to_nodes(struct All_variables *E);
@@ -64,9 +63,6 @@
             parallel_process_termination();
         }
 
-        input_int("reset_initial_composition",
-                  &(E->composition.ireset_initial_composition),"0",m);
-
     }
 
 
@@ -93,12 +89,8 @@
 
 void composition_setup(struct All_variables *E)
 {
+    allocate_composition_memory(E);
 
-    if (E->composition.on) {
-        allocate_composition_memory(E);
-        init_composition(E);
-    }
-
     return;
 }
 
@@ -131,12 +123,6 @@
 
         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");
@@ -209,7 +195,7 @@
 }
 
 
-static void init_composition(struct All_variables *E)
+void init_composition(struct All_variables *E)
 {
   if (E->composition.ichemical_buoyancy && 
       E->composition.ibuoy_type) {
@@ -364,22 +350,10 @@
 static void init_bulk_composition(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 idum0, idum1;
 
-
-    FILE *fp;
-
-
     /* ival=0 returns integral not average */
 
     volume = return_bulk_value_d(E,E->composition.comp_node,ival);
@@ -387,24 +361,6 @@
     E->composition.bulk_composition = volume;
     E->composition.initial_bulk_composition = volume;
 
-
-    /* If retarting tracers, the initital bulk composition is read from file */
-    if (E->trace.ic_method == 2 &&
-        !E->composition.ireset_initial_composition) {
-
-        sprintf(output_file,"%s.comp_el.%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",
-               &idum0,&idum1,&rdum1,&rdum2,&rdum3);
-
-        E->composition.initial_bulk_composition = rdum2;
-        fclose(fp);
-
-    }
-
     return;
 }
 

Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -126,7 +126,6 @@
 
     char output_file[200];
     void get_neighboring_caps();
-    void initialize_tracers();
     void analytical_test();
 
     /* Some error control */
@@ -216,17 +215,14 @@
     make_regular_grid(E);
 
 
-    initialize_tracers(E);
-
-
     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);
+    if (E->composition.on)
+        composition_setup(E);
 
     return;
 }
@@ -1861,7 +1857,7 @@
         }
     if (E->trace.ic_method==2)
         {
-            fprintf(E->trace.fpt,"Restarting Tracers\n");
+            fprintf(E->trace.fpt,"Reading individual tracer files\n");
         }
 
     fprintf(E->trace.fpt,"Number of tracer flavors: %d\n", E->trace.nflavors);

Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -38,7 +38,7 @@
 
 #include "initial_temperature.h"
 void debug_tic(struct All_variables *);
-void restart_tic_from_file(struct All_variables *);
+void read_tic_from_file(struct All_variables *);
 
 #ifdef USE_GZDIR
 void restart_tic_from_gzdir_file(struct All_variables *);
@@ -57,6 +57,9 @@
   /* When tic_method is 0 (default), the temperature is a linear profile +
      perturbation at some layers.
 
+     When tic_method is -1, the temperature is read in from the
+     [datafile_old].velo.[rank].[solution_cycles_init] files.
+
      When tic_method is 1, the temperature is isothermal (== bottom b.c.) +
      uniformly cold plate (thickness specified by 'half_space_age').
 
@@ -109,7 +112,8 @@
       E->convection.load_depth[0] = (noz+1)/2;
     }
 
-  } else if (E->convection.tic_method == 1) {
+  }
+  else if (E->convection.tic_method == 1) {
 
     input_float("half_space_age", &(E->convection.half_space_age), "40.0,1e-3,nomax", m);
 
@@ -142,14 +146,25 @@
 
 
 
+/* This function is replaced by CitcomS.Components.IC.initTemperature()*/
 void convection_initial_temperature(struct All_variables *E)
 {
-  if (E->control.restart)
-    restart_tic(E);
-  else if (E->control.post_p)
-    restart_tic(E);
+  void report();
+
+  report(E,"Initialize temperature field");
+
+  if (E->control.lith_age)
+    lith_age_construct_tic(E);
+  else if (E->convection.tic_method == -1) {
+#ifdef USE_GZDIR
+      if(strcmp(E->output.format, "ascii-gz") == 0)
+          restart_tic_from_gzdir_file(E);
+      else
+#endif
+          read_tic_from_file(E);
+  }
   else
-    construct_tic(E);
+    (E->solver.construct_tic_from_input)(E);
 
   /* Note: it is the callee's responsibility to conform tbc. */
   /* like a call to temperatures_conform_bcs(E); */
@@ -161,36 +176,6 @@
 }
 
 
-
-void restart_tic(struct All_variables *E)
-{
-  /*
-  if (E->control.lith_age)
-    lith_age_restart_tic(E);
-  else
-  */
-#ifdef USE_GZDIR
-  if(strcmp(E->output.format, "ascii-gz") == 0)
-    restart_tic_from_gzdir_file(E);
-  else
-#endif
-    restart_tic_from_file(E);
-
-  return;
-}
-
-
-void construct_tic(struct All_variables *E)
-{
-  if (E->control.lith_age)
-    lith_age_construct_tic(E);
-  else
-    (E->solver.construct_tic_from_input)(E);
-
-  return;
-}
-
-
 void debug_tic(struct All_variables *E)
 {
   int m, j;
@@ -208,10 +193,12 @@
 
 
 
-void restart_tic_from_file(struct All_variables *E)
+void read_tic_from_file(struct All_variables *E)
 {
+  void temperatures_conform_bcs();
+
   int ii, ll, mm;
-  float restart_elapsed_time;
+  float tt;
   int i, m;
   char output_file[255], input_s[1000];
   FILE *fp;
@@ -227,10 +214,10 @@
   }
 
   if (E->parallel.me==0)
-    fprintf(E->fp,"Reading %s for restarted temperature\n",output_file);
+    fprintf(E->fp,"Reading %s for initial temperature\n",output_file);
 
   fgets(input_s,1000,fp);
-  sscanf(input_s,"%d %d %f",&ll,&mm,&restart_elapsed_time);
+  sscanf(input_s,"%d %d %f",&ll,&mm,&tt);
 
   for(m=1;m<=E->sphere.caps_per_proc;m++) {
     fgets(input_s,1000,fp);
@@ -239,15 +226,15 @@
       fgets(input_s,1000,fp);
       sscanf(input_s,"%g %g %g %f",&(v1),&(v2),&(v3),&(g));
 
-      /*  E->sphere.cap[m].V[1][i] = d;
-	  E->sphere.cap[m].V[1][i] = e;
-	  E->sphere.cap[m].V[1][i] = f;  */
+      /* Truncate the temperature to be within (0,1). */
+      /* This might not be desirable in some situations. */
       E->T[m][i] = max(0.0,min(g,1.0));
     }
   }
   fclose (fp);
 
   temperatures_conform_bcs(E);
+
   return;
 }
 

Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -61,7 +61,6 @@
 void construct_ien(struct All_variables*);
 void construct_lm(struct All_variables*);
 void construct_masks(struct All_variables*);
-void construct_mat_group(struct All_variables*);
 void construct_shape_functions(struct All_variables*);
 void construct_sub_element(struct All_variables*);
 void construct_surf_det (struct All_variables*);
@@ -71,7 +70,6 @@
 void lith_age_init(struct All_variables *E);
 void mass_matrix(struct All_variables*);
 void output_init(struct All_variables*);
-void read_mat_from_file(struct All_variables*);
 void set_elapsed_time(struct All_variables*);
 void set_sphere_harmonics (struct All_variables*);
 void set_starting_age(struct All_variables*);
@@ -107,7 +105,7 @@
 
     allocate_velocity_vars(E);
 
-    get_initial_elapsed_time(E);  /* Get elapsed time from restart run*/
+    get_initial_elapsed_time(E);  /* Set elapsed time */
     set_starting_age(E);  /* set the starting age to elapsed time, if desired */
     set_elapsed_time(E);         /* reset to elapsed time to zero, if desired */
 
@@ -148,25 +146,13 @@
 	(E->problem_tracer_setup)(E);
     }
 
-    if(E->control.mat_control) {
-      if(E->parallel.me ==0) fprintf(stderr,"IN Instructions.c\n");
-      fflush(stderr);
-      read_mat_from_file(E);
-     }
-    else
-      construct_mat_group(E);
-
 }
 
 
 void read_instructions(struct All_variables *E, char *filename)
 {
-
-    void common_initial_fields();
     void read_initial_settings();
     void global_default_values();
-    void general_stokes_solver_setup();
-    void initial_mesh_solver_setup();
 
     void setup_parser();
     void shutdown_parser();
@@ -187,24 +173,67 @@
     global_default_values(E);
     read_initial_settings(E);
 
+    shutdown_parser(E);
+
+    return;
+}
+
+
+/* This function is replaced by CitcomS.Solver._setup() */
+void initial_setup(struct All_variables *E)
+{
+    void general_stokes_solver_setup();
+    void initial_mesh_solver_setup();
+
     initial_mesh_solver_setup(E);
 
     general_stokes_solver_setup(E);
 
+    (E->next_buoyancy_field_init)(E);
+
+
     if (E->parallel.me==0) fprintf(stderr,"time=%f\n",
                                    CPU_time0()-E->monitor.cpu_time_at_start);
 
+    return;
+}
+
+
+void initialize_material(struct All_variables *E)
+{
+    void construct_mat_group();
+    void read_mat_from_file();
+
+    if(E->control.mat_control)
+        read_mat_from_file(E);
+    else
+        construct_mat_group(E);
+}
+
+
+/* This function is replaced by CitcomS.Components.IC.launch()*/
+void initial_conditions(struct All_variables *E)
+{
+    void initialize_tracers();
+    void init_composition();
+    void common_initial_fields();
+
+    initialize_material(E);
+
+    if (E->control.tracer==1) {
+        initialize_tracers(E);
+
+        if (E->composition.on)
+            init_composition(E);
+    }
+
     (E->problem_initial_fields)(E);   /* temperature/chemistry/melting etc */
     common_initial_fields(E);  /* velocity/pressure/viscosity (viscosity must be done LAST) */
 
-    shutdown_parser(E);
-
     return;
 }
 
 
-
-
 void read_initial_settings(struct All_variables *E)
 {
   void set_convection_defaults();
@@ -372,6 +401,7 @@
   input_float("tole_compressibility",&(E->control.tole_comp),"0.0",m);
 
   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_boolean("precond",&(E->control.precondition),"off",m);
@@ -920,6 +950,7 @@
    ============================================================= */
 
 
+/* This function is replaced by CitcomS.Components.IC.launch()*/
 void common_initial_fields(E)
     struct All_variables *E;
 {
@@ -938,7 +969,8 @@
 
     return;
 
-   }
+}
+
 /* ========================================== */
 
 void initial_pressure(E)
@@ -963,9 +995,6 @@
         E->sphere.cap[m].V[1][i]=0.0;
         E->sphere.cap[m].V[2][i]=0.0;
         E->sphere.cap[m].V[3][i]=0.0;
-        E->sphere.cap[m].Vprev[1][i]=0.0;
-        E->sphere.cap[m].Vprev[2][i]=0.0;
-        E->sphere.cap[m].Vprev[3][i]=0.0;
         }
 
     return;
@@ -1106,7 +1135,8 @@
       parallel_process_termination();
   }
 
-  if (E->control.restart ||
+  if (E->control.restart || E->control.post_p ||
+      E->convection.tic_method == -1 ||
       (E->control.tracer && E->trace.ic_method == 2)) {
       found = strchr(E->control.data_prefix_old, '/');
       if (found) {
@@ -1209,7 +1239,8 @@
     snprintf(E->control.data_file, 200, "%s/%s", E->control.data_dir,
 	     E->control.data_prefix);
 
-    if (E->control.restart ||
+    if (E->control.restart || E->control.post_p ||
+        E->convection.tic_method == -1 ||
         (E->control.tracer && E->trace.ic_method == 2)) {
 	expand_datadir(E, E->control.data_dir_old);
 	snprintf(E->control.old_P_file, 200, "%s/%s", E->control.data_dir_old,
@@ -1264,7 +1295,7 @@
   if (E->fp_out)
     fclose(E->fp_out);
 
-  /* 
+  /*
      remove VTK geo file in case we used that for IO (we'll only do
      this for one processor, since this IO option requires shared
      filesystems anyway

Modified: mc/3D/CitcomS/trunk/lib/Lith_age.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Lith_age.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/Lith_age.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -1,6 +1,6 @@
 /*
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- * 
+ *
  *<LicenseText>
  *
  * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,7 +22,7 @@
  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  *
  *</LicenseText>
- * 
+ *
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  */
 
@@ -36,7 +36,7 @@
 #include "lith_age.h"
 
 float find_age_in_MY();
-void lith_age_restart_conform_tbc(struct All_variables *E);
+void lith_age_update_tbc(struct All_variables *E);
 
 
 void lith_age_input(struct All_variables *E)
@@ -104,7 +104,7 @@
 }
 
 
-void lith_age_restart_tic(struct All_variables *E)
+void lith_age_construct_tic(struct All_variables *E)
 {
   int i, j, k, m, node, nodeg;
   int nox, noy, noz, gnox, gnoy, gnoz;
@@ -134,15 +134,16 @@
 	    }
 	}
 
-  /* modify temperature BC to be concorded with restarted T */
-  lith_age_restart_conform_tbc(E);
+  /* modify temperature BC to be concorded with read in T */
+  lith_age_update_tbc(E);
+
   temperatures_conform_bcs(E);
 
   return;
 }
 
 
-void lith_age_restart_conform_tbc(struct All_variables *E)
+void lith_age_update_tbc(struct All_variables *E)
 {
   int i, j, k, m, node;
   int nox, noy, noz;
@@ -173,14 +174,6 @@
 }
 
 
-
-void lith_age_construct_tic(struct All_variables *E)
-{
-  lith_age_restart_tic(E);
-  return;
-}
-
-
 void lith_age_temperature_bound_adj(struct All_variables *E, int lv)
 {
   int j,node,nno;

Modified: mc/3D/CitcomS/trunk/lib/Makefile.am
===================================================================
--- mc/3D/CitcomS/trunk/lib/Makefile.am	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/Makefile.am	2007-08-16 21:36:18 UTC (rev 7837)
@@ -60,6 +60,8 @@
 	advection_diffusion.h \
 	advection.h \
 	BC_util.c \
+	Checkpoints.c \
+	checkpoints.h \
 	Citcom_init.c \
 	citcom_init.h \
 	Composition_related.c \

Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/Output.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -552,4 +552,3 @@
 
   return;
 }
-

Modified: mc/3D/CitcomS/trunk/lib/Problem_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Problem_related.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/Problem_related.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -1,6 +1,6 @@
 /*
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- * 
+ *
  *<LicenseText>
  *
  * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,7 +22,7 @@
  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  *
  *</LicenseText>
- * 
+ *
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  */
 #include <math.h>
@@ -70,8 +70,9 @@
     char output_file[255],input_s[1000];
 
     E->monitor.elapsed_time = 0.0;
-    if ((E->control.restart || E->control.post_p))    {
 
+    if (E->convection.tic_method == -1) {
+
 #ifdef USE_GZDIR		/* gzdir output */
       if(strcmp(E->output.format, "ascii-gz") == 0){
 	if(E->output.gzdir.vtk_io)
@@ -100,8 +101,8 @@
       if(rezip)
 	gzip_file(output_file);
 #endif
-    } /* end control.restart */
-    
+    } /* end if tic_method */
+
     return;
 }
 

Modified: mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -56,7 +56,6 @@
 
 void regional_tracer_setup(struct All_variables *E)
 {
-    void initialize_tracers();
 
     char output_file[255];
 
@@ -138,12 +137,7 @@
 
     make_mesh_ijk(E);
 
-
-    initialize_tracers(E);
-
-
     composition_setup(E);
-    tracer_post_processing(E);
 
     return;
 }
@@ -163,7 +157,7 @@
         fprintf(E->trace.fpt,"Reading tracer file %s\n",E->trace.tracer_file);
     }
     if (E->trace.ic_method==2) {
-        fprintf(E->trace.fpt,"Restarting Tracers\n");
+        fprintf(E->trace.fpt,"Read individual tracer files\n");
     }
 
     fprintf(E->trace.fpt,"Number of tracer flavors: %d\n", E->trace.nflavors);

Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -58,6 +58,8 @@
 void expand_later_array(struct All_variables *E, int j);
 void expand_tracer_arrays(struct All_variables *E, int j);
 void tracer_post_processing(struct All_variables *E);
+void allocate_tracer_arrays(struct All_variables *E,
+                            int j, int number_of_tracers);
 
 
 static void find_tracers(struct All_variables *E);
@@ -68,12 +70,10 @@
 static void generate_random_tracers(struct All_variables *E,
                                     int tracers_cap, int j);
 static void read_tracer_file(struct All_variables *E);
-static void restart_tracers(struct All_variables *E);
+static void read_old_tracer_file(struct All_variables *E);
 static void check_sum(struct All_variables *E);
 static int isum_tracers(struct All_variables *E);
 static void init_tracer_flavors(struct All_variables *E);
-static void allocate_tracer_arrays(struct All_variables *E,
-                                   int j, int number_of_tracers);
 static void reduce_tracer_arrays(struct All_variables *E);
 static void put_away_later(struct All_variables *E, int j, int it);
 static void eject_tracer(struct All_variables *E, int j, int it);
@@ -144,10 +144,10 @@
 	  parallel_process_termination();
 	  break;
 	}
-   
 
 
 
+
         if(E->parallel.nprocxy == 12)
             full_tracer_input(E);
 
@@ -629,7 +629,7 @@
     else if (E->trace.ic_method==1)
         read_tracer_file(E);
     else if (E->trace.ic_method==2)
-        restart_tracers(E);
+        read_old_tracer_file(E);
     else {
         fprintf(E->trace.fpt,"Not ready for other inputs yet\n");
         fflush(E->trace.fpt);
@@ -866,7 +866,7 @@
             /* XXX: if E->trace.number_of_extra_quantities is greater than 1 */
             /* this part has to be changed... */
             else {
-                fprintf(E->trace.fpt,"ERROR(restart tracers)-huh?\n");
+                fprintf(E->trace.fpt,"ERROR(read tracer file)-huh?\n");
                 fflush(E->trace.fpt);
                 exit(10);
             }
@@ -931,12 +931,12 @@
 }
 
 
-/************** RESTART TRACERS ******************************************/
+/************** READ OLD TRACER FILE *************************************/
 /*                                                                       */
-/* This function restarts tracers written from previous calculation      */
+/* This function read tracers written from previous calculation          */
 /* and the tracers are read as seperate files for each processor domain. */
 
-static void restart_tracers(struct All_variables *E)
+static void read_old_tracer_file(struct All_variables *E)
 {
 
     char output_file[200];
@@ -956,7 +956,7 @@
     FILE *fp1;
 
     if (E->trace.number_of_extra_quantities>99) {
-        fprintf(E->trace.fpt,"ERROR(restart_tracers)-increase size of extra[]\n");
+        fprintf(E->trace.fpt,"ERROR(read_old_tracer_file)-increase size of extra[]\n");
         fflush(E->trace.fpt);
         parallel_process_termination();
     }
@@ -972,7 +972,7 @@
     }else{
       sprintf(output_file,"%s.tracer.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
       if ( (fp1=fopen(output_file,"r"))==NULL) {
-        fprintf(E->trace.fpt,"ERROR(restart tracers)-file not found %s\n",output_file);
+        fprintf(E->trace.fpt,"ERROR(read_old_tracer_file)-gziped file not found %s\n",output_file);
         fflush(E->trace.fpt);
         exit(10);
       }
@@ -980,13 +980,13 @@
 #else
     sprintf(output_file,"%s.tracer.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
     if ( (fp1=fopen(output_file,"r"))==NULL) {
-        fprintf(E->trace.fpt,"ERROR(restart tracers)-file not found %s\n",output_file);
+        fprintf(E->trace.fpt,"ERROR(read_old_tracer_file)-file not found %s\n",output_file);
         fflush(E->trace.fpt);
         exit(10);
     }
 #endif
 
-    fprintf(stderr,"Restarting Tracers from %s\n",output_file);
+    fprintf(stderr,"Read old tracers from %s\n",output_file);
     fflush(stderr);
 
 
@@ -997,7 +997,7 @@
 
         /* some error control */
         if (E->trace.number_of_extra_quantities+3 != ncolumns) {
-            fprintf(E->trace.fpt,"ERROR(restart tracers)-wrong # of columns\n");
+            fprintf(E->trace.fpt,"ERROR(read_old_tracer_file)-wrong # of columns\n");
             fflush(E->trace.fpt);
             exit(10);
         }
@@ -1018,7 +1018,7 @@
             /* XXX: if E->trace.number_of_extra_quantities is greater than 1 */
             /* this part has to be changed... */
             else {
-                fprintf(E->trace.fpt,"ERROR(restart tracers)-huh?\n");
+                fprintf(E->trace.fpt,"ERROR(read_old_tracer_file)-huh?\n");
                 fflush(E->trace.fpt);
                 exit(10);
             }
@@ -1163,11 +1163,11 @@
       /* any tracer above z_interface is of flavor 0    */
       /* any tracer below z_interface is of flavor 1    */
       for (j=1;j<=E->sphere.caps_per_proc;j++) {
-	
+
 	number_of_tracers = E->trace.ntracers[j];
 	for (kk=1;kk<=number_of_tracers;kk++) {
 	  rad = E->trace.basicq[j][2][kk];
-	  
+
 	  if (rad<=E->trace.z_interface) E->trace.extraq[j][0][kk]=1.0;
 	  if (rad>E->trace.z_interface) E->trace.extraq[j][0][kk]=0.0;
 	}
@@ -1186,7 +1186,7 @@
 
 
     default:
-      
+
       fprintf(stderr,"ic_method_for_flavors %i undefined\n",E->trace.ic_method_for_flavors);
       parallel_process_termination();
       break;
@@ -1318,8 +1318,8 @@
 /*                                                                            */
 /* This function allocates memories to tracer arrays.                         */
 
-static void allocate_tracer_arrays(struct All_variables *E,
-                                   int j, int number_of_tracers)
+void allocate_tracer_arrays(struct All_variables *E,
+                            int j, int number_of_tracers)
 {
 
     int kk;

Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -234,21 +234,9 @@
       fflush(E->fp_out);
     }
 
+    /* interpolate from gauss quadrature points to node points for output */
     visc_from_gint_to_nodes(E,evisc,visc,E->mesh.levmax);
 
-    /*visc_from_nodes_to_gint(E,visc,evisc,E->mesh.levmax);*/
-
-    /*
-        visc_to_node_interpolate(E,evisc,visc);
-     */
-
-    /*    for(m=1;m<=E->sphere.caps_per_proc;m++) {
-          for(i=1;i<=E->lmesh.nel;i++)
-          if (i%E->lmesh.elz==0) {
-          fprintf(E->fp_out,"%.4e %.4e %.4e %5d %2d\n",E->eco[m][i].centre[1],E->eco[m][i].centre[2],log10(evisc[m][(i-1)*vpts+1]),i,E->mat[m][i]);
-
-	  }
-          }  */
     return;
 }
 

Added: mc/3D/CitcomS/trunk/lib/checkpoints.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/checkpoints.h	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/checkpoints.h	2007-08-16 21:36:18 UTC (rev 7837)
@@ -0,0 +1,31 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *<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 output_checkpoint(struct All_variables *E);
+void read_checkpoint(struct All_variables *E);

Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h	2007-08-16 21:36:18 UTC (rev 7837)
@@ -487,7 +487,7 @@
 
     float Q0;
     float Q0ER;
-  
+
     int precondition;
     int keep_going;
     int v_steps_low;
@@ -511,6 +511,7 @@
     int total_iteration_cycles;
     int total_v_solver_calls;
 
+    int checkpoint_frequency;
     int record_every;
     int record_all_until;
 
@@ -622,7 +623,6 @@
     int on;
 
     double buoyancy_ratio;
-    int ireset_initial_composition;
     int ibuoy_type;
 
     double *comp_el[13];

Modified: mc/3D/CitcomS/trunk/lib/initial_temperature.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/initial_temperature.h	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/initial_temperature.h	2007-08-16 21:36:18 UTC (rev 7837)
@@ -1,6 +1,6 @@
 /*
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- * 
+ *
  *<LicenseText>
  *
  * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,10 +22,9 @@
  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  *
  *</LicenseText>
- * 
+ *
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  */
 void convection_initial_temperature(struct All_variables *);
 void tic_input(struct All_variables *);
-void restart_tic(struct All_variables *);
 void construct_tic(struct All_variables *);

Modified: mc/3D/CitcomS/trunk/lib/lith_age.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/lith_age.h	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/lib/lith_age.h	2007-08-16 21:36:18 UTC (rev 7837)
@@ -1,6 +1,6 @@
 /*
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- * 
+ *
  *<LicenseText>
  *
  * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,12 +22,11 @@
  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  *
  *</LicenseText>
- * 
+ *
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  */
 void lith_age_input(struct All_variables *E);
 void lith_age_init(struct All_variables *E);
-void lith_age_restart_tic(struct All_variables *E);
 void lith_age_construct_tic(struct All_variables *E) ;
 void lith_age_temperature_bound_adj(struct All_variables *E, int lv);
 void lith_age_conform_tbc(struct All_variables *E);

Modified: mc/3D/CitcomS/trunk/module/bindings.c
===================================================================
--- mc/3D/CitcomS/trunk/module/bindings.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/module/bindings.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -169,16 +169,21 @@
 
     /* from initial_conditions.h */
 
+    {pyCitcom_ic_initialize_material__name__,
+     pyCitcom_ic_initialize_material,
+     METH_VARARGS,
+     pyCitcom_ic_initialize_material__doc__},
+
+    {pyCitcom_ic_init_tracer_composition__name__,
+     pyCitcom_ic_init_tracer_composition,
+     METH_VARARGS,
+     pyCitcom_ic_init_tracer_composition__doc__},
+
     {pyCitcom_ic_constructTemperature__name__,
      pyCitcom_ic_constructTemperature,
      METH_VARARGS,
      pyCitcom_ic_constructTemperature__doc__},
 
-    {pyCitcom_ic_restartTemperature__name__,
-     pyCitcom_ic_restartTemperature,
-     METH_VARARGS,
-     pyCitcom_ic_restartTemperature__doc__},
-
     {pyCitcom_ic_initPressure__name__,
      pyCitcom_ic_initPressure,
      METH_VARARGS,
@@ -194,6 +199,11 @@
      METH_VARARGS,
      pyCitcom_ic_initViscosity__doc__},
 
+    {pyCitcom_ic_readCheckpoint__name__,
+     pyCitcom_ic_readCheckpoint,
+     METH_VARARGS,
+     pyCitcom_ic_readCheckpoint__doc__},
+
     /* from mesher.h */
 
     {pyCitcom_full_sphere_launch__name__,
@@ -223,6 +233,11 @@
      METH_VARARGS,
      pyCitcom_output_time__doc__},
 
+    {pyCitcom_output_checkpoint__name__,
+     pyCitcom_output_checkpoint,
+     METH_VARARGS,
+     pyCitcom_output_checkpoint__doc__},
+
     /* from setProperties.h */
 
     {pyCitcom_Advection_diffusion_set_properties__name__,

Modified: mc/3D/CitcomS/trunk/module/initial_conditions.c
===================================================================
--- mc/3D/CitcomS/trunk/module/initial_conditions.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/module/initial_conditions.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -32,62 +32,84 @@
 #include "global_defs.h"
 
 
-void construct_tic(struct All_variables*);
-void debug_tic(struct All_variables*);
+void initialize_material(struct All_variables*);
+void initialize_tracers(struct All_variables*);
+void init_composition(struct All_variables*);
 void initial_pressure(struct All_variables*);
 void initial_velocity(struct All_variables*);
 void initial_viscosity(struct All_variables*);
 void report(struct All_variables*, char* str);
-void restart_tic(struct All_variables*);
+void read_checkpoint(struct All_variables*);
 
 
-char pyCitcom_ic_constructTemperature__doc__[] = "";
-char pyCitcom_ic_constructTemperature__name__[] = "constructTemperature";
+char pyCitcom_ic_initialize_material__doc__[] = "";
+char pyCitcom_ic_initialize_material__name__[] = "initialize_material";
 
-PyObject * pyCitcom_ic_constructTemperature(PyObject *self, PyObject *args)
+PyObject * pyCitcom_ic_initialize_material(PyObject *self, PyObject *args)
 {
     PyObject *obj;
     struct All_variables* E;
 
-    if (!PyArg_ParseTuple(args, "O:constructTemperature", &obj))
+    if (!PyArg_ParseTuple(args, "O:initialize_material", &obj))
         return NULL;
 
     E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
 
-    report(E,"Initialize temperature field");
-    construct_tic(E);
+    initialize_material(E);
 
-    if(E->control.verbose)
-        debug_tic(E);
-    
     Py_INCREF(Py_None);
     return Py_None;
 }
 
 
 
-char pyCitcom_ic_restartTemperature__doc__[] = "";
-char pyCitcom_ic_restartTemperature__name__[] = "restartTemperature";
+char pyCitcom_ic_init_tracer_composition__doc__[] = "";
+char pyCitcom_ic_init_tracer_composition__name__[] = "init_tracer_composition";
 
-PyObject * pyCitcom_ic_restartTemperature(PyObject *self, PyObject *args)
+PyObject * pyCitcom_ic_init_tracer_composition(PyObject *self, PyObject *args)
 {
     PyObject *obj;
     struct All_variables* E;
 
-    if (!PyArg_ParseTuple(args, "O:restartTemperature", &obj))
+    if (!PyArg_ParseTuple(args, "O:init_tracer_composition", &obj))
         return NULL;
 
     E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
 
-    report(E,"Initialize temperature field");
-    restart_tic(E);
+    if (E->control.tracer==1) {
+        initialize_tracers(E);
 
+        if (E->composition.on)
+            init_composition(E);
+    }
+
     Py_INCREF(Py_None);
     return Py_None;
 }
 
 
 
+char pyCitcom_ic_constructTemperature__doc__[] = "";
+char pyCitcom_ic_constructTemperature__name__[] = "constructTemperature";
+
+PyObject * pyCitcom_ic_constructTemperature(PyObject *self, PyObject *args)
+{
+    PyObject *obj;
+    struct All_variables* E;
+
+    if (!PyArg_ParseTuple(args, "O:constructTemperature", &obj))
+        return NULL;
+
+    E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
+
+    (E->problem_initial_fields)(E);
+
+    Py_INCREF(Py_None);
+    return Py_None;
+}
+
+
+
 char pyCitcom_ic_initPressure__doc__[] = "";
 char pyCitcom_ic_initPressure__name__[] = "initPressure";
 
@@ -154,6 +176,27 @@
 
 
 
+char pyCitcom_ic_readCheckpoint__doc__[] = "";
+char pyCitcom_ic_readCheckpoint__name__[] = "readCheckpoint";
+
+PyObject * pyCitcom_ic_readCheckpoint(PyObject *self, PyObject *args)
+{
+    PyObject *obj;
+    struct All_variables* E;
+
+    if (!PyArg_ParseTuple(args, "O:readCheckpoint", &obj))
+        return NULL;
+
+    E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
+
+    read_checkpoint(E);
+
+    Py_INCREF(Py_None);
+    return Py_None;
+}
+
+
+
 /* $Id$ */
 
 /* End of file */

Modified: mc/3D/CitcomS/trunk/module/initial_conditions.h
===================================================================
--- mc/3D/CitcomS/trunk/module/initial_conditions.h	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/module/initial_conditions.h	2007-08-16 21:36:18 UTC (rev 7837)
@@ -29,16 +29,21 @@
 #define pyCitcom_initial_conditions_h
 
 
+extern char pyCitcom_ic_initialize_material__name__[];
+extern char pyCitcom_ic_initialize_material__doc__[];
+PyObject * pyCitcom_ic_initialize_material(PyObject *, PyObject *);
+
+
+extern char pyCitcom_ic_init_tracer_composition__name__[];
+extern char pyCitcom_ic_init_tracer_composition__doc__[];
+PyObject * pyCitcom_ic_init_tracer_composition(PyObject *, PyObject *);
+
+
 extern char pyCitcom_ic_constructTemperature__name__[];
 extern char pyCitcom_ic_constructTemperature__doc__[];
 PyObject * pyCitcom_ic_constructTemperature(PyObject *, PyObject *);
 
 
-extern char pyCitcom_ic_restartTemperature__name__[];
-extern char pyCitcom_ic_restartTemperature__doc__[];
-PyObject * pyCitcom_ic_restartTemperature(PyObject *, PyObject *);
-
-
 extern char pyCitcom_ic_initPressure__name__[];
 extern char pyCitcom_ic_initPressure__doc__[];
 PyObject * pyCitcom_ic_initPressure(PyObject *, PyObject *);
@@ -53,6 +58,12 @@
 extern char pyCitcom_ic_initViscosity__doc__[];
 PyObject * pyCitcom_ic_initViscosity(PyObject *, PyObject *);
 
+
+extern char pyCitcom_ic_readCheckpoint__name__[];
+extern char pyCitcom_ic_readCheckpoint__doc__[];
+PyObject * pyCitcom_ic_readCheckpoint(PyObject *, PyObject *);
+
+
 #endif
 
 /* $Id$ */

Modified: mc/3D/CitcomS/trunk/module/outputs.c
===================================================================
--- mc/3D/CitcomS/trunk/module/outputs.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/module/outputs.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -102,6 +102,28 @@
 }
 
 
+char pyCitcom_output_checkpoint__doc__[] = "";
+char pyCitcom_output_checkpoint__name__[] = "output_checkpoint";
+
+PyObject * pyCitcom_output_checkpoint(PyObject *self, PyObject *args)
+{
+    void read_checkpoint(struct All_variables*);
+    PyObject *obj;
+    struct All_variables* E;
+
+    if (!PyArg_ParseTuple(args, "O:output_checkpoint", &obj))
+        return NULL;
+
+    E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
+
+    output_checkpoint(E);
+
+
+    Py_INCREF(Py_None);
+    return Py_None;
+}
+
+
 /* $Id$ */
 
 /* End of file */

Modified: mc/3D/CitcomS/trunk/module/outputs.h
===================================================================
--- mc/3D/CitcomS/trunk/module/outputs.h	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/module/outputs.h	2007-08-16 21:36:18 UTC (rev 7837)
@@ -40,6 +40,10 @@
 extern char pyCitcom_output_time__doc__[];
 PyObject * pyCitcom_output_time(PyObject *, PyObject *);
 
+extern char pyCitcom_output_checkpoint__name__[];
+extern char pyCitcom_output_checkpoint__doc__[];
+PyObject * pyCitcom_output_checkpoint(PyObject *, PyObject *);
+
 #endif
 
 /* $Id$ */

Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/module/setProperties.c	2007-08-16 21:36:18 UTC (rev 7837)
@@ -465,6 +465,7 @@
     E->advection.min_timesteps = 1;
     E->advection.max_timesteps = 1;
     E->advection.max_total_timesteps = 1;
+    E->control.checkpoint_frequency = 1;
     E->control.record_every = 1;
     E->control.record_all_until = 1;
 
@@ -612,7 +613,6 @@
     if (E->composition.ichemical_buoyancy==1) {
         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);
     }
 
 

Added: mc/3D/CitcomS/trunk/tests/checkpoint/restart.cfg
===================================================================
--- mc/3D/CitcomS/trunk/tests/checkpoint/restart.cfg	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/tests/checkpoint/restart.cfg	2007-08-16 21:36:18 UTC (rev 7837)
@@ -0,0 +1,13 @@
+
+#[CitcomS]
+#steps = 4
+
+[CitcomS.solver]
+datafile = yyy
+datadir_old = .
+datafile_old = zzz
+
+[CitcomS.solver.ic]
+restart = on
+solution_cycles_init = 2
+

Added: mc/3D/CitcomS/trunk/tests/checkpoint/run.sh
===================================================================
--- mc/3D/CitcomS/trunk/tests/checkpoint/run.sh	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/tests/checkpoint/run.sh	2007-08-16 21:36:18 UTC (rev 7837)
@@ -0,0 +1,8 @@
+#!/bin/sh
+
+## Run an ordinary job then restart it. The end result should be the same
+
+../../bin/citcoms tracer.cfg && ../../bin/citcoms tracer.cfg restart.cfg && diff *.tracer.0.5 && diff *.velo.0.5
+
+## clean up
+rm yyy.* zzz.* pid*.cfg

Added: mc/3D/CitcomS/trunk/tests/checkpoint/tracer.cfg
===================================================================
--- mc/3D/CitcomS/trunk/tests/checkpoint/tracer.cfg	2007-08-16 20:10:58 UTC (rev 7836)
+++ mc/3D/CitcomS/trunk/tests/checkpoint/tracer.cfg	2007-08-16 21:36:18 UTC (rev 7837)
@@ -0,0 +1,33 @@
+
+[CitcomS]
+steps = 5
+
+[CitcomS.controller]
+monitoringFrequency = 1
+checkpointFrequency = 1
+
+[CitcomS.solver]
+datafile = zzz
+
+[CitcomS.solver.mesher]
+nprocx =  1
+nprocy =  1
+nodex  =  5
+nodey  =  5
+nodez  =  5
+
+[CitcomS.solver.tracer]
+tracer = on
+tracer_ic_method = 0
+tracers_per_element = 20
+
+tracer_flavors = 2
+ic_method_for_flavors = 0
+z_interface = 0.7
+
+chemical_buoyancy = 1
+buoyancy_ratio = 0.4
+
+regular_grid_deltheta = 1.0
+regular_grid_delphi = 1.0
+



More information about the cig-commits mailing list