From b8802390193d4ff018636e4726027bdfb8f6528f Mon Sep 17 00:00:00 2001 From: David Minton Date: Sat, 10 Feb 2024 18:07:21 -0500 Subject: [PATCH] Restructured python files and fixed bugs to make the code more compatible with ifx --- .gitignore | 1 + cmake/Modules/SetFortranFlags.cmake | 16 +- ctem/ctem/NPF.py | 104 ------ ctem/ctem/__init__.py | 1 - ctem/ctem/craterproduction.py | 216 ----------- ctem/ctem/driver.py | 383 -------------------- ctem/ctem/util.py | 443 ----------------------- ctem/ctem/viewer3d.py | 124 ------- ctem/setup.py | 6 - src/crater/crater_subpixel_diffusion.f90 | 8 +- src/ejecta/ejecta_emplace.f90 | 4 +- src/init/init_dist.f90 | 12 +- src/io/io_ejecta_table.f90 | 2 +- src/io/io_input.f90 | 2 +- src/io/io_read_craterlist.f90 | 6 +- src/io/io_read_prod.f90 | 2 +- src/io/io_read_vdist.f90 | 6 +- 17 files changed, 31 insertions(+), 1305 deletions(-) delete mode 100644 ctem/ctem/NPF.py delete mode 100644 ctem/ctem/__init__.py delete mode 100644 ctem/ctem/craterproduction.py delete mode 100644 ctem/ctem/driver.py delete mode 100644 ctem/ctem/util.py delete mode 100644 ctem/ctem/viewer3d.py delete mode 100644 ctem/setup.py diff --git a/.gitignore b/.gitignore index 6ce9617f..70c42081 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,7 @@ !distclean.cmake !README.md !version.txt +!ctem/ !ctem/**.py !cmake/Modules/*.cmake !*.bash diff --git a/cmake/Modules/SetFortranFlags.cmake b/cmake/Modules/SetFortranFlags.cmake index edf73442..f42d8d75 100644 --- a/cmake/Modules/SetFortranFlags.cmake +++ b/cmake/Modules/SetFortranFlags.cmake @@ -533,10 +533,10 @@ IF (CMAKE_BUILD_TYPE STREQUAL "RELEASE" OR CMAKE_BUILD_TYPE STREQUAL "PROFILE") ) # Tells the compiler to link to certain libraries in the Intel oneAPI Math Kernel Library (oneMKL). - SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" - Fortran "/Qmkl:cluster" # Intel Windows - "/Qmkl" # Intel Windows - ) + # SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" + # Fortran "/Qmkl:cluster" # Intel Windows + # "/Qmkl" # Intel Windows + # ) # Enables additional interprocedural optimizations for a single file compilation SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" @@ -574,10 +574,10 @@ IF (CMAKE_BUILD_TYPE STREQUAL "RELEASE" OR CMAKE_BUILD_TYPE STREQUAL "PROFILE") ) # Tells the compiler to link to certain libraries in the Intel oneAPI Math Kernel Library (oneMKL). - SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" - Fortran "-qmkl=cluster" # Intel - "-qmkl" # Intel - ) + # SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" + # Fortran "-qmkl=cluster" # Intel + # "-qmkl" # Intel + # ) # Enables additional interprocedural optimizations for a single file compilation SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" diff --git a/ctem/ctem/NPF.py b/ctem/ctem/NPF.py deleted file mode 100644 index cb7ee363..00000000 --- a/ctem/ctem/NPF.py +++ /dev/null @@ -1,104 +0,0 @@ -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/ctem/__init__.py b/ctem/ctem/__init__.py deleted file mode 100644 index 2e8ce2a0..00000000 --- a/ctem/ctem/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from ctem.driver import * \ No newline at end of file diff --git a/ctem/ctem/craterproduction.py b/ctem/ctem/craterproduction.py deleted file mode 100644 index 4d81b687..00000000 --- a/ctem/ctem/craterproduction.py +++ /dev/null @@ -1,216 +0,0 @@ -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/ctem/driver.py b/ctem/ctem/driver.py deleted file mode 100644 index 2646d1db..00000000 --- a/ctem/ctem/driver.py +++ /dev/null @@ -1,383 +0,0 @@ -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/ctem/util.py b/ctem/ctem/util.py deleted file mode 100644 index 00922421..00000000 --- a/ctem/ctem/util.py +++ /dev/null @@ -1,443 +0,0 @@ -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/ctem/viewer3d.py b/ctem/ctem/viewer3d.py deleted file mode 100644 index 8aacf84e..00000000 --- a/ctem/ctem/viewer3d.py +++ /dev/null @@ -1,124 +0,0 @@ -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() - - - diff --git a/ctem/setup.py b/ctem/setup.py deleted file mode 100644 index f7003601..00000000 --- a/ctem/setup.py +++ /dev/null @@ -1,6 +0,0 @@ -from setuptools import setup, find_packages - -setup(name='ctem', - version='0.1', - author='David A. Minton', - packages=find_packages()) diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90 index cdba3408..63a69bfe 100644 --- a/src/crater/crater_subpixel_diffusion.f90 +++ b/src/crater/crater_subpixel_diffusion.f90 @@ -59,9 +59,6 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin) real(DP) :: rayfmult = (5)**(-4.0_DP / (1.2_DP)) real(DP) :: l1 - rray = 11.95*crater%frad**1.32 - l1 = 5.32*crater%frad**1.27 - ! Create box for soften calculation (will be no bigger than the grid itself) do j = 0,user%gridsize + 1 do i = 0,user%gridsize + 1 @@ -187,7 +184,10 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin) allocate(diffdistribution(imin:imax,jmin:jmax)) allocate(ejdistribution(imin:imax,jmin:jmax)) - call ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution,l1) + rray = 11.95_DP*crater%frad**1.32 + l1 = 5.32_DP*crater%frad**1.27 + + call ejecta_ray_pattern(user,surf,crater,inc,imin,imax,jmin,jmax,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution,l1) ! Loop over affected matrix area !!$OMP PARALLEL DO DEFAULT(SHARED) IF(inc > INCPAR) & !!$OMP FIRSTPRIVATE(jmin,jmax,imin,imax) & diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index 5e6de56c..9cc16e41 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -254,7 +254,9 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ klo = int((log(lrad) - log(crater%ejrad)) / domain%ejbres) do n = 1,MAXLOOP call ejecta_interpolate(crater,domain,distance,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,erad=erad,melt=melt) - if ((n > 1).and.((abs(ebh0 - ebh) / ebh0) < domain%small)) exit + if (n > 1) then + if ((abs(ebh0 - ebh) / ebh0) < domain%small) exit + endif ebh0 = ebh erad = exp(erad) diff --git a/src/init/init_dist.f90 b/src/init/init_dist.f90 index fc17a1bd..72e96af5 100644 --- a/src/init/init_dist.f90 +++ b/src/init/init_dist.f90 @@ -33,7 +33,7 @@ subroutine init_dist(user,domain) ! Executable code ! Get size of crater production function array - open(unit=LUN,file=user%sfdfile,status="old",iostat=ierr) + open(unit=LUN,file=trim(adjustl(user%sfdfile)),status="old",iostat=ierr) if (ierr /= 0) then write(*,*) "Unable to open file ", trim(user%sfdfile) stop @@ -47,12 +47,12 @@ subroutine init_dist(user,domain) end do if (domain%pnum == 0) then - write(*,*) "No valid entries in ",trim(user%sfdfile) + write(*,*) "No valid entries in ",trim(adjustl(user%sfdfile)) end if close(LUN) ! Get size of velocity distribution array - open(unit=LUN,file=user%velfile,status="old",iostat=ierr) + open(unit=LUN,file=trim(adjustl(user%velfile)),status="old",iostat=ierr) if (ierr /= 0) then write(*,*) "Unable to open file ",trim(user%velfile) stop @@ -65,15 +65,15 @@ subroutine init_dist(user,domain) domain%vnum = domain%vnum + 1 end do if (domain%vnum == 0) then - write(*,*) "No valid entries in ",trim(user%velfile) + write(*,*) "No valid entries in ",trim(adjustl(user%velfile)) end if close(LUN) ! Get size of real crater list array if (user%doquasimc) then - open(unit=LUN,file=rcfile,status="old",iostat=ierr) + open(unit=LUN,file=trim(adjustl(rcfile)),status="old",iostat=ierr) if (ierr /= 0) then - write(*,*) "Unable to open file ", trim(rcfile) + write(*,*) "Unable to open file ", trim(adjustl(rcfile)) stop end if diff --git a/src/io/io_ejecta_table.f90 b/src/io/io_ejecta_table.f90 index f5e1287a..342af640 100644 --- a/src/io/io_ejecta_table.f90 +++ b/src/io/io_ejecta_table.f90 @@ -35,7 +35,7 @@ subroutine io_ejecta_table(crater,domain,ejb,ejtble,filename) integer(I4B) :: k ! Executable code - open(LUN, FILE=filename, status='replace') + open(LUN, FILE=trim(adjustl(filename)), status='replace') write(LUN,'("# trad = ",ES12.5, " frad = ",ES12.5)') crater%rad,crater%frad write(LUN,'("# ejrim = ",ES12.5, " ejdis = ",ES12.5," imp = ",ES12.5)') crater%ejrim,crater%ejdis,crater%imp write(LUN,'(A63)') '# "r (m)" "h (m)" "v (m/s)" "ang (deg)" "erad (m)"' diff --git a/src/io/io_input.f90 b/src/io/io_input.f90 index b89b2612..4a528246 100644 --- a/src/io/io_input.f90 +++ b/src/io/io_input.f90 @@ -97,7 +97,7 @@ subroutine io_input(infile,user) user%dotopodiffusion = .true. write(user%sfdfile,*) trim(adjustl(SFDFILE)) - open(unit=LUN,file=infile,status="old",iostat=ierr) + open(unit=LUN,file=trim(adjustl(infile)),status="old",iostat=ierr) if (ierr /= 0) then write(*,*) "Unable to open file ",trim(infile) stop diff --git a/src/io/io_read_craterlist.f90 b/src/io/io_read_craterlist.f90 index ec17ff58..35b4e38d 100644 --- a/src/io/io_read_craterlist.f90 +++ b/src/io/io_read_craterlist.f90 @@ -35,16 +35,16 @@ subroutine io_read_craterlist(rclist, user, domain) ! Read the next crater from the craterlist - open(unit=LUN,file=rcfile,status='old',iostat=ierr) + open(unit=LUN,file=trim(adjustl(rcfile)),status='old',iostat=ierr) if (ierr /= 0) then - write(*,*) "Unable to open file ",trim(rcfile) + write(*,*) "Unable to open file ",trim(adjustl(rcfile)) stop end if do i=1,domain%rcnum read(LUN,*,iostat=ierr) rclist(1:6,i) if (ierr/=0) then - write(*,*) "Unable to read file ",trim(rcfile) + write(*,*) "Unable to read file ",trim(adjustl(rcfile)) stop end if end do diff --git a/src/io/io_read_prod.f90 b/src/io/io_read_prod.f90 index 2eb47516..8a845980 100644 --- a/src/io/io_read_prod.f90 +++ b/src/io/io_read_prod.f90 @@ -35,7 +35,7 @@ subroutine io_read_prod(prod,user,domain) ! Executable code ! Read in the size-frequency distribution file - open(unit=LUN,file=user%sfdfile,status='old',iostat=ierr) + open(unit=LUN,file=trim(adjustl(user%sfdfile)),status='old',iostat=ierr) if (ierr /= 0) then write(*,*) "Unable to open file ",trim(user%sfdfile) stop diff --git a/src/io/io_read_vdist.f90 b/src/io/io_read_vdist.f90 index a15dc8dd..84309feb 100644 --- a/src/io/io_read_vdist.f90 +++ b/src/io/io_read_vdist.f90 @@ -35,16 +35,16 @@ subroutine io_read_vdist(vdist,user,domain) ! Executable code ! Read in velocity distribution file - open(unit=LUN,file=user%velfile,status='old',iostat=ierr) + open(unit=LUN,file=trim(adjustl(user%velfile)),status='old',iostat=ierr) if (ierr /= 0) then - write(*,*) "Unable to open file ",trim(user%velfile) + write(*,*) "Unable to open file ",trim(adjustl(user%velfile)) stop end if do i=1,domain%vnum read(LUN,*,iostat=ierr) vdist(1:3,i) if (ierr/=0) then - write(*,*) "Unable to read file ",trim(user%velfile) + write(*,*) "Unable to read file ",trim(adjustl(user%velfile)) stop end if end do