Improved robustness of the Python driver. Separated pre and post-processing steps more thoroughly so that dedicated scripts can be written to do those tasks independently

python/ctem/ctem/driver.py

@@ -43,7 +43,9 @@ def __init__(self, param_file="ctem.in", isnew=True):
             'sfdfile': None
         }
         self.user = util.read_user_input(self.user)
-
+
+        # This is containsall files generated by the main Fortran CTEM program plus the ctem.dat file that
+        # communicates the state of the simulation from the Python driver to the Fortran program
         self.output_filenames = {
             'dat': 'ctem.dat',
             'dem': 'surface_dem.dat',
@@ -59,8 +61,12 @@ def __init__(self, param_file="ctem.in", isnew=True):
             'impmass' : 'impactmass.dat',
             'fracdone' : 'fracdone.dat',
             'regodepth' : 'regolithdepth.dat',
+            'ejmax' : 'ejecta_table_max.dat',
+            'ejmin' : 'ejecta_table_min.dat',
+            'testprof' : 'testprofile.dat',
+            'craterscale' : 'craterscale.dat'
         }
-        
+
         for k, v in self.output_filenames.items():
             self.output_filenames[k] = os.path.join(currentdir, v)
 
@@ -79,23 +85,25 @@ def __init__(self, param_file="ctem.in", isnew=True):
             self.ph1 = util.read_formatted_ascii(sfdfile, skip_lines=0)
 
 
-        # Starting new or old run?
-        if (self.user['restart'].upper() == 'F' and isnew):
-            print('Starting a new run')
+        # If this is a new simulation, run and post-process results. Otherwise, don't do anything more
+        if isnew:
+            # Starting new or old run?
+            if (self.user['restart'].upper() == 'F'):
+                print('Starting a new run')
 
-            util.create_dir_structure(self.user)
-            # Delete any old output files
-            for k, v in self.output_filenames.items():
-                if os.path.isfile(v):
-                    os.remove(v)
+                util.create_dir_structure(self.user)
+                # Delete any old output files
+                for k, v in self.output_filenames.items():
+                    if os.path.isfile(v):
+                        os.remove(v)
 
-            # Scale the production function to the simulation domain
-            self.scale_production()
-
-            util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
-        else:
-            print('Continuing a previous run')
-            self.process_interval(isnew=False)
+                # Scale the production function to the simulation domain
+                self.scale_production()
+                
+                util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
+            else:
+                print('Continuing a previous run')
+                self.read_output()
 
         return
 
@@ -136,31 +144,16 @@ def run(self):
                 #Execute Fortran code
                 self.compute_one_interval()
 
-                # Process output files
-                self.process_interval()
-
-                # Redirect output files to storage
-                self.redirect_outputs(['odist', 'ocum', 'pdist', 'tdist'], 'dist')
-                if (self.user['savetruelist'].upper() == 'T'):
-                    self.redirect_outputs(['tcum'], 'dist')
-                self.redirect_outputs(['impmass'], 'misc')
+                # Read in output files
+                self.read_output()
 
-            # Display results
-            print(self.user['ncount'], '  Displaying results')
-
-            # Write surface dem, surface ejecta, shaded relief, and rplot data
-            util.image_dem(self.user, self.surface_dem)
-            if (self.user['saverego'].upper() == 'T'):
-                util.image_regolith(self.user, self.surface_ejc)
-            if (self.user['saveshaded'].upper() == 'T'):
-                util.image_shaded_relief(self.user, self.surface_dem)
-            if (self.user['savepres'].upper() == 'T'):
-                util.create_rplot(self.user, self.odist, self.pdist, self.tdist, self.ph1)
+            # Process the o utput files
+            self.process_output()
 
             # Update ncount
             self.user['ncount'] = self.user['ncount'] + 1
 
-            if ((self.user['runtype'].upper() == 'STATISTICAL') or (self.user['ncount'] == 1)):
+            if ((self.user['runtype'].upper() == 'STATISTICAL') or (self.user['ncount'] == 1)): # Reset the simulation
                 self.user['restart'] = 'F'
                 self.user['curyear'] = 0.0
                 self.user['totalimpacts'] = 0
@@ -174,7 +167,7 @@ def run(self):
             else:
                 self.user['restart'] = 'T'
 
-            # Write ctem.dat file
+            # Write ctem.dat file for next interval
             util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
 
         return
