Skip to content

Commit

Permalink
Pulled driver.py from the allocatable_arrays branch
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Mar 10, 2023
1 parent aa5b2d8 commit 52160b5
Showing 1 changed file with 102 additions and 38 deletions.
140 changes: 102 additions & 38 deletions python/ctem/ctem/driver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,12 @@
import shutil
from ctem import util
import sys
import pandas
from ctem import craterproduction
from scipy.interpolate import interp1d
from ctem import __file__ as _pyfile
from pathlib import Path
import warnings

class Simulation:
"""
Expand All @@ -19,7 +25,6 @@ def __init__(self, param_file="ctem.in", isnew=True):
'popupconsole': None,
'saveshaded': None,
'saverego': None,
'savepres': None,
'savetruelist': None,
'seedn': 1,
'totalimpacts': 0,
Expand All @@ -33,18 +38,31 @@ def __init__(self, param_file="ctem.in", isnew=True):
'gridsize': -1,
'seed': 0,
'maxcrat': 1.0,
'shadedminhdefault': 0,
'shadedmaxhdefault': 0,
'shadedminhdefault': 1,
'shadedmaxhdefault': 1,
'shadedminh': 0.0,
'shadedmaxh': 0.0,
'workingdir': currentdir,
'ctemfile': os.path.join(currentdir, param_file),
'ctemfile': param_file,
'impfile': None,
'sfdcompare': None,
'sfdfile': None
'quasimc': None,
'sfdfile' : None,
'realcraterlist': None,
}

# Get the location of the CTEM executable
self.ctem_executable = Path(_pyfile).parent.parent.parent.parent / "build" / "src" / "CTEM"
if not self.ctem_executable.exists():
print(f"CTEM driver not found at {self.ctem_executable}. Trying current directory.")
self.ctem_executable = Path(currentdir) / "CTEM"
if not self.ctem_executable.exists():
warnings.warn(f"Cannot find the CTEM driver {str(self.ctem_executable)}", stacklevel=2)
self.ctem_executable = 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 = {
Expand All @@ -54,6 +72,11 @@ def __init__(self, param_file="ctem.in", isnew=True):
'ejc': 'surface_ejc.dat',
'pos': 'surface_pos.dat',
'time': 'surface_time.dat',
'stack': 'surface_stacknum.dat',
'rego': 'surface_rego.dat',
'melt': 'surface_melt.dat',
'comp': 'surface_comp.dat',
'age' : 'surface_age.dat',
'ocum': 'ocumulative.dat',
'odist': 'odistribution.dat',
'pdist': 'pdistribution.dat',
Expand All @@ -65,12 +88,22 @@ def __init__(self, param_file="ctem.in", isnew=True):
'ejmax' : 'ejecta_table_max.dat',
'ejmin' : 'ejecta_table_min.dat',
'testprof' : 'testprofile.dat',
'craterscale' : 'craterscale.dat'
'craterscale' : 'craterscale.dat',
'craterlist' : 'craterlist.dat',
'sfdfile' : 'production.dat'
}
if self.user['sfdfile'] is not None: # Override the default sfdfile name if the user supplies an alternative
self.output_filenames['sfdfile'] = self.user['sfdfile']

for k, v in self.output_filenames.items():
self.output_filenames[k] = os.path.join(currentdir, v)


self.directories = ['dist', 'misc', 'surf']
if self.user['saveshaded'].upper() == 'T':
self.directories.append('shaded')
if self.user['saverego'].upper() == 'T':
self.directories.append('rego')

# Set up data arrays
self.seedarr = np.zeros(100, dtype=int)
self.seedarr[0] = self.user['seed']
Expand All @@ -79,6 +112,7 @@ def __init__(self, param_file="ctem.in", isnew=True):
self.tdist = np.zeros([1, 6])
self.surface_dem = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float)
self.surface_ejc = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float)
self.ph1 = None

if self.user['sfdcompare'] is not None:
# Read sfdcompare file
Expand All @@ -92,15 +126,50 @@ def __init__(self, param_file="ctem.in", isnew=True):
if (self.user['restart'].upper() == 'F'):
print('Starting a new run')

util.create_dir_structure(self.user)
util.create_dir_structure(self.user, self.directories)
# 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()


# Setup Quasi-MC run

if (self.user['quasimc'] == 'T'):

#Read list of real craters
print("quasi-MC mode is ON")
print("Generating the crater scaling data in CTEM")
rclist = util.read_formatted_ascii(self.user['realcraterlist'], skip_lines = 0)
tempfile = os.path.join(currentdir, 'temp.in')

# Generate craterlist.dat
shutil.copy2(self.user['ctemfile'], tempfile )

#Write a temporary input file to generate the necessary quasimc files
util.write_temp_input(self.user, tempfile)
util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
self.compute_one_interval(ctemin=tempfile)
os.remove(tempfile)

#Interpolate craterscale.dat to get impactor sizes from crater sizes given
df = pandas.read_csv(self.output_filenames['craterscale'], sep='\s+')
df['log(Dc)'] = np.log(df['Dcrat(m)'])
df['log(Di)'] = np.log(df['#Dimp(m)'])
xnew = df['log(Dc)'].values
ynew = df['log(Di)'].values
interp = interp1d(xnew, ynew, fill_value='extrapolate')
rclist[:,0] = np.exp(interp(np.log(rclist[:,0])))

#Convert age in Ga to "interval time"
rclist[:,5] = (self.user['interval'] * self.user['numintervals']) - craterproduction.Tscale(rclist[:,5], 'NPF_Moon')
rclist = rclist[rclist[:,5].argsort()]

#Export to dat file
util.write_realcraters(self.output_filenames['craterlist'], rclist)

util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
else:
print('Continuing a previous run')
Expand All @@ -124,7 +193,7 @@ def scale_production(self):
self.production[:, 1] = self.production[:, 1] * area * self.user['interval']

# Write the scaled production function to file
util.write_production(self.user, self.production)
util.write_production(self.output_filenames['sfdfile'], self.production)

return

Expand All @@ -143,12 +212,12 @@ def run(self):
self.redirect_outputs(['dat'], 'misc')

#Execute Fortran code
self.compute_one_interval()
self.compute_one_interval(self.user['ctemfile'])

# Read in output files
self.read_output()

# Process the o utput files
# Process the output files
self.process_output()

# Update ncount
Expand All @@ -173,27 +242,28 @@ def run(self):

return

def compute_one_interval(self):
def compute_one_interval(self, ctemin="ctem.in"):
"""
Executes the Fortran code to generate one interval of simulation output.
"""
# Create crater population and display CTEM progress on screen
print(self.user['ncount'], ' Calling FORTRAN routine')
try:
with subprocess.Popen([os.path.join(self.user['workingdir'], 'CTEM')],
p = subprocess.Popen([os.path.join(self.user['workingdir'], self.ctem_executable), ctemin],
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
universal_newlines=True) as p:
for line in p.stdout:
if "%" in line:
print(line.replace('\n','\r'), end='')
else:
print(line,end='')
res = p.communicate()
if p.returncode != 0:
for line in res[1]:
print(line, end='')
raise Exception ("CTEM Failure")
universal_newlines=True,
shell=False)
for line in p.stdout:
if "%" in line:
print(line.replace('\n','\r'), end='')
else:
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()
Expand Down Expand Up @@ -222,6 +292,7 @@ def read_output(self):

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


def process_output(self):
"""
Expand All @@ -231,16 +302,14 @@ def process_output(self):
# Display results
print(self.user['ncount'], ' Generating surface images and plots')

# Write surface dem, surface ejecta, shaded relief, and rplot data
# Write surface dem, surface ejecta and shaded relief 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
Expand All @@ -261,6 +330,8 @@ def process_output(self):
if (self.user['savetruelist'].upper() == 'T'):
self.redirect_outputs(['tcum'], 'dist')
self.redirect_outputs(['impmass'], 'misc')
if (self.user['saverego'].upper() == 'T') :
self.redirect_outputs(['stack','rego','age','melt','comp'], 'rego')



Expand Down Expand Up @@ -289,23 +360,16 @@ def cleanup(self):
"""
# 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)
util.destroy_dir_structure(self.user, self.directories)
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()
sim.run()
sim.run()

0 comments on commit 52160b5

Please sign in to comment.