From 46ece71b72ec8c3ed353f8349d4122c2cae5d259 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 9 Mar 2022 17:15:39 -0500 Subject: [PATCH] Updated quasimc with new temp file system --- examples/quasimc-global/ctem.in | 19 +++++------ examples/quasimc-global/misc/.gitignore | 1 - examples/quasimc-global/surf/.gitignore | 3 -- python/ctem/ctem/driver.py | 45 +++++++++---------------- python/ctem/ctem/util.py | 5 +-- src/globals/module_globals.f90 | 7 ++-- src/io/io_input.f90 | 11 +++--- src/main/CTEM.f90 | 4 +-- 8 files changed, 36 insertions(+), 59 deletions(-) mode change 100755 => 100644 examples/quasimc-global/ctem.in delete mode 100644 examples/quasimc-global/misc/.gitignore delete mode 100644 examples/quasimc-global/surf/.gitignore diff --git a/examples/quasimc-global/ctem.in b/examples/quasimc-global/ctem.in old mode 100755 new mode 100644 index ca50a7fa..c5585032 --- a/examples/quasimc-global/ctem.in +++ b/examples/quasimc-global/ctem.in @@ -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 @@ -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) diff --git a/examples/quasimc-global/misc/.gitignore b/examples/quasimc-global/misc/.gitignore deleted file mode 100644 index 773a6df9..00000000 --- a/examples/quasimc-global/misc/.gitignore +++ /dev/null @@ -1 +0,0 @@ -*.dat diff --git a/examples/quasimc-global/surf/.gitignore b/examples/quasimc-global/surf/.gitignore deleted file mode 100644 index a17efa9a..00000000 --- a/examples/quasimc-global/surf/.gitignore +++ /dev/null @@ -1,3 +0,0 @@ -*.jpg -*.png - diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index 7371770f..52174cc8 100644 --- a/python/ctem/ctem/driver.py +++ b/python/ctem/ctem/driver.py @@ -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 = { @@ -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) @@ -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+') @@ -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 @@ -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() @@ -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 diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py index b20bc08a..5d2c1700 100644 --- a/python/ctem/ctem/util.py +++ b/python/ctem/ctem/util.py @@ -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): @@ -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 diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 6157dc9d..f29b7600 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -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' @@ -284,7 +284,6 @@ 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' @@ -292,7 +291,9 @@ module module_globals 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... diff --git a/src/io/io_input.f90 b/src/io/io_input.f90 index 57995ed7..9c56298e 100644 --- a/src/io/io_input.f90 +++ b/src/io/io_input.f90 @@ -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 ", & @@ -48,7 +48,6 @@ subroutine io_input(infile,user) "TRHO_B ", & "MAT ", & "PRHO ", & - "SFDFILE ", & "VELFILE "/) integer(I4B), parameter :: seismic_numrequired=5 @@ -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 @@ -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) diff --git a/src/main/CTEM.f90 b/src/main/CTEM.f90 index db31e955..da0ea979 100644 --- a/src/main/CTEM.f90 +++ b/src/main/CTEM.f90 @@ -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)