@@ -196,7 +189,7 @@ def compute_one_interval(self):
 
         return
 
-    def process_interval(self, isnew=True):
+    def read_output(self):
         """
         Reads in all of the files generated by the Fortran code after one interval and redirects them to appropriate
         folders for storage after appending the interval number to the filename.
@@ -214,25 +207,50 @@ def process_interval(self, isnew=True):
         self.pdist = util.read_formatted_ascii(self.output_filenames['pdist'], skip_lines=1)
 
         # Read impact mass from file
-        impact_mass = util.read_impact_mass(self.output_filenames['impmass'])
+        self.impact_mass = util.read_impact_mass(self.output_filenames['impmass'])
 
         # Read ctem.dat file
         util.read_datfile(self.user, self.output_filenames['dat'], self.seedarr)
+
+    def process_output(self):
+        """
+        Processes output to update masses, time, computes regolith depth, generates all plots, and redirects files
+        into intermediate storage.
+        """
+        # Display results
+        print(self.user['ncount'], '  Generating surface images and plots')
 
-        # Save copy of crater distribution files
-        # Update user: mass, curyear, regolith properties
-        self.user['masstot'] = self.user['masstot'] + impact_mass
-
-        self.user['curyear'] = self.user['curyear'] + self.user['fracdone'] * self.user['interval']
-        template = "%(fracdone)9.6f %(curyear)19.12E\n"
-        with open(self.output_filenames['fracdone'], 'w') as fp_frac:
-            fp_frac.write(template % self.user)
-
-        reg_text = "%19.12E %19.12E %19.12E %19.12E\n" % (self.user['curyear'],
-                                                          np.mean(self.surface_ejc), np.amax(self.surface_ejc),
-                                                          np.amin(self.surface_ejc))
-        with open(self.output_filenames['regodepth'], 'w') as fp_reg:
-            fp_reg.write(reg_text)
+        # Write surface dem, surface ejecta, shaded relief, and rplot data
+        util.image_dem(self.user, self.surface_dem)
+        if (self.user['saverego'].upper() == 'T'):
+            util.image_regolith(self.user, self.surface_ejc)
+        if (self.user['saveshaded'].upper() == 'T'):
+            util.image_shaded_relief(self.user, self.surface_dem)
+
+        if self.user['ncount'] > 0: # These aren't available yet from the initial conditions
+            if (self.user['savepres'].upper() == 'T'):
+                util.create_rplot(self.user, self.odist, self.pdist, self.tdist, self.ph1)
+
+            # Save copy of crater distribution files
+            # Update user: mass, curyear, regolith properties
+            self.user['masstot'] = self.user['masstot'] + self.impact_mass
+
+            self.user['curyear'] = self.user['curyear'] + self.user['fracdone'] * self.user['interval']
+            template = "%(fracdone)9.6f %(curyear)19.12E\n"
+            with open(self.output_filenames['fracdone'], 'w') as fp_frac:
+                fp_frac.write(template % self.user)
+
+            reg_text = "%19.12E %19.12E %19.12E %19.12E\n" % (self.user['curyear'],
+                                                              np.mean(self.surface_ejc), np.amax(self.surface_ejc),
+                                                              np.amin(self.surface_ejc))
+            with open(self.output_filenames['regodepth'], 'w') as fp_reg:
+                fp_reg.write(reg_text)
+            # Redirect output files to storage
+            self.redirect_outputs(['odist', 'ocum', 'pdist', 'tdist'], 'dist')
+            if (self.user['savetruelist'].upper() == 'T'):
+                self.redirect_outputs(['tcum'], 'dist')
+            self.redirect_outputs(['impmass'], 'misc')
+
 
 
     def redirect_outputs(self, filekeys, foldername):
@@ -253,6 +271,29 @@ def redirect_outputs(self, filekeys, foldername):
             forig = self.output_filenames[k]
             fdest = os.path.join(self.user['workingdir'], foldername, f"{k}_{self.user['ncount']:06d}.dat")
             shutil.copy2(forig, fdest)
+
+    def cleanup(self):
+        """
+        Deletes all files and folders generated by CTEM.
+        """
+        # This is a list of files generated by the main Fortran program
+        print("Deleting all files generated by CTEM")
+        util.destroy_dir_structure(self.user)
+        for key, filename in self.output_filenames.items():
+            print(f"Deleting file {filename}")
+            try:
+                os.remove(filename)
+            except OSError as error:
+                print(error)
+        # Delete scaled production file
+        prodfile = self.user['sfdfile']
+        print(f"Deleting file {prodfile}")
+        try:
+            os.remove(prodfile)
+        except OSError as error:
+            print(error)
+
+        return
 
 if __name__ == '__main__':
     sim = Simulation()

python/ctem/ctem/util.py

@@ -8,17 +8,33 @@
 # Set pixel scaling common for image writing, at 1 pixel/ array element
 dpi = 72.0
 
+# These are directories that are created by CTEM in order to store intermediate results
+directories = ['dist', 'misc', 'rego', 'rplot', 'surf', 'shaded']
+
 def create_dir_structure(user):
-    # Create directories for various output files if they do not already exist
-    directories = ['dist', 'misc', 'rego', 'rplot', 'surf', 'shaded']
-
-    for directory in directories:
-        dir_test = os.path.join(user['workingdir'], directory)
-        if not os.path.isdir(dir_test):
-            os.makedirs(dir_test)
+    """
+    Create directories for various output files if they do not already exist
+    """
+    for dirname in directories:
+        dirpath = os.path.join(user['workingdir'], dirname)
+        if not os.path.isdir(dirpath):
+            print(f"Creating directory {dirpath}")
+            os.makedirs(dirpath)
 
     return
 
+def destroy_dir_structure(user):
+    """
+    Deletes directories generated by create_dir_structure
+    """
+    for dirname in directories:
+        dirpath = os.path.join(user['workingdir'], dirname)
+        if os.path.isdir(dirpath):
+            print(f"Deleting directory {dirpath}")
+            shutil.rmtree(dirpath, ignore_errors=True)
+
+    return
+
 
 def create_rplot(user, odist, pdist, tdist, ph1):
     # Parameters: empirical saturation limit and dfrac
@@ -239,7 +255,11 @@ def read_datfile(user, datfile, seedarr):
 def read_formatted_ascii(filename, skip_lines):
     # Generalized ascii text reader
     # For use with sfdcompare, production, odist, tdist, pdist data files
-    data = np.genfromtxt(filename, skip_header=skip_lines)
+    try:
+        data = np.genfromtxt(filename, skip_header=skip_lines)
+    except:
+        print(f"Error reading {filename}")
+        data = None
     return data
 
 

python/ctem/ctem/viewer3d.py

@@ -6,6 +6,7 @@ class Polysurface(ctem.Simulation):
     """A model of a self-gravitating small body"""
     def __init__(self):
         ctem.Simulation.__init__(self, isnew=False) # Initialize the data structures, but doesn't start a new run
+        self.read_output(isnew=False)
         #used for Open3d module
         self.mesh = open3d.geometry.TriangleMesh()
         return

python/ctem/tests/run_and_reprocess/ctem_cleanup.py

@@ -0,0 +1,6 @@
+"""
+Deletes all output generated by a previous simulation
+"""
+import ctem
+sim = ctem.Simulation(isnew=False)
+sim.cleanup()

python/ctem/tests/run_and_reprocess/ctem_compute_one_interval.py

@@ -0,0 +1,7 @@
+"""
+This is an example of a standalone script. Files must be generated by a preprocessing script (this is an atypical useage
+and is only here for testing purposes).
+"""
+import ctem
+sim = ctem.Simulation(isnew=False)
+sim.compute_one_interval()

python/ctem/tests/run_and_reprocess/ctem_driver.py

@@ -1,3 +1,7 @@
+"""
+This is an example of a standard "driver" script. This will read in the input files, pre-process them for the main
+Fortran code and post-process them.
+"""
 import ctem
 sim = ctem.Simulation()
 sim.run()

python/ctem/tests/run_and_reprocess/ctem_postprocess.py

@@ -0,0 +1,10 @@
+"""
+This is an example of a post-processing script. This will read in the input and output files and process them like as in
+a typical CTEM simulation
+
+"""
+import ctem
+sim = ctem.Simulation(isnew=False)
+sim.read_output()
+sim.user['ncount'] += 1
+sim.process_output()

python/ctem/tests/run_and_reprocess/ctem_preprocess.py

@@ -0,0 +1,5 @@
+"""
+This is an example of a pre-processing script. This will read in the input files, but won't execute a simulation.
+"""
+import ctem
+sim = ctem.Simulation()