Skip to content

Commit

Permalink
Updated quasimc with new temp file system
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Mar 9, 2022
1 parent 36b0bf2 commit 46ece71
Show file tree
Hide file tree
Showing 8 changed files with 36 additions and 59 deletions.
19 changes: 9 additions & 10 deletions examples/quasimc-global/ctem.in
100755 → 100644
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,35 +2,35 @@


! 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
testimp 15000.0 ! Diameter of test impactor (m)
testflag F ! Set to T to create a single crater with user-defined impactor properties
testimp 15000.0 ! Diameter of test impactor (m)
testvel 18.3e3 ! Velocity of test crater (m/s)
testang 45.0 ! Impact angle of test crater (deg) - Default 90.0
testxoffset 0.0 ! x-axis offset of crater center from grid center (m) - Default 0.0
testyoffset 0.0 ! y-axis offset of crater center from grid center (m) - Default 0.0
tallyonly F ! Tally the craters without generating any craters
testtally F
quasimc T ! MC run constrained by non-MC 'real' craters given in a list
quasimc T ! MC run constrained by non-MC 'real' craters given in a list
realcraterlist craterlist.in ! list of 'real' craters for Quasi-MC runs



! IDL driver in uts
interval 6.280
numintervals 100 ! Total number of intervals (total time = interval * numintervals) <--when runtype is 'single'
interval 6.280
numintervals 100 ! Total number of interval 1 !s (total time = interval 1 ! * numintervals 1 !) <--when runtype is 'single'
restart F ! Restart a previous run
impfile NPFextrap.dat ! Impactor SFD rate file (col 1: Dimp (m), col 2: ! impactors > D (m**(-2) y**(-1))
popupconsole F ! Pop up console window every output interval
popupconsole F ! Pop up console window every output interval 1 !
saveshaded F ! Output shaded relief images
saverego F ! Output regolith map images
savepres F ! Output simplified console display images (presentation-compatible images)
savetruelist T ! Save the true cumulative crater distribution for each interval (large file size)
savetruelist T ! Save the true cumulative crater distribution for each interval 1 ! (large file size)
sfdcompare LOLASethCraterCatalogv8gt20-binned.dat ! File name for the SFD comparison file used in the console display
shadedminh -85.0 ! Minimum height for shaded relief map (m) (Default - automatic)
shadedmaxh 85.0 ! Maximum height for shaded relief map (m) (Default - automatic)
runtype single ! Run type: options are normal / statistical
! single: craters accumulate in successive intervals
! statistical: surface is reset between intervals
! single: craters accumulate in successive interval 1 !s
! statistical: surface is reset between interval 1 !s

! CTEM required inputs
seed 76535 ! Random number generator seed
Expand All @@ -52,7 +52,6 @@ ybar_r 0.00e6 ! Regolith strength (Pa)
gaccel 1.62e0 ! Gravitational acceleration at target (m/s**2)
trad 1737.35e3 ! Target radius (m)
prho 2500.0e0 ! Projectile density (kg/m**3)
sfdfile production.dat ! Impactor SFD file (col 1: Dimp (m), col 2: ! impactors > D
velfile lunar-MBA-impactor-velocities.dat ! Impactor velocity dist file

! Seismic shaking input (required if seismic shaking is set to T, otherwise ignored)
Expand Down
1 change: 0 additions & 1 deletion examples/quasimc-global/misc/.gitignore
View file

This file was deleted.

3 changes: 0 additions & 3 deletions examples/quasimc-global/surf/.gitignore
View file

This file was deleted.

45 changes: 15 additions & 30 deletions python/ctem/ctem/driver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,16 +40,16 @@ def __init__(self, param_file="ctem.in", isnew=True):
'shadedminh': 0.0,
'shadedmaxh': 0.0,
'workingdir': currentdir,
'ctemfile': os.path.join(currentdir, param_file),
'ctemfile': param_file,
'impfile': None,
'sfdcompare': None,
'sfdfile': None,
'quasimc': None,
'sfdfile' : None,
'realcraterlist': None,
}

self.user = util.read_user_input(self.user)

# This is containsall files generated by the main Fortran CTEM program plus the ctem.dat file that
# communicates the state of the simulation from the Python driver to the Fortran program
self.output_filenames = {
Expand All @@ -71,8 +71,11 @@ def __init__(self, param_file="ctem.in", isnew=True):
'ejmin' : 'ejecta_table_min.dat',
'testprof' : 'testprofile.dat',
'craterscale' : 'craterscale.dat',
'craterlist' : 'craterlist.dat'
'craterlist' : 'craterlist.dat',
'sfdfile' : 'production.dat'
}
if self.user['sfdfile'] is not None: # Override the default sfdfile name if the user supplies an alternative
self.output_filenames['sfdfile'] = self.user['sfdfile']

for k, v in self.output_filenames.items():
self.output_filenames[k] = os.path.join(currentdir, v)
Expand Down Expand Up @@ -120,29 +123,18 @@ def __init__(self, param_file="ctem.in", isnew=True):

#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)
self.tempfile = os.path.join(currentdir, 'temp.in')
tempfile = os.path.join(currentdir, 'temp.in')

# Generate craterlist.dat
with open(self.user['ctemfile']) as perm, open(self.tempfile, "w") as temp:
for line in perm:
temp.write(line)
shutil.copy2(self.user['ctemfile'], tempfile )

#Write a temporary input file to generate the necessary quasimc files
self.permfile = self.user['ctemfile']
self.temp = util.write_temp_input(self.tempfile)
self.user['ctemfile'] = self.tempfile
self.user = util.read_user_input(self.user)
self.permin = os.path.join(currentdir, 'perm.in')
self.ctemin = os.path.join(currentdir, 'ctem.in')
os.rename(self.permfile, self.permin)
os.rename(self.tempfile, self.ctemin)
util.write_temp_input(tempfile)
util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
self.compute_one_interval()

os.replace(self.permin, self.ctemin)
self.user['ctemfile'] = self.permfile
self.user = util.read_user_input(self.user)
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+')
Expand Down Expand Up @@ -183,7 +175,7 @@ def scale_production(self):
self.production[:, 1] = self.production[:, 1] * area * self.user['interval']

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

return

Expand All @@ -202,7 +194,7 @@ def run(self):
self.redirect_outputs(['dat'], 'misc')

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

# Read in output files
self.read_output()
Expand Down Expand Up @@ -352,13 +344,6 @@ def cleanup(self):
os.remove(filename)
except OSError as error:
print(error)
# Delete scaled production file
prodfile = self.user['sfdfile']
print(f"Deleting file {prodfile}")
try:
os.remove(prodfile)
except OSError as error:
print(error)

return

Expand Down
5 changes: 1 addition & 4 deletions python/ctem/ctem/util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -243,8 +243,6 @@ def read_user_input(user):
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['sfdfile'] is None):
print('Invalid value for or missing variable SFDFILE in ' + inputfile)
if (user['impfile'] is None):
print('Invalid value for or missing variable IMPFILE in ' + inputfile)
if (user['popupconsole'] is None):
Expand Down Expand Up @@ -336,8 +334,7 @@ def write_datfile(user, filename, seedarr):
return


