Updated examples and removed the variable discontinuous from user

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

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()

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
+
+

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)

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

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

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

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)

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