From 82cbe5357990818a4d26f934bbed767212c04c86 Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Wed, 9 Nov 2022 12:54:53 -0500 Subject: [PATCH 1/5] reordered calls to nf90_inq_varid in netcdf.f90 --- src/netcdf/netcdf.f90 | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 0cd27f737..0c444669a 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -366,12 +366,10 @@ module subroutine netcdf_open(self, param, readonly) call check( nf90_inq_varid(self%ncid, TIME_DIMNAME, self%time_varid), "netcdf_open nf90_inq_varid time_varid" ) call check( nf90_inq_varid(self%ncid, ID_DIMNAME, self%id_varid), "netcdf_open nf90_inq_varid id_varid" ) - call check( nf90_inq_varid(self%ncid, NPL_VARNAME, self%npl_varid), "netcdf_open nf90_inq_varid npl_varid" ) - call check( nf90_inq_varid(self%ncid, NTP_VARNAME, self%ntp_varid), "netcdf_open nf90_inq_varid ntp_varid" ) - if (param%integrator == SYMBA) call check( nf90_inq_varid(self%ncid, NPLM_VARNAME, self%nplm_varid), "netcdf_open nf90_inq_varid nplm_varid" ) call check( nf90_inq_varid(self%ncid, NAME_VARNAME, self%name_varid), "netcdf_open nf90_inq_varid name_varid" ) call check( nf90_inq_varid(self%ncid, PTYPE_VARNAME, self%ptype_varid), "netcdf_open nf90_inq_varid ptype_varid" ) call check( nf90_inq_varid(self%ncid, STATUS_VARNAME, self%status_varid), "netcdf_open nf90_inq_varid status_varid" ) + call check( nf90_inq_varid(self%ncid, GMASS_VARNAME, self%Gmass_varid), "netcdf_open nf90_inq_varid Gmass_varid" ) if ((param%out_form == XV) .or. (param%out_form == XVEL)) then call check( nf90_inq_varid(self%ncid, XHX_VARNAME, self%xhx_varid), "netcdf_open nf90_inq_varid xhx_varid" ) @@ -409,7 +407,6 @@ module subroutine netcdf_open(self, param, readonly) call check( nf90_inq_varid(self%ncid, OMEGA_VARNAME, self%omega_varid), "netcdf_open nf90_inq_varid omega_varid" ) call check( nf90_inq_varid(self%ncid, CAPM_VARNAME, self%capm_varid), "netcdf_open nf90_inq_varid capm_varid" ) end if - call check( nf90_inq_varid(self%ncid, GMASS_VARNAME, self%Gmass_varid), "netcdf_open nf90_inq_varid Gmass_varid" ) if (param%lrhill_present) call check( nf90_inq_varid(self%ncid, RHILL_VARNAME, self%rhill_varid), "netcdf_open nf90_inq_varid rhill_varid" ) From 516ae33973372d55b760c21edb4ea1020155b6ec Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Wed, 9 Nov 2022 12:55:29 -0500 Subject: [PATCH 2/5] commented out readonly bits so that new variables can be added --- src/netcdf/netcdf.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 0c444669a..9f2e9fba0 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -353,9 +353,9 @@ module subroutine netcdf_open(self, param, readonly) character(len=NF90_MAX_NAME) :: str_dim_name mode = NF90_WRITE - if (present(readonly)) then - if (readonly) mode = NF90_NOWRITE - end if + !if (present(readonly)) then + ! if (readonly) mode = NF90_NOWRITE + !end if call check( nf90_open(param%outfile, mode, self%ncid), "netcdf_open nf90_open" ) From 0a859905fd364825e9b7dddc2df9be6b3f4f9306 Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Wed, 9 Nov 2022 12:56:27 -0500 Subject: [PATCH 3/5] check if npl, ntp, or nplm exist and if not, calculate and define them --- src/netcdf/netcdf.f90 | 42 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 9f2e9fba0..e66498c35 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -351,6 +351,9 @@ module subroutine netcdf_open(self, param, readonly) ! Internals integer(I4B) :: mode, status character(len=NF90_MAX_NAME) :: str_dim_name + integer(I4B) :: idmax + real(DP), dimension(:), allocatable :: gmtemp + logical, dimension(:), allocatable :: tpmask, plmask, plmmask mode = NF90_WRITE !if (present(readonly)) then @@ -371,6 +374,45 @@ module subroutine netcdf_open(self, param, readonly) call check( nf90_inq_varid(self%ncid, STATUS_VARNAME, self%status_varid), "netcdf_open nf90_inq_varid status_varid" ) call check( nf90_inq_varid(self%ncid, GMASS_VARNAME, self%Gmass_varid), "netcdf_open nf90_inq_varid Gmass_varid" ) + if ((nf90_inq_varid(self%ncid, NPL_VARNAME, self%npl_varid) /= nf90_noerr) .or. & + (nf90_inq_varid(self%ncid, NTP_VARNAME, self%ntp_varid) /= nf90_noerr) .or. & + ((nf90_inq_varid(self%ncid, NPLM_VARNAME, self%nplm_varid) /= nf90_noerr) .and. (param%integrator == SYMBA))) then + call check( nf90_inquire_dimension(self%ncid, self%id_dimid, len=idmax), "netcdf_open nf90_inquire_dimension id_dimid" ) + allocate(gmtemp(idmax)) + call check( nf90_get_var(self%ncid, self%Gmass_varid, gmtemp, start=[1,1]), "netcdf_open nf90_getvar Gmass_varid" ) + allocate(tpmask(idmax)) + allocate(plmask(idmax)) + allocate(plmmask(idmax)) + plmask(:) = gmtemp(:) == gmtemp(:) + tpmask(:) = .not. plmask(:) + plmask(1) = .false. ! This is the central body + select type (param) + class is (symba_parameters) + plmmask(:) = gmtemp(:) > param%GMTINY .and. plmask(:) + end select + if ((nf90_inq_varid(self%ncid, NPL_VARNAME, self%npl_varid) /= nf90_noerr)) then + call check( nf90_redef(self%ncid), "netcdf_open nf90_redef npl_varid") + call check( nf90_def_var(self%ncid, NPL_VARNAME, NF90_INT, self%time_dimid, self%npl_varid), "netcdf_open nf90_def_var npl_varid" ) + call check( nf90_enddef(self%ncid), "netcdf_open nf90_enddef npl_varid") + call check( nf90_put_var(self%ncid, self%npl_varid, count(plmask(:)), start=[1]), "netcdf_open nf90_put_var npl_varid" ) + call check( nf90_inq_varid(self%ncid, NPL_VARNAME, self%npl_varid), "netcdf_open nf90_inq_varid npl_varid" ) + end if + if (nf90_inq_varid(self%ncid, NTP_VARNAME, self%ntp_varid) /= nf90_noerr) then + call check( nf90_redef(self%ncid), "netcdf_open nf90_redef ntp_varid") + call check( nf90_def_var(self%ncid, NTP_VARNAME, NF90_INT, self%time_dimid, self%ntp_varid), "netcdf_open nf90_def_var ntp_varid" ) + call check( nf90_enddef(self%ncid), "netcdf_open nf90_enddef ntp_varid") + call check( nf90_put_var(self%ncid, self%ntp_varid, count(tpmask(:)), start=[1]), "netcdf_open nf90_put_var ntp_varid" ) + call check( nf90_inq_varid(self%ncid, NTP_VARNAME, self%ntp_varid), "netcdf_open nf90_inq_varid ntp_varid" ) + end if + if ((nf90_inq_varid(self%ncid, NPLM_VARNAME, self%nplm_varid) /= nf90_noerr) .and. (param%integrator == SYMBA)) then + call check( nf90_redef(self%ncid), "netcdf_open nf90_redef nplm_varid") + call check( nf90_def_var(self%ncid, NPLM_VARNAME, NF90_INT, self%time_dimid, self%nplm_varid), "netcdf_open nf90_def_var nplm_varid" ) + call check( nf90_enddef(self%ncid), "netcdf_open nf90_enddef nplm_varid") + call check( nf90_put_var(self%ncid, self%nplm_varid, count(plmmask(:)), start=[1]), "netcdf_open nf90_put_var nplm_varid" ) + call check( nf90_inq_varid(self%ncid, NPLM_VARNAME, self%nplm_varid), "netcdf_open nf90_inq_varid nplm_varid" ) + end if + end if + if ((param%out_form == XV) .or. (param%out_form == XVEL)) then call check( nf90_inq_varid(self%ncid, XHX_VARNAME, self%xhx_varid), "netcdf_open nf90_inq_varid xhx_varid" ) call check( nf90_inq_varid(self%ncid, XHY_VARNAME, self%xhy_varid), "netcdf_open nf90_inq_varid xhy_varid" ) From d592f0e6a1d8f5ae30823276d9408beaf463ab79 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 9 Nov 2022 17:08:44 -0500 Subject: [PATCH 4/5] Made improvements to type handling by the Swiftest xarray methods --- python/swiftest/swiftest/init_cond.py | 1 + python/swiftest/swiftest/io.py | 80 ++++++++------------ python/swiftest/swiftest/simulation_class.py | 5 ++ 3 files changed, 38 insertions(+), 48 deletions(-) diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index 0166e78a8..bf2720a3c 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/python/swiftest/swiftest/init_cond.py @@ -417,4 +417,5 @@ def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, ds_float = da_float.to_dataset(dim="vec") ds_str = da_str.to_dataset(dim="vec") ds = xr.combine_by_coords([ds_float, ds_str,info_ds]) + return ds \ No newline at end of file diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 28303f8f8..dc184b175 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -49,6 +49,7 @@ # This defines Xarray Dataset variables that are strings, which must be processed due to quirks in how NetCDF-Fortran # handles strings differently than Python's Xarray. string_varnames = ["name", "particle_type", "status", "origin_type"] +int_varnames = ["id", "ntp", "npl", "nplm", "discard_body_id", "collision_id"] def bool2yesno(boolval): """ @@ -847,54 +848,14 @@ def swiftest2xr(param, verbose=True): ------- xarray dataset """ - if ((param['OUT_TYPE'] == 'REAL8') or (param['OUT_TYPE'] == 'REAL4')): - dims = ['time', 'id', 'vec'] - cb = [] - pl = [] - tp = [] - cbn = None - try: - with FortranFile(param['BIN_OUT'], 'r') as f: - for t, cbid, cbnames, cvec, clab, \ - npl, plid, plnames, pvec, plab, \ - ntp, tpid, tpnames, tvec, tlab in swiftest_stream(f, param): - # Prepare frames by adding an extra axis for the time coordinate - cbframe = np.expand_dims(cvec, axis=0) - plframe = np.expand_dims(pvec, axis=0) - tpframe = np.expand_dims(tvec, axis=0) - - - # Create xarray DataArrays out of each body type - cbxr = xr.DataArray(cbframe, dims=dims, coords={'time': t, 'id': cbid, 'vec': clab}) - cbxr = cbxr.assign_coords(name=("id", cbnames)) - plxr = xr.DataArray(plframe, dims=dims, coords={'time': t, 'id': plid, 'vec': plab}) - plxr = plxr.assign_coords(name=("id", plnames)) - tpxr = xr.DataArray(tpframe, dims=dims, coords={'time': t, 'id': tpid, 'vec': tlab}) - tpxr = tpxr.assign_coords(name=("id", tpnames)) - - cb.append(cbxr) - pl.append(plxr) - tp.append(tpxr) - - sys.stdout.write('\r' + f"Reading in time {t[0]:.3e}") - sys.stdout.flush() - except IOError: - print(f"Error encountered reading in {param['BIN_OUT']}") - - cbda = xr.concat(cb, dim='time') - plda = xr.concat(pl, dim='time') - tpda = xr.concat(tp, dim='time') - - cbds = cbda.to_dataset(dim='vec') - plds = plda.to_dataset(dim='vec') - tpds = tpda.to_dataset(dim='vec') - if verbose: print('\nCreating Dataset') - ds = xr.combine_by_coords([cbds, plds, tpds]) - elif ((param['OUT_TYPE'] == 'NETCDF_DOUBLE') or (param['OUT_TYPE'] == 'NETCDF_FLOAT')): + if ((param['OUT_TYPE'] == 'NETCDF_DOUBLE') or (param['OUT_TYPE'] == 'NETCDF_FLOAT')): if verbose: print('\nCreating Dataset from NetCDF file') ds = xr.open_dataset(param['BIN_OUT'], mask_and_scale=False) - ds = clean_string_values(ds) + if param['OUT_TYPE'] == "NETCDF_DOUBLE": + ds = fix_types(ds,ftype=np.float64) + elif param['OUT_TYPE'] == "NETCDF_FLOAT": + ds = fix_types(ds,ftype=np.float32) else: print(f"Error encountered. OUT_TYPE {param['OUT_TYPE']} not recognized.") return None @@ -931,7 +892,7 @@ def string_converter(da): """ if da.dtype == np.dtype(object): da = da.astype(' Date: Wed, 9 Nov 2022 17:57:07 -0500 Subject: [PATCH 5/5] Restructured when the param file is read in and fixed typo in bool2str --- examples/Basic_Simulation/param.in | 46 ----------------- python/swiftest/swiftest/io.py | 52 +++++++++----------- python/swiftest/swiftest/simulation_class.py | 21 ++++---- 3 files changed, 34 insertions(+), 85 deletions(-) delete mode 100644 examples/Basic_Simulation/param.in diff --git a/examples/Basic_Simulation/param.in b/examples/Basic_Simulation/param.in deleted file mode 100644 index bec3573de..000000000 --- a/examples/Basic_Simulation/param.in +++ /dev/null @@ -1,46 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -! VERSION Swiftest parameter input -T0 0.0 -TSTOP 10 -DT 0.005 -ISTEP_OUT 200 -ISTEP_DUMP 200 -OUT_FORM XVEL -OUT_TYPE NETCDF_DOUBLE -OUT_STAT REPLACE -IN_TYPE ASCII -PL_IN pl.in -TP_IN tp.in -CB_IN cb.in -BIN_OUT out.nc -CHK_QMIN 0.004650467260962157 -CHK_RMIN 0.004650467260962157 -CHK_RMAX 1000.0 -CHK_EJECT 1000.0 -CHK_QMIN_COORD HELIO -CHK_QMIN_RANGE 0.004650467260962157 1000.0 -MU2KG 1.988409870698051e+30 -TU2S 31557600.0 -DU2M 149597870700.0 -IN_FORM EL -EXTRA_FORCE NO -BIG_DISCARD NO -CHK_CLOSE YES -RHILL_PRESENT YES -FRAGMENTATION YES -ROTATION YES -TIDES NO -ENERGY YES -GR YES -INTERACTION_LOOPS ADAPTIVE -ENCOUNTER_CHECK ADAPTIVE -GMTINY 1e-06 -MIN_GMFRAG 1e-09 diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index dc184b175..50c84751f 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -46,6 +46,12 @@ "YARKOVSKY", "YORP"] +int_param = ["ISTEP_OUT", "ISTEP_DUMP"] +float_param = ["T0", "TSTART", "TSTOP", "DT", "CHK_RMIN", "CHK_RMAX", "CHK_EJECT", "CHK_QMIN", "DU2M", "MU2KG", + "TU2S", "MIN_GMFRAG", "GMTINY"] + +upper_str_param = ["OUT_TYPE","OUT_FORM","OUT_STAT","IN_TYPE","IN_FORM"] + # This defines Xarray Dataset variables that are strings, which must be processed due to quirks in how NetCDF-Fortran # handles strings differently than Python's Xarray. string_varnames = ["name", "particle_type", "status", "origin_type"] @@ -107,10 +113,10 @@ def str2bool(input_str): valid_false = ["NO", "N", "F", "FALSE", ".FALSE."] if input_str.upper() in valid_true: return True - elif input_str.lower() in valid_false: + elif input_str.upper() in valid_false: return False else: - raise ValueError(f"{input_str} cannot is not recognized as boolean") + raise ValueError(f"{input_str} is not recognized as boolean") @@ -147,7 +153,8 @@ def read_swiftest_param(param_file_name, param, verbose=True): A dictionary containing the entries in the user parameter file """ param['! VERSION'] = f"Swiftest parameter input from file {param_file_name}" - + + # Read param.in file if verbose: print(f'Reading Swiftest file {param_file_name}') try: @@ -158,38 +165,23 @@ def read_swiftest_param(param_file_name, param, verbose=True): if fields[0][0] != '!': key = fields[0].upper() param[key] = fields[1] - #for key in param: - # if (key == fields[0].upper()): param[key] = fields[1] # Special case of CHK_QMIN_RANGE requires a second input if fields[0].upper() == 'CHK_QMIN_RANGE': alo = real2float(fields[1]) ahi = real2float(fields[2]) param['CHK_QMIN_RANGE'] = f"{alo} {ahi}" - - param['ISTEP_OUT'] = int(param['ISTEP_OUT']) - param['ISTEP_DUMP'] = int(param['ISTEP_DUMP']) - param['OUT_TYPE'] = param['OUT_TYPE'].upper() - param['OUT_FORM'] = param['OUT_FORM'].upper() - param['OUT_STAT'] = param['OUT_STAT'].upper() - param['IN_TYPE'] = param['IN_TYPE'].upper() - param['IN_FORM'] = param['IN_FORM'].upper() - param['T0'] = real2float(param['T0']) - param['TSTART'] = real2float(param['TSTART']) - param['TSTOP'] = real2float(param['TSTOP']) - param['DT'] = real2float(param['DT']) - param['CHK_RMIN'] = real2float(param['CHK_RMIN']) - param['CHK_RMAX'] = real2float(param['CHK_RMAX']) - param['CHK_EJECT'] = real2float(param['CHK_EJECT']) - param['CHK_QMIN'] = real2float(param['CHK_QMIN']) - param['DU2M'] = real2float(param['DU2M']) - param['MU2KG'] = real2float(param['MU2KG']) - param['TU2S'] = real2float(param['TU2S']) - param['INTERACTION_LOOPS'] = param['INTERACTION_LOOPS'].upper() - param['ENCOUNTER_CHECK'] = param['ENCOUNTER_CHECK'].upper() - if 'GMTINY' in param: - param['GMTINY'] = real2float(param['GMTINY']) - if 'MIN_GMFRAG' in param: - param['MIN_GMFRAG'] = real2float(param['MIN_GMFRAG']) + + for uc in upper_str_param: + if uc in param: + param[uc] = param[uc].upper() + + for i in int_param: + if i in param and type(i) != int: + param[i] = int(param[i]) + + for f in float_param: + if f in param and type(f) is str: + param[f] = real2float(param[f]) for b in bool_param: if b in param: param[b] = str2bool(param[b]) diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index a7f7f7a6d..a7b8ed27c 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -240,6 +240,15 @@ def __init__(self, self.verbose = verbose self.restart = restart + # If the parameter file is in a different location than the current working directory, we will need + # to use it to properly open bin files + self.sim_dir = os.path.dirname(os.path.realpath(param_file)) + if read_param: + if os.path.exists(param_file): + self.read_param(param_file, codename=codename, verbose=self.verbose) + else: + print(f"{param_file} not found.") + self.set_distance_range(rmin=rmin, rmax=rmax, verbose = False) self.set_unit_system(MU=MU, DU=DU, TU=TU, @@ -279,14 +288,7 @@ def __init__(self, verbose = False ) - # If the parameter file is in a different location than the current working directory, we will need - # to use it to properly open bin files - self.sim_dir = os.path.dirname(os.path.realpath(param_file)) - if read_param: - if os.path.exists(param_file): - self.read_param(param_file, codename=codename, verbose=self.verbose) - else: - print(f"{param_file} not found.") + if read_old_output_file: binpath = os.path.join(self.sim_dir,self.param['BIN_OUT']) @@ -1487,7 +1489,8 @@ def read_param(self, param_file, codename="Swiftest", verbose=True): self.ds : xarray dataset """ if codename == "Swiftest": - self.param = io.read_swiftest_param(param_file, self.param, verbose=verbose) + param_old = self.param.copy() + self.param = io.read_swiftest_param(param_file, param_old, verbose=verbose) self.codename = "Swiftest" elif codename == "Swifter": self.param = io.read_swifter_param(param_file, verbose=verbose)