def write_production(user, production):
filename = user['sfdfile']
def write_production(filename, production):
np.savetxt(filename, production, fmt='%1.8e', delimiter=' ')

return
Expand Down
7 changes: 4 additions & 3 deletions src/globals/module_globals.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -270,7 +270,7 @@ module module_globals
integer(I4B) :: pbarpos
character(len=PBARSIZE) :: pbarchar

! Grid array file names
! Default file names
character(*),parameter :: DIAMFILE = 'surface_diam.dat'
character(*),parameter :: TIMEFILE = 'surface_time.dat'
character(*),parameter :: EJCOVFILE = 'surface_ejc.dat'
Expand All @@ -284,15 +284,16 @@ module module_globals
character(*),parameter :: POROFILE = 'porosity_porosity.dat'
character(*),parameter :: DEPTHFILE = 'porosity_depth.dat'
character(*),parameter :: POROTHICKFILE = 'porosity_thickness.dat'
!character(*),parameter :: THICKFILE = 'surface_crustal_thickness.dat'
character(*),parameter :: POSFILE = 'surface_pos.dat'
character(*),parameter :: TDISTFILE = 'tdistribution.dat'
character(*),parameter :: TLISTFILE = 'tcumulative.dat'
character(*),parameter :: ODISTFILE = 'odistribution.dat'
character(*),parameter :: OLISTFILE = 'ocumulative.dat'
character(*),parameter :: PDISTFILE = 'pdistribution.dat'
character(*),parameter :: CRTSCLFILE = 'craterscale.dat'
character(*),parameter :: DATFILE = 'ctem.dat'
character(*),parameter :: SFDFILE = 'production.dat'
character(*),parameter :: USERFILE = 'ctem.in'
character(*),parameter :: DATFILE = 'ctem.dat' !
character(*),parameter :: MASSFILE = 'impactmass.dat'
character(*),parameter :: RCFILE = 'craterlist.dat' !not sure if this is where this line should go, but putting it here for now...

