diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index e81571e0..32c7f4be 100644 --- a/python/ctem/ctem/driver.py +++ b/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() diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py index d0d383e5..7d363f2f 100644 --- a/python/ctem/ctem/util.py +++ b/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 diff --git a/python/ctem/ctem/viewer3d.py b/python/ctem/ctem/viewer3d.py index e3d66cb7..b7c3267d 100644 --- a/python/ctem/ctem/viewer3d.py +++ b/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 diff --git a/python/ctem/tests/run_and_reprocess/ctem_cleanup.py b/python/ctem/tests/run_and_reprocess/ctem_cleanup.py new file mode 100644 index 00000000..d9de5e24 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/python/ctem/tests/run_and_reprocess/ctem_compute_one_interval.py b/python/ctem/tests/run_and_reprocess/ctem_compute_one_interval.py new file mode 100644 index 00000000..fe9ecb39 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/python/ctem/tests/run_and_reprocess/ctem_driver.py b/python/ctem/tests/run_and_reprocess/ctem_driver.py index 50392285..04b5970d 100644 --- a/python/ctem/tests/run_and_reprocess/ctem_driver.py +++ b/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() \ No newline at end of file diff --git a/python/ctem/tests/run_and_reprocess/ctem_postprocess.py b/python/ctem/tests/run_and_reprocess/ctem_postprocess.py new file mode 100644 index 00000000..5d56dac5 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/python/ctem/tests/run_and_reprocess/ctem_preprocess.py b/python/ctem/tests/run_and_reprocess/ctem_preprocess.py new file mode 100644 index 00000000..4f982486 --- /dev/null +++ b/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() \ No newline at end of file