diff --git a/examples/global-lunar-bombardment/ctem.in b/examples/global-lunar-bombardment/ctem.in index 355c22d3..470aaee4 100755 --- a/examples/global-lunar-bombardment/ctem.in +++ b/examples/global-lunar-bombardment/ctem.in @@ -70,6 +70,7 @@ Kd1 0.0001 psi 2.000 fe 4.00 ejecta_truncation 4.0 +doregotrack F dorays T superdomain F dorealistic F diff --git a/examples/morphology_test_cases/global/craterproduction.py b/examples/morphology_test_cases/global/craterproduction.py new file mode 100644 index 00000000..4d81b687 --- /dev/null +++ b/examples/morphology_test_cases/global/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/examples/morphology_test_cases/global/ctem.in b/examples/morphology_test_cases/global/ctem.in index c8dc45fb..43a9e4a1 100755 --- a/examples/morphology_test_cases/global/ctem.in +++ b/examples/morphology_test_cases/global/ctem.in @@ -67,6 +67,11 @@ Kd1 0.0001 psi 2.000 fe 4.00 ejecta_truncation 10.0 +doregotrack F dorays F superdomain F dorealistic T +quasimc T +realcraterlist craterlist.in ! list of 'real' craters for Quasi-MC runs + + diff --git a/examples/morphology_test_cases/global/ctem_driver.py b/examples/morphology_test_cases/global/ctem_driver.py index bcab88ec..897be89e 100755 --- a/examples/morphology_test_cases/global/ctem_driver.py +++ b/examples/morphology_test_cases/global/ctem_driver.py @@ -8,14 +8,18 @@ #August 2016 #Import general purpose modules + import numpy import os import subprocess import shutil +import pandas +from scipy.interpolate import interp1d #Import CTEM modules import ctem_io_readers import ctem_io_writers +import craterproduction #craterproduction had to be cp'd to example dir #Create and initialize data dictionaries for parameters and options from CTEM.in notset = '-NOTSET-' @@ -49,7 +53,9 @@ 'datfile': 'ctem.dat', 'impfile': notset, 'sfdcompare': notset, - 'sfdfile': notset} + 'sfdfile': notset, + 'quasimc': notset, + 'realcraterlist': notset} #Read ctem.in to initialize parameter values based on user input ctem_io_readers.read_ctemin(parameters,notset) @@ -71,6 +77,51 @@ impfile = parameters['workingdir'] + parameters['impfile'] prodfunction = ctem_io_readers.read_formatted_ascii(impfile, skip_lines = 0) +if (parameters['quasimc'] == 'T'): + + #Read list of real craters + print("quasi-MC mode is ON") + craterlistfile = parameters['workingdir'] + parameters['realcraterlist'] + rclist = ctem_io_readers.read_formatted_ascii(craterlistfile, skip_lines = 0) + + #Interpolate craterscale.dat to get impactor sizes from crater sizes given + df = pandas.read_csv('craterscale.dat', sep='\s+') + df['log(Dc)'] = numpy.log(df['Dcrat(m)']) + df['log(Di)'] = numpy.log(df['#Dimp(m)']) + xnew = df['log(Dc)'].values + ynew = df['log(Di)'].values + interp = interp1d(xnew, ynew, fill_value='extrapolate') + rclist[:,0] = numpy.exp(interp(numpy.log(rclist[:,0]))) + + #Convert latitude and longitude to y- and x-offset + for lat in range(0, len(rclist[:,3])): + if numpy.abs(rclist[lat,3]) > 90.0: + print("non-physical latitude on line %i of craterlist.in. Please enter a value between -90 and 90 degrees." %(lat+1)) + quit() + else: + rclist[lat,3] = rclist[lat,3] * ((parameters['pix'] * (parameters['gridsize']/2)) / 90.0) + + for lon in range(0, len(rclist[:,4])): + if numpy.abs(rclist[lon,4]) > 360.0: + print("Non-physical longitude on line %i of craterlist.in. Please enter a value between -360 and 360 degrees." %(lon+1)) + quit() + else: + if rclist[lon,4] < -180.0: + rclist[lon,4] = 360.0 - numpy.abs(rclist[lon,4]) + elif rclist[lon,4] > 180.0: + rclist[lon,4] = -(360.0 - rclist[lon,4]) + else: + rclist[lon,4] = rclist[lon,4] + + rclist[:,4] = rclist[:,4] * ((parameters['pix'] * (parameters['gridsize']/2)) / 180.0) + + #Convert age in Ga to "interval time" + rclist[:,5] = (parameters['interval'] * parameters['numintervals']) - craterproduction.Tscale(rclist[:,5], 'NPF_Moon') + rclist = rclist[rclist[:,5].argsort()] + + #Export to dat file for Fortran use + ctem_io_writers.write_realcraters(parameters, rclist) + #Create impactor production population area = (parameters['gridsize'] * parameters['pix'])**2 production = numpy.copy(prodfunction) diff --git a/examples/morphology_test_cases/global/ctem_io_readers.py b/examples/morphology_test_cases/global/ctem_io_readers.py index a58ba9ab..c825fdc9 100644 --- a/examples/morphology_test_cases/global/ctem_io_readers.py +++ b/examples/morphology_test_cases/global/ctem_io_readers.py @@ -42,6 +42,8 @@ def read_ctemin(parameters,notset): if ('savetruelist' == fields[0].lower()): parameters['savetruelist']=fields[1] if ('runtype' == fields[0].lower()): parameters['runtype']=fields[1] if ('restart' == fields[0].lower()): parameters['restart']=fields[1] + if ('quasimc' == fields[0].lower()): parameters['quasimc']=fields[1] + if ('realcraterlist' == fields[0].lower()): parameters['realcraterlist']=fields[1] if ('shadedminh' == fields[0].lower()): parameters['shadedminh'] = float(fields[1]) parameters['shadedminhdefault'] = 0 diff --git a/examples/morphology_test_cases/global/ctem_io_writers.py b/examples/morphology_test_cases/global/ctem_io_writers.py index 5d591d9a..e7a2192c 100644 --- a/examples/morphology_test_cases/global/ctem_io_writers.py +++ b/examples/morphology_test_cases/global/ctem_io_writers.py @@ -271,4 +271,12 @@ def create_rplot(parameters,odist,pdist,tdist,ph1): matplotlib.pyplot.savefig(filename) + return + +def write_realcraters(parameters, realcraters): + """writes file of real craters for use in quasi-MC runs""" + + filename = parameters['workingdir'] + 'craterlist.dat' + numpy.savetxt(filename, realcraters, fmt='%1.8e', delimiter='\t') + return \ No newline at end of file diff --git a/examples/morphology_test_cases/global/runjob b/examples/morphology_test_cases/global/runjob deleted file mode 100755 index 121d4be4..00000000 --- a/examples/morphology_test_cases/global/runjob +++ /dev/null @@ -1,14 +0,0 @@ -#!/bin/sh -l -#PBS -q daminton -#PBS -l walltime=48:00:00 -#PBS -l nodes=1:ppn=6,naccesspolicy=singleuser -#PBS -N globaltest -#PBS -j oe -source /etc/profile -module load intel -module load idl -cd $PBS_O_WORKDIR -export OMP_NUM_THREADS=6 -date > time.out -idl < start.in > term.out -date >> time.out diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index 3d005073..68dfd2bd 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -126,7 +126,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r ! Executable code - call regolith_melt_zone(user,crater,crater%imp,crater%impvel,rm,dm) + if (user%doregotrack) call regolith_melt_zone(user,crater,crater%imp,crater%impvel,rm,dm) crater%vdepth = crater%ejrim + crater%floordepth crater%vrim = crater%ejrim + crater%rimheight diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index d3a4cca2..42948217 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -211,7 +211,6 @@ module module_globals real(DP) :: ejecta_truncation ! Set the number of crater diameters to truncate the ejecta logical :: dorays ! Set T to use ray model logical :: superdomain ! Set T to include the superdomain - logical :: discontinuous ! Regolith tracking variables logical :: doregotrack ! Set T to use the regolith tracking model (EXPERIMENTAL) diff --git a/src/io/io_input.f90 b/src/io/io_input.f90 index 394f39d4..57995ed7 100644 --- a/src/io/io_input.f90 +++ b/src/io/io_input.f90 @@ -93,7 +93,6 @@ subroutine io_input(infile,user) user%fe = 5.0_DP user%dorealistic = .false. user%doquasimc = .false. - user%discontinuous = .true. user%ejecta_truncation = 10.0_DP open(unit=LUN,file=infile,status="old",iostat=ierr) @@ -331,11 +330,6 @@ subroutine io_input(infile,user) call io_get_token(line, ilength, ifirst, ilast, ierr) token = line(ifirst:ilast) read(token, *) user%ejecta_truncation - case ("DISCONTINUOUS") - ifirst = ilast + 1 - call io_get_token(line, ilength, ifirst, ilast, ierr) - token = line(ifirst:ilast) - read(token, *) user%discontinuous ! Porosity model case ("POROSITYFLG") ifirst = ilast + 1 diff --git a/test/ctem.in b/test/ctem.in index 813fb8a1..51a2fb60 100644 --- a/test/ctem.in +++ b/test/ctem.in @@ -70,7 +70,6 @@ doregotrack T ! Regolith layers track flag. Se soften_factor 1.6 soften_slope 1.8 soften_size 200.0 -discontinuous T ! Testing input. These are used to perform non-Monte Carlo tests. testflag F ! Set to T to create a single crater with user-defined impactor properties