Expand Down
11 changes: 5 additions & 6 deletions src/io/io_input.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ subroutine io_input(infile,user)
integer(I4B), parameter :: LUN = 7
integer(I4B) :: ierr, ilength, ifirst, ilast,i
character(STRMAX) :: line, token
integer(I4B), parameter :: numrequired=18
integer(I4B), parameter :: numrequired=17
character(STRMAX),dimension(numrequired),parameter :: requiredvar = (/"GRIDSIZE ",&
"NUMLAYERS", &
"PIX ", &
Expand All @@ -48,7 +48,6 @@ subroutine io_input(infile,user)
"TRHO_B ", &
"MAT ", &
"PRHO ", &
"SFDFILE ", &
"VELFILE "/)

integer(I4B), parameter :: seismic_numrequired=5
Expand Down Expand Up @@ -94,6 +93,7 @@ subroutine io_input(infile,user)
user%dorealistic = .false.
user%doquasimc = .false.
user%ejecta_truncation = 10.0_DP
write(user%sfdfile,*) trim(adjustl(SFDFILE))

open(unit=LUN,file=infile,status="old",iostat=ierr)
if (ierr /= 0) then
Expand Down Expand Up @@ -212,14 +212,13 @@ subroutine io_input(infile,user)
ifirst = ilast + 1
call io_get_token(line, ilength, ifirst, ilast, ierr)
token = line(ifirst:ilast)
read(token, *) user%sfdfile
read(token, *) user%velfile
ismissing(17)=.false.
case (trim(requiredvar(18)))
case ("SFDFILE")
ifirst = ilast + 1
call io_get_token(line, ilength, ifirst, ilast, ierr)
token = line(ifirst:ilast)
read(token,*) user%velfile
ismissing(18)=.false.
read(token, *) user%sfdfile
case ("DEPLIMIT")
ifirst = ilast + 1
call io_get_token(line, ilength, ifirst, ilast, ierr)
Expand Down
4 changes: 2 additions & 2 deletions src/main/CTEM.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,9 +69,9 @@ program CTEM
call io_splash()
narg = command_argument_count()
if (narg /= 0) then
call get_command_argument(1, infile, status = ierr)
call get_command_argument(1, infile, status = ierr) ! Use first argument as the user input file name
else
infile="ctem.in"
infile = USERFILE ! No arguments, so use the default file name for the user inputs
end if
write(*,*) 'Reading input file ',trim(adjustl(infile))
call io_input(infile,user)
Expand Down

0 comments on commit 46ece71

Please sign in to comment.