Skip to content

Commit

Permalink
Merge branch 'realistic'
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Feb 24, 2022
2 parents d7630cc + 0295dbb commit 3f4b32e
Show file tree
Hide file tree
Showing 17 changed files with 2,117 additions and 106 deletions.
2 changes: 2 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,3 +69,5 @@ examples/global-lunar-bombardment/scale.ipynb
*.code-workspace

*.png

python/ctem/ctem.egg-info/
4 changes: 2 additions & 2 deletions examples/global-lunar-bombardment/ctem.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,9 +34,9 @@ runtype single ! Run type: options are normal /

! CTEM required inputs
seed 76535 ! Random number generator seed
gridsize 2000 ! Size of grid in pixels
gridsize 200 ! Size of grid in pixels
numlayers 10 ! Number of perched layers
pix 3.08e3 ! Pixel size (m)
pix 3.08e4 ! Pixel size (m)
mat rock ! Material (rock or ice)
! Bedrock scaling parameters
mu_b 0.55e0 ! Experimentally derived parameter for bedrock crater scaling law
Expand Down
29 changes: 0 additions & 29 deletions python/ctem/ctem/dist_reader.py
View file

This file was deleted.

190 changes: 118 additions & 72 deletions python/ctem/ctem/driver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,8 @@
import os
import subprocess
import shutil
from ctem import io
from ctem import util
import sys

class Simulation:
"""
Expand Down Expand Up @@ -42,8 +43,10 @@ def __init__(self, param_file="ctem.in", isnew=True):
'sfdcompare': None,
'sfdfile': None
}
self.user = io.read_user_input(self.user)

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',
Expand All @@ -59,8 +62,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)

Expand All @@ -71,37 +78,33 @@ def __init__(self, param_file="ctem.in", isnew=True):
self.pdist = np.zeros([1, 6])
self.tdist = np.zeros([1, 6])
self.surface_dem = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float)
self.regolith = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float)
self.surface_ejc = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float)

if self.user['sfdcompare'] is not None:
# Read sfdcompare file
sfdfile = os.path.join(self.user['workingdir'], self.user['sfdcompare'])
self.ph1 = io.read_formatted_ascii(sfdfile, skip_lines=0)
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')

io.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)
# 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')

# Scale the production function to the simulation domain
self.scale_production()

if (self.user['runtype'].upper() == 'STATISTICAL'):
self.user['ncount'] = 1
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)

# Write ctem.dat file
io.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
# Scale the production function to the simulation domain
self.scale_production()

util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
else:
self.user['ncount'] = 0
else:
print('Continuing a previous run')
self.process_interval(isnew=False)
print('Continuing a previous run')
self.read_output()

return

Expand All @@ -113,15 +116,15 @@ def scale_production(self):

# Read production function file
impfile = os.path.join(self.user['workingdir'], self.user['impfile'])
prodfunction = io.read_formatted_ascii(impfile, skip_lines=0)
prodfunction = util.read_formatted_ascii(impfile, skip_lines=0)

# Create impactor production population
area = (self.user['gridsize'] * self.user['pix']) ** 2
self.production = np.copy(prodfunction)
self.production[:, 1] = self.production[:, 1] * area * self.user['interval']

# Write the scaled production function to file
io.write_production(self.user, self.production)
util.write_production(self.user, self.production)

return

Expand All @@ -142,40 +145,31 @@ def run(self):
#Execute Fortran code
self.compute_one_interval()

# Process output files
self.process_interval()
# 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
io.image_dem(self.user, self.surface_dem)
if (self.user['saverego'].upper() == 'T'):
io.image_regolith(self.user, self.surface_ejc)
if (self.user['saveshaded'].upper() == 'T'):
io.image_shaded_relief(self.user, self.surface_dem)
if (self.user['savepres'].upper() == 'T'):
io.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
self.user['masstot'] = 0.0

# Delete tdistribution file, if it exists
# Delete tdistribution file, if it exists so that we don't accumulate true craters in statistical mode
tdist_file = self.user['workingdir'] + 'tdistribution.dat'
if os.path.isfile(tdist_file):
os.remove(tdist_file)

else:
self.user['restart'] = 'T'

# Write ctem.dat file
io.write_datfile(self.user, self.output_filenames['dat'],self.seedarr)
# Write ctem.dat file for next interval
util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)

return

Expand All @@ -185,15 +179,25 @@ def compute_one_interval(self):
"""
# Create crater population and display CTEM progress on screen
print(self.user['ncount'], ' Calling FORTRAN routine')
with subprocess.Popen([os.path.join(self.user['workingdir'], 'CTEM')],
try:
p = subprocess.Popen([os.path.join(self.user['workingdir'], 'CTEM')],
stdout=subprocess.PIPE,
universal_newlines=True) as p:
stderr=subprocess.PIPE,
universal_newlines=True)
for line in p.stdout:
print(line, end='')
res = p.communicate()
if p.returncode != 0:
for line in res[1]:
print(line, end='')
raise Exception ("CTEM Failure")
except:
print("Error executing main CTEM program")
sys.exit()

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.
Expand All @@ -202,41 +206,60 @@ def process_interval(self, isnew=True):
print(self.user['ncount'], ' Reading Fortran output')

# Read surface dem(shaded relief) and ejecta data files
self.surface_dem = io.read_unformatted_binary(self.output_filenames['dem'], self.user['gridsize'])
self.surface_ejc = io.read_unformatted_binary(self.output_filenames['ejc'], self.user['gridsize'])
self.surface_dem = util.read_unformatted_binary(self.output_filenames['dem'], self.user['gridsize'])
self.surface_ejc = util.read_unformatted_binary(self.output_filenames['ejc'], self.user['gridsize'])

# Read odistribution, tdistribution, and pdistribution files
self.odist = io.read_formatted_ascii(self.output_filenames['odist'], skip_lines=1)
self.tdist = io.read_formatted_ascii(self.output_filenames['tdist'], skip_lines=1)
self.pdist = io.read_formatted_ascii(self.output_filenames['pdist'], skip_lines=1)
self.odist = util.read_formatted_ascii(self.output_filenames['odist'], skip_lines=1)
self.tdist = util.read_formatted_ascii(self.output_filenames['tdist'], skip_lines=1)
self.pdist = util.read_formatted_ascii(self.output_filenames['pdist'], skip_lines=1)

# Read impact mass from file
impact_mass = io.read_impact_mass(self.output_filenames['impmass'])
self.impact_mass = util.read_impact_mass(self.output_filenames['impmass'])

# Read ctem.dat file
io.read_datfile(self.user, self.output_filenames['dat'], self.seedarr)

# 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)
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')

if isnew:
# 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')
return



def redirect_outputs(self, filekeys, foldername):
"""
Expand All @@ -256,6 +279,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()
Expand Down
Loading

0 comments on commit 3f4b32e

Please sign in to comment.