diff --git a/ctem/NPF.py b/ctem/NPF.py new file mode 100644 index 00000000..cb7ee363 --- /dev/null +++ b/ctem/NPF.py @@ -0,0 +1,104 @@ +import numpy as np +# The Neukum production function +# Ivanov, Neukum, and Hartmann (2001) SSR v. 96 pp. 55-86 +def N1lunar(T): + return 5.44e-14 * (np.exp(6.93*T) - 1) + 8.38e-4 * T + + +def Nlunar(D): + # Lunar crater SFD + aL00 = -3.0876 + aL01 = -3.557528 + aL02 = 0.781027 + aL03 = 1.021521 + aL04 = -0.156012 + aL05 = -0.444058 + aL06 = 0.019977 + aL07 = 0.086850 + aL08 = -0.005874 + aL09 = -0.006809 + aL10 = 8.25e-4 + aL11 = 5.54e-5 + + return np.where((D > 0.01) & (D < 1000.0), + 10**(aL00 + + aL01 * np.log10(D) ** 1 + + aL02 * np.log10(D) ** 2 + + aL03 * np.log10(D) ** 3 + + aL04 * np.log10(D) ** 4 + + aL05 * np.log10(D) ** 5 + + aL06 * np.log10(D) ** 6 + + aL07 * np.log10(D) ** 7 + + aL08 * np.log10(D) ** 8 + + aL09 * np.log10(D) ** 9 + + aL10 * np.log10(D) ** 10 + + aL11 * np.log10(D) ** 11),float('nan')) + +def Rproj(D): + #Projectile SFD + aP00 = 0.0 + aP01 = -1.375458 + aP02 = 1.272521e-1 + aP03 = -1.282166 + aP04 = -3.074558e-1 + aP05 = 4.149280e-1 + aP06 = 1.910668e-1 + aP07 = -4.260980e-2 + aP08 = -3.976305e-2 + aP09 = -3.180179e-3 + aP10 = 2.799369e-3 + aP11 = 6.892223e-4 + aP12 = 2.614385e-6 + aP13 = -1.416178e-5 + aP14 = -1.191124e-6 + + return np.where((D > 1e-4) & (D < 300), + 10**(aP00 + + aP01 * np.log10(D)**1 + + aP02 * np.log10(D)**2 + + aP03 * np.log10(D)**3 + + aP04 * np.log10(D)**4 + + aP05 * np.log10(D)**5 + + aP06 * np.log10(D)**6 + + aP07 * np.log10(D)**7 + + aP08 * np.log10(D)**8 + + aP09 * np.log10(D)**9 + + aP10 * np.log10(D)**10 + + aP11 * np.log10(D)**11 + + aP12 * np.log10(D)**12 + + aP13 * np.log10(D)**13 + + aP14 * np.log10(D)**14),float('nan')) + +def N1mars(T): + # Mars crater SFD + # Ivanov (2001) SSR v. 96 pp. 87-104 + return 2.68e-14 * (np.exp(6.93 * T) - 1) + 4.13e-4 * T + +def Nmars(D): + aM00 = -3.384 + aM01 = -3.197 + aM02 = 1.257 + aM03 = 0.7915 + aM04 = -0.4861 + aM05 = -0.3630 + aM06 = 0.1016 + aM07 = 6.756e-2 + aM08 = -1.181e-2 + aM09 = -4.753e-3 + aM10 = 6.233e-4 + aM11 = 5.805e-5 + + return np.where((D > 0.01) & (D < 1000.), + 10**(aM00 + + aM01 * np.log10(D)**1 + + aM02 * np.log10(D)**2 + + aM03 * np.log10(D)**3 + + aM04 * np.log10(D)**4 + + aM05 * np.log10(D)**5 + + aM06 * np.log10(D)**6 + + aM07 * np.log10(D)**7 + + aM08 * np.log10(D)**8 + + aM09 * np.log10(D)**9 + + aM10 * np.log10(D)**10 + + aM11 * np.log10(D)**11), + np.full_like(D, np.nan)) diff --git a/ctem/__init__.py b/ctem/__init__.py new file mode 100644 index 00000000..2e8ce2a0 --- /dev/null +++ b/ctem/__init__.py @@ -0,0 +1 @@ +from ctem.driver import * \ No newline at end of file diff --git a/ctem/craterproduction.py b/ctem/craterproduction.py new file mode 100644 index 00000000..4d81b687 --- /dev/null +++ b/ctem/craterproduction.py @@ -0,0 +1,216 @@ +import numpy as np +from scipy import optimize + +# The Neukum production function for the Moon and Mars +#Lunar PF from: Ivanov, Neukum, and Hartmann (2001) SSR v. 96 pp. 55-86 +#Mars PF from: Ivanov (2001) SSR v. 96 pp. 87-104 +sfd_coef = { + "NPF_Moon" : [-3.0876,-3.557528,+0.781027,+1.021521,-0.156012,-0.444058,+0.019977,+0.086850,-0.005874,-0.006809,+8.25e-4, +5.54e-5], + "NPF_Mars" : [-3.384, -3.197, +1.257, +0.7915, -0.4861, -0.3630, +0.1016, +6.756e-2,-1.181e-2,-4.753e-3,+6.233e-4,+5.805e-5] + } +sfd_range = { + "NPF_Moon" : [0.01,1000], + "NPF_Mars" : [0.015,362] +} +#Exponential time constant (Ga) +tau = 6.93 +Nexp = 5.44e-14 + +time_coef = { + "NPF_Moon" : Nexp, + "NPF_Mars" : Nexp * 10**(sfd_coef.get("NPF_Mars")[0]) / 10**(sfd_coef.get("NPF_Moon")[0]) +} + +def N1(T, pfmodel): + """ + Return the N(1) value as a function of time for a particular production function model + + Parameters + ---------- + T : numpy array + Time in units of Ga + pfmodel : string + The production function model to use + Currently options are 'NPF_Moon' or 'NPF_Mars' + + Returns + ------- + N1 : numpy array + The number of craters per square kilometer greater than 1 km in diameter + """ + T_valid_range = np.where((T <= 4.5) & (T >= 0.0), T, np.nan) + return time_coef.get(pfmodel) * (np.exp(tau * T_valid_range) - 1.0) + 10 ** (sfd_coef.get(pfmodel)[0]) * T_valid_range + +def CSFD(Dkm,planet): + """ + Return the cumulative size-frequency distribution for a particular production function model + + Parameters + ---------- + Dkm : numpy array + Diameters in units of km + pfmodel : string + The production function model to use + Currently options are 'NPF_Moon' or 'NPF_Mars' + + Returns + ------- + CSFD : numpy array + The number of craters per square kilometer greater than Dkm in diameter at T=1 Ga + """ + logCSFD = sum(co * np.log10(Dkm) ** i for i, co in enumerate(sfd_coef.get(planet))) + return 10**(logCSFD) + +def DSFD(Dkm,planet): + """ + Return the differential size-frequency distribution for a particular production function model + + Parameters + ---------- + Dkm : numpy array + Diameters in units of km + pfmodel : string + The production function model to use + Currently options are 'NPF_Moon' or 'NPF_Mars' + + Returns + ------- + DSFD : numpy array + The differential number of craters (dN/dD) per square kilometer greater than Dkm in diameter at T = 1 Ga + """ + dcoef = sfd_coef.get(planet)[1:] + logDSFD = sum(co * np.log10(Dkm) ** i for i, co in enumerate(dcoef)) + return 10**(logDSFD) * CSFD(Dkm,planet) / Dkm + +def Tscale(T,planet): + """ + Return the number density of craters at time T relative to time T = 1 Ga + + Parameters + ---------- + T : numpy array + Time in units of Ga + pfmodel : string + The production function model to use + Currently options are 'NPF_Moon' or 'NPF_Mars' + + Returns + ------- + Tscale : numpy array + N1(T) / CSFD(Dkm = 1.0) + """ + return N1(T,planet) / CSFD(1.0,planet) + +def pf_csfd(T,Dkm,planet): + """ + Return the cumulative size-frequency distribution for a particular production function model as a function of Time + + Parameters + ---------- + T : numpy array + Time in units of Ga + Dkm : numpy array + Diameters in units of km + pfmodel : string + The production function model to use + Currently options are 'NPF_Moon' or 'NPF_Mars' + + Returns + ------- + pf_csfd : numpy array + The cumulative number of craters per square kilometer greater than Dkm in diameter at time T T + """ + D_valid_range = np.where((Dkm >= sfd_range.get(planet)[0]) & (Dkm <= sfd_range.get(planet)[1]),Dkm,np.nan) + return CSFD(D_valid_range,planet) * Tscale(T,planet) + +def pf_dsfd(T,Dkm,planet): + """ + Return the differential size-frequency distribution for a particular production function model as a function of Time + + Parameters + ---------- + T : numpy array + Time in units of Ga + Dkm : numpy array + Diameters in units of km + pfmodel : string + The production function model to use + Currently options are 'NPF_Moon' or 'NPF_Mars' + + Returns + ------- + pf_dsfd : numpy array + The cumulative number of craters per square kilometer greater than Dkm in diameter at time T T + """ + D_valid_range = np.where((Dkm >= sfd_range.get(planet)[0]) & (Dkm <= sfd_range.get(planet)[1]), Dkm, np.nan) + return DSFD(D_valid_range, planet) * Tscale(T, planet) + +def T_from_scale(TS,planet): + """ + Return the time in Ga for the given number density of craters relative to that at 1 Ga. + This is the inverse of Tscale + + Parameters + ---------- + TS : numpy array + number density of craters relative to that at 1 Ga + pfmodel : string + The production function model to use + Currently options are 'NPF_Moon' or 'NPF_Mars' + + Returns + ------- + T_from_scale : numpy array + The time in Ga + """ + def func(S): + return Tscale(S,planet) - TS + return optimize.fsolve(func, np.full_like(TS,4.4),xtol=1e-10) + + +if __name__ == "__main__": + import matplotlib.pyplot as plt + import matplotlib.ticker as ticker + print("Tests go here") + print(f"T = 1 Ga, N(1) = {pf_csfd(1.0,1.00,'NPF_Moon')}") + print(f"T = 4.2 Ga, N(1) = {pf_csfd(4.2,1.00,'NPF_Moon')}") + print("Tscale test: Should return all 1s") + Ttest = np.logspace(-4,np.log10(4.4),num=100) + Tres = T_from_scale(Tscale(Ttest,'NPF_Mars'),'NPF_Mars') + print(Ttest / Tres) + #for i,t in enumerate(Ttest): + # print(t,Tscale(t,'Moon'),Tres[i]) + + CSFDfig = plt.figure(1, figsize=(8, 7)) + ax = {'NPF_Moon': CSFDfig.add_subplot(121), + 'NPF_Mars': CSFDfig.add_subplot(122)} + + tvals = [0.01,1.0,4.0] + x_min = 1e-3 + x_max = 1e3 + y_min = 1e-9 + y_max = 1e3 + Dvals = np.logspace(np.log10(x_min), np.log10(x_max)) + for key in ax: + ax[key].title.set_text(key) + ax[key].set_xscale('log') + ax[key].set_yscale('log') + ax[key].set_ylabel('$\mathregular{N_{>D} (km^{-2})}$') + ax[key].set_xlabel('Diameter (km)') + ax[key].set_xlim(x_min, x_max) + ax[key].set_ylim(y_min, y_max) + ax[key].yaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=20)) + ax[key].yaxis.set_minor_locator(ticker.LogLocator(base=10.0, subs=np.arange(2,10), numticks=100)) + ax[key].xaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=20)) + ax[key].xaxis.set_minor_locator(ticker.LogLocator(base=10.0, subs=np.arange(2,10), numticks=100)) + ax[key].grid(True,which="minor",ls="-",lw=0.5,zorder=5) + ax[key].grid(True,which="major",ls="-",lw=1,zorder=10) + for t in tvals: + prod = pf_csfd(t,Dvals,key) + ax[key].plot(Dvals, prod, '-', color='black', linewidth=1.0, zorder=50) + labeli = 15 + ax[key].text(Dvals[labeli],prod[labeli],f"{t:.2f} Ga", ha="left", va="top",rotation=-72) + + plt.tick_params(axis='y', which='minor') + plt.tight_layout() + plt.show() diff --git a/ctem/driver.py b/ctem/driver.py new file mode 100644 index 00000000..2646d1db --- /dev/null +++ b/ctem/driver.py @@ -0,0 +1,383 @@ +import numpy as np +import os +import subprocess +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: + """ + This is a class that defines the basic CTEM simulation object. It initializes a dictionary of user and reads + it in from a file (default name is ctem.in). It also creates the directory structure needed to store simulation + files if necessary. + """ + def __init__(self, param_file="ctem.in", isnew=True): + currentdir = os.getcwd() + self.user = { + 'restart': None, + 'runtype': None, + 'popupconsole': None, + 'saveshaded': None, + 'saverego': None, + 'savetruelist': None, + 'seedn': 1, + 'totalimpacts': 0, + 'ncount': 0, + 'curyear': 0.0, + 'fracdone': 1.0, + 'masstot': 0.0, + 'interval': 0.0, + 'numintervals': 0, + 'pix': -1.0, + 'gridsize': -1, + 'seed': 0, + 'maxcrat': 1.0, + 'shadedminhdefault': 1, + 'shadedmaxhdefault': 1, + 'shadedminh': 0.0, + 'shadedmaxh': 0.0, + 'workingdir': currentdir, + 'ctemfile': param_file, + 'impfile': None, + 'sfdcompare': None, + 'quasimc': None, + 'sfdfile' : None, + 'realcraterlist': None, + } + + # Get the location of the CTEM executable + self.ctem_executable = Path(currentdir) / "CTEM" + if not self.ctem_executable.exists(): + print(f"CTEM driver not found at {self.ctem_executable}. Trying package directory..") + self.ctem_executable = Path(_pyfile).parent.parent.parent / "bin" / "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 = { + 'dat': 'ctem.dat', + 'dem': 'surface_dem.dat', + 'diam': 'surface_diam.dat', + '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', + 'tcum' : 'tcumulative.dat', + 'tdist' : 'tdistribution.dat', + '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', + 'craterlist' : 'craterlist.dat', + 'sfdfile' : 'production.dat', + 'distfrac' : 'surface_distfrac.dat', + 'ejm' : 'surface_ejm.dat', + 'ejmf' : 'surface_ejmf.dat', + 'meltdist' : 'surface_meltdist.dat', + 'meltfrac' : 'surface_meltfrac.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'] + self.odist = np.zeros([1, 6]) + 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.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 + sfdfile = os.path.join(self.user['workingdir'], self.user['sfdcompare']) + self.ph1 = util.read_formatted_ascii(sfdfile, skip_lines=0) + + + # 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, 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" + if (self.user['runtype'].upper() == 'STATISTICAL'): + rclist[:,5] = (self.user['interval']) - craterproduction.Tscale(rclist[:,5], 'NPF_Moon') + else: + 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') + self.read_output() + + return + + def scale_production(self): + """ + Scales the production function to the simulation domain and saves the scaled production function to a file that + will be read in by the main Fortran program. + """ + + # Read production function file + impfile = os.path.join(self.user['workingdir'], self.user['impfile']) + 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 + util.write_production(self.output_filenames['sfdfile'], self.production) + + return + + def run(self): + """ + Runs a complete simulation over all intervals specified in the input user. This method loops over all + intervals and process outputs into images, and redirect output files to their apporpriate folders. + + This method replaces most of the functionality of the original ctem_driver.py + """ + + print('Beginning loops') + while (self.user['ncount'] <= self.user['numintervals']): + if (self.user['ncount'] > 0): + # Move ctem.dat + self.redirect_outputs(['dat'], 'misc') + + #Execute Fortran code + self.compute_one_interval(self.user['ctemfile']) + + # Read in output files + self.read_output() + + # Process the output files + self.process_output() + + # Update ncount + self.user['ncount'] = 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 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 for next interval + util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr) + + return + + 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: + p = subprocess.Popen([os.path.join(self.user['workingdir'], self.ctem_executable), ctemin], + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + 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() + + return + + 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. + """ + # Read Fortran output + print(self.user['ncount'], ' Reading Fortran output') + + # Read surface dem(shaded relief) and ejecta data files + 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 = 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 + 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') + + # 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 + + # 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') + if (self.user['saverego'].upper() == 'T') : + self.redirect_outputs(['stack','rego','age','melt','comp','ejm','meltdist'], 'rego') + + + + def redirect_outputs(self, filekeys, foldername): + """ + Copies a set of output files from the working directory into a subfolder. + Takes as input the dictionaries of user parameters and output file names, the dictionary keys to the file names + you wish to redirect, and the redirection destination folder name. The new name will be the key name + zero padded + interval number. + + Example: + calling redirect_outputs(['odist','tdist'], 'dist') when user['ncount'] is 1 + should do the following: + copy 'odistribution.dat' to 'dist/odist_000001.dat' + copy 'tdistribution.dat' to 'dist/tdist_000001.dat' + """ + + for k in filekeys: + 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, 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) + + return + +if __name__ == '__main__': + sim = Simulation() + sim.run() \ No newline at end of file diff --git a/ctem/util.py b/ctem/util.py new file mode 100644 index 00000000..00922421 --- /dev/null +++ b/ctem/util.py @@ -0,0 +1,443 @@ +import numpy as np +import os +import shutil +from matplotlib.colors import LightSource +import matplotlib.cm as cm +import matplotlib.pyplot as plt +import re +from tempfile import mkstemp +from scipy.io import FortranFile + +# Set pixel scaling common for image writing, at 1 pixel/ array element +dpi = 72.0 +class CTEMLightSource(LightSource): + """ + Override one function in LightSource to prevent the contrast from being rescaled. + """ + def shade_normals(self, normals, fraction=1.): + """ + Calculate the illumination intensity for the normal vectors of a + surface using the defined azimuth and elevation for the light source. + + Imagine an artificial sun placed at infinity in some azimuth and + elevation position illuminating our surface. The parts of the surface + that slope toward the sun should brighten while those sides facing away + should become darker. + + Changes by David Minton: The matplotlib version of this rescales the intensity + in a way that causes the brightness level of the images to change between + simulation outputs. This makes movies of the surface evolution appear to flicker. + + Parameters + ---------- + fraction : number, optional + Increases or decreases the contrast of the hillshade. Values + greater than one will cause intermediate values to move closer to + full illumination or shadow (and clipping any values that move + beyond 0 or 1). Note that this is not visually or mathematically + the same as vertical exaggeration. + + Returns + ------- + ndarray + A 2D array of illumination values between 0-1, where 0 is + completely in shadow and 1 is completely illuminated. + """ + + intensity = normals.dot(self.direction) + + # Apply contrast stretch + imin, imax = intensity.min(), intensity.max() + intensity *= fraction + + # Rescale to 0-1, keeping range before contrast stretch + # If constant slope, keep relative scaling (i.e. flat should be 0.5, + # fully occluded 0, etc.) + # if (imax - imin) > 1e-6: + # # Strictly speaking, this is incorrect. Negative values should be + # # clipped to 0 because they're fully occluded. However, rescaling + # # in this manner is consistent with the previous implementation and + # # visually appears better than a "hard" clip. + # intensity -= imin + # intensity /= (imax - imin) + intensity = np.clip(intensity, 0, 1) + + return intensity + + + +# These are directories that are created by CTEM in order to store intermediate results + +def create_dir_structure(user, directories): + """ + 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, directories): + """ + 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 image_dem(user, DEM): + dpi = 300.0 # 72.0 + pix = user['pix'] + gridsize = user['gridsize'] + ve = 1.0 + azimuth = 300.0 # user['azimuth'] + solar_angle = 20.0 # user['solar_angle'] + + ls = CTEMLightSource(azdeg=azimuth, altdeg=solar_angle) + dem_img = ls.hillshade(np.flip(DEM,axis=0), vert_exag=ve, dx=pix, dy=pix, fraction=1) + + # Generate image to put into an array + height = gridsize / dpi + width = gridsize / dpi + fig = plt.figure(figsize=(width, height), dpi=dpi) + ax = plt.axes([0, 0, 1, 1]) + ax.imshow(dem_img, interpolation="nearest", cmap='gray', vmin=0.0, vmax=1.0) + plt.axis('off') + # Save image to file + filename = os.path.join(user['workingdir'],'surf',"surf%06d.png" % user['ncount']) + plt.savefig(filename, dpi=dpi, bbox_inches=0) + plt.close() + + return + + +def image_regolith(user, regolith): + # Create scaled regolith image + minref = user['pix'] * 1.0e-4 + maxreg = np.amax(regolith) + minreg = np.amin(regolith) + if (minreg < minref): minreg = minref + if (maxreg < minref): maxreg = (minref + 1.0e3) + regolith_scaled = np.copy(regolith) + np.place(regolith_scaled, regolith_scaled < minref, minref) + regolith_scaled = 254.0 * ( + (np.log(regolith_scaled) - np.log(minreg)) / (np.log(maxreg) - np.log(minreg))) + + # Save image to file + filename = os.path.join(user['workingdir'],'rego', "rego%06d.png" % user['ncount']) + height = user['gridsize'] / dpi + width = height + fig = plt.figure(figsize=(width, height), dpi=dpi) + fig.figimage(regolith_scaled, cmap=cm.nipy_spectral, origin='lower') + plt.savefig(filename) + plt.close() + + return + + +def image_shaded_relief(user, DEM): + dpi = 300.0 # 72.0 + pix = user['pix'] + gridsize = user['gridsize'] + ve = 1.0 + mode = 'overlay' + azimuth = 300.0 # user['azimuth'] + solar_angle = 20.0 # user['solar_angle'] + + ls = CTEMLightSource(azdeg=azimuth, altdeg=solar_angle) + cmap = cm.cividis + + # If min and max appear to be reversed, then fix them + if (user['shadedminh'] > user['shadedmaxh']): + temp = user['shadedminh'] + user['shadedminh'] = user['shadedmaxh'] + user['shadedmaxh'] = temp + else: + user['shadedminh'] = user['shadedminh'] + user['shadedmaxh'] = user['shadedmaxh'] + + # If no shadedmin/max user are read in from ctem.dat, determine the values from the data + if (user['shadedminhdefault'] == 1): + shadedminh = np.amin(DEM) + else: + shadedminh = user['shadedminh'] + if (user['shadedmaxhdefault'] == 1): + shadedmaxh = np.amax(DEM) + else: + shadedmaxh = user['shadedmaxh'] + + dem_img = ls.shade(np.flip(DEM,axis=0), cmap=cmap, blend_mode=mode, fraction=1.0, + vert_exag=ve, dx=pix, dy=pix, + vmin=shadedminh, vmax=shadedmaxh) + + # Generate image to put into an array + height = gridsize / dpi + width = gridsize / dpi + fig = plt.figure(figsize=(width, height), dpi=dpi) + ax = plt.axes([0, 0, 1, 1]) + ax.imshow(dem_img, interpolation="nearest", vmin=0.0, vmax=1.0) + plt.axis('off') + # Save image to file + filename = os.path.join(user['workingdir'],'shaded',"shaded%06d.png" % user['ncount']) + plt.savefig(filename, dpi=dpi, bbox_inches=0) + plt.close() + + return user + + +def read_datfile(user, datfile, seedarr): + # Read and parse ctem.dat file + + # Read ctem.dat file + print('Reading input file ' + datfile) + fp = open(datfile, 'r') + lines = fp.readlines() + fp.close() + + # Parse file lines and update parameter fields + fields = lines[0].split() + if len(fields) > 0: + user['totalimpacts'] = real2float(fields[0]) + user['ncount'] = int(fields[1]) + user['curyear'] = real2float(fields[2]) + user['restart'] = fields[3] + user['fracdone'] = real2float(fields[4]) + user['masstot'] = real2float(fields[5]) + + # Parse remainder of file to build seed array + nlines = len(lines) + index = 1 + while (index < nlines): + fields = lines[index].split() + seedarr[index - 1] = real2float(fields[0]) + index += 1 + + user['seedn'] = index - 1 + + return + + +def read_formatted_ascii(filename, skip_lines): + # Generalized ascii text reader + # For use with sfdcompare, production, odist, tdist, pdist data files + try: + data = np.genfromtxt(filename, skip_header=skip_lines) + except: + print(f"Error reading {filename}") + data = None + return data + + +def read_impact_mass(filename): + # Read impact mass file + + fp = open(filename, 'r') + line = fp.readlines() + fp.close() + + fields = line[0].split() + if (len(fields) > 0): + mass = real2float(fields[0]) + else: + mass = 0 + + return mass + + +# Write production function to file production.dat +# This file format does not exactly match that generated from IDL. Does it work? +def read_user_input(user): + # Read and parse ctem.in file + inputfile = user['ctemfile'] + + # Read ctem.in file + print('Reading input file ' + user['ctemfile']) + with open(inputfile, 'r') as fp: + lines = fp.readlines() + + # Process file text + for line in lines: + fields = line.split() + if len(fields) > 0: + if ('pix' == fields[0].lower()): user['pix'] = real2float(fields[1]) + if ('gridsize' == fields[0].lower()): user['gridsize'] = int(fields[1]) + if ('seed' == fields[0].lower()): user['seed'] = int(fields[1]) + if ('sfdfile' == fields[0].lower()): user['sfdfile'] = os.path.join(user['workingdir'],fields[1]) + if ('impfile' == fields[0].lower()): user['impfile'] = os.path.join(user['workingdir'],fields[1]) + if ('maxcrat' == fields[0].lower()): user['maxcrat'] = real2float(fields[1]) + if ('sfdcompare' == fields[0].lower()): user['sfdcompare'] = os.path.join(user['workingdir'], fields[1]) + if ('realcraterlist' == fields[0].lower()): user['realcraterlist'] = os.path.join(user['workingdir'], fields[1]) + if ('interval' == fields[0].lower()): user['interval'] = real2float(fields[1]) + if ('numintervals' == fields[0].lower()): user['numintervals'] = int(fields[1]) + if ('popupconsole' == fields[0].lower()): user['popupconsole'] = fields[1] + if ('saveshaded' == fields[0].lower()): user['saveshaded'] = fields[1] + if ('saverego' == fields[0].lower()): user['saverego'] = fields[1] + if ('savepres' == fields[0].lower()): user['savepres'] = fields[1] + if ('savetruelist' == fields[0].lower()): user['savetruelist'] = fields[1] + if ('runtype' == fields[0].lower()): user['runtype'] = fields[1] + if ('restart' == fields[0].lower()): user['restart'] = fields[1] + if ('quasimc' == fields[0].lower()): user['quasimc'] = fields[1] + if ('shadedminh' == fields[0].lower()): + user['shadedminh'] = real2float(fields[1]) + user['shadedminhdefault'] = 0 + if ('shadedmaxh' == fields[0].lower()): + user['shadedmaxh'] = real2float(fields[1]) + user['shadedmaxhdefault'] = 0 + + # Test values for further processing + if (user['interval'] <= 0.0): + print('Invalid value for or missing variable INTERVAL in ' + inputfile) + if (user['numintervals'] <= 0): + print('Invalid value for or missing variable NUMINTERVALS in ' + inputfile) + if (user['pix'] <= 0.0): + print('Invalid value for or missing variable PIX in ' + inputfile) + if (user['gridsize'] <= 0): + print('Invalid value for or missing variable GRIDSIZE in ' + inputfile) + if (user['seed'] == 0): + print('Invalid value for or missing variable SEED in ' + inputfile) + if (user['impfile'] is None): + print('Invalid value for or missing variable IMPFILE in ' + inputfile) + if (user['popupconsole'] is None): + print('Invalid value for or missing variable POPUPCONSOLE in ' + inputfile) + if (user['saveshaded'] is None): + print('Invalid value for or missing variable SAVESHADED in ' + inputfile) + if (user['saverego'] is None): + print('Invalid value for or missing variable SAVEREGO in ' + inputfile) + if (user['savepres'] is None): + print('Invalid value for or missing variable SAVEPRES in ' + inputfile) + if (user['savetruelist'] is None): + print('Invalid value for or missing variable SAVETRUELIST in ' + inputfile) + if (user['runtype'] is None): + print('Invalid value for or missing variable RUNTYPE in ' + inputfile) + if (user['restart'] is None): + print('Invalid value for or missing variable RESTART in ' + inputfile) + + return user + + +def read_unformatted_binary(filename, gridsize, kind='DP'): + # Read unformatted binary files created by Fortran + # For use with surface ejecta and surface dem data files + if kind == 'DP': + dt = np.dtype('f8') + elif kind == 'SP': + dt = np.dtype('f4') + elif kind == 'I4B': + dt = np.dtype(' count: + break + + fout.writelines(fin.readlines()) + + + fin.close() + fout.close() + + shutil.move(name, source) + + return + +def write_datfile(user, filename, seedarr): + # Write various user and random number seeds into ctem.dat file + fp = open(filename, 'w') + + template = "%(totalimpacts)17d %(ncount)12d %(curyear)19.12E %(restart)s %(fracdone)9.6f %(masstot)19.12E\n" + fp.write(template % user) + + # Write random number seeds to the file + for index in range(user['seedn']): + fp.write("%12d\n" % seedarr[index]) + + fp.close() + + return + + +def write_production(filename, production): + np.savetxt(filename, production, fmt='%1.8e', delimiter=' ') + + return + + +def write_realcraters(filename, realcraters): + """Writes file of real craters for use in quasi-MC runs""" + + np.savetxt(filename, realcraters, fmt='%1.8e', delimiter='\t') + + return + +def write_temp_input(user, filename): + """Makes changes to a temporary input file for use when generating craterlist.dat for quasimc runs""" + + sed('testflag', 'testflag T!', filename) + sed('testimp', f'testimp {user["pix"]*1e-3} !', filename) # Make a tiny test crater. We don't care about the crater itself, just that we run CTEM once to get all of the converted files + sed('quasimc', 'quasimc F!', filename) + sed('interval', 'interval 1 !', filename) + sed('numinterval 1 !s', 'numintervals 1 !', filename) + + return \ No newline at end of file diff --git a/ctem/viewer3d.py b/ctem/viewer3d.py new file mode 100644 index 00000000..8aacf84e --- /dev/null +++ b/ctem/viewer3d.py @@ -0,0 +1,124 @@ +import numpy as np +import open3d +import ctem + +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() + #used for Open3d module + self.mesh = open3d.geometry.TriangleMesh() + return + + def ctem2blockmesh(self): + # Build mesh grid + s = self.user['gridsize'] + pix = self.user['pix'] + ygrid, xgrid = np.mgrid[-s/2:s/2, -s/2:s/2] * pix + y0, x0 = ygrid - pix/2, xgrid - pix/2 + y1, x1 = ygrid - pix/2, xgrid + pix/2 + y2, x2 = ygrid + pix/2, xgrid + pix/2 + y3, x3 = ygrid + pix/2, xgrid - pix/2 + + xvals = np.append( + np.append( + np.append(x0.flatten(), + x1.flatten()), + x2.flatten()), + x3.flatten()) + yvals = np.append( + np.append( + np.append(y0.flatten(), + y1.flatten()), + y2.flatten()), + y3.flatten()) + zvals = np.append( + np.append( + np.append(self.surface_dem.flatten(), + self.surface_dem.flatten()), + self.surface_dem.flatten()), + self.surface_dem.flatten()) + verts = np.array((xvals, yvals, zvals)).T + nface_triangles = 10 + faces = np.full([nface_triangles*s**2, 3], -1, dtype=np.int64) + for j in range(s): + for i in range(s): + i0 = s*j + i + i1 = i0 + s**2 + i2 = i1 + s**2 + i3 = i2 + s**2 + + fidx = np.arange(nface_triangles*i0,nface_triangles*(i0+1)) + # Make the two top triangles + faces[fidx[0],:] = [i0, i1, i2] + faces[fidx[1],:] = [i0, i2, i3] + # Make the two west side triangles + if i > 0: + faces[fidx[2],:] = [i0, i3, i2-1] + faces[fidx[3],:] = [i0, i2-1, i1-1] + # Make the two south side triangles + if j > 0: + faces[fidx[4],:] = [i1, i0, i3-s ] + faces[fidx[5],:] = [i1, i3-s, i2-s] + # Make the two east side triangles + if i < (s - 1): + faces[fidx[6],:] = [i2, i1, i0+1] + faces[fidx[7],:] = [i2, i0+1, i3+1] + #Make the two north side triangles + if j < (s -1): + faces[fidx[8],:] = [i3, i2, i1+s] + faces[fidx[9],:] = [i3, i1+s, i0+s] + nz = faces[:,0] != -1 + f2 = faces[nz,:] + self.mesh.vertices = open3d.utility.Vector3dVector(verts) + self.mesh.triangles = open3d.utility.Vector3iVector(f2) + self.mesh.compute_vertex_normals() + self.mesh.compute_triangle_normals() + + return + + def ctem2trimesh(self): + # Build mesh grid + s = self.user['gridsize'] + pix = self.user['pix'] + ygrid, xgrid = np.mgrid[-s/2:s/2, -s/2:s/2] * pix + + xvals = xgrid.flatten() + yvals = ygrid.flatten() + zvals = self.surface_dem.flatten() + verts = np.array((xvals, yvals, zvals)).T + faces = np.full([2 * (s-1)**2, 3], -1, dtype=np.int64) + for j in range(s - 1): + for i in range(s - 1): + idx = (s - 1) * j + i + faces[idx, :] = [j * s + i, j * s + i + 1, (j + 1) * s + i + 1] + idx += (s - 1) ** 2 + faces[idx, :] = [(j + 1) * s + i + 1, (j + 1) * s + i, j * s + i ] + self.mesh.vertices = open3d.utility.Vector3dVector(verts) + self.mesh.triangles = open3d.utility.Vector3iVector(faces) + self.mesh.compute_vertex_normals() + self.mesh.compute_triangle_normals() + + return + + def visualize(self): + vis = open3d.visualization.Visualizer() + vis.create_window() + vis.add_geometry(self.mesh) + opt = vis.get_render_option() + opt.background_color = np.asarray([0, 0, 0]) + + self.mesh.paint_uniform_color([0.5, 0.5, 0.5]) + vis.run() + vis.destroy_window() + + +if __name__ == '__main__': + sim = Polysurface() + sim.surface_dem *= 10 + sim.ctem2trimesh() + sim.visualize() + + +