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/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..50c84751f 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -46,9 +46,16 @@ "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"] +int_varnames = ["id", "ntp", "npl", "nplm", "discard_body_id", "collision_id"] def bool2yesno(boolval): """ @@ -106,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") @@ -146,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: @@ -157,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]) @@ -847,54 +840,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 +884,7 @@ def string_converter(da): """ if da.dtype == np.dtype(object): da = da.astype(' 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" ) @@ -409,7 +449,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" )