From 8a00b7ab0870aa4c700c33defad402195b82c1ed Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 18 Oct 2022 13:07:21 -0400 Subject: [PATCH 1/5] Added a new example that adds Neptune and Pluto to a Helio run. I had to do a lot of maintenence on the initial conditions generators in the Python side of things in order to get it to generate properly --- examples/helio_neptune_pluto/init_cond.py | 42 +++ python/swiftest/swiftest/init_cond.py | 337 +++++++++---------- python/swiftest/swiftest/io.py | 5 +- python/swiftest/swiftest/simulation_class.py | 33 +- src/modules/swiftest_classes.f90 | 37 ++ 5 files changed, 266 insertions(+), 188 deletions(-) create mode 100755 examples/helio_neptune_pluto/init_cond.py diff --git a/examples/helio_neptune_pluto/init_cond.py b/examples/helio_neptune_pluto/init_cond.py new file mode 100755 index 000000000..c3c6c96c7 --- /dev/null +++ b/examples/helio_neptune_pluto/init_cond.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 +import swiftest + +sim = swiftest.Simulation() +sim.param['IN_TYPE'] = "NETCDF_DOUBLE" +sim.param['NC_IN'] = "init_cond.nc" +sim.param['IN_FORM'] = "XV" +sim.param['BIN_OUT'] = "bin.nc" +sim.param['OUT_FORM'] = "XVEL" +sim.param['OUT_TYPE'] = "NETCDF_DOUBLE" + +sim.param['MU2KG'] = swiftest.MSun +sim.param['TU2S'] = swiftest.YR2S +sim.param['DU2M'] = swiftest.AU2M +sim.param['T0'] = 0.0 +sim.param['DT'] = 5.0 +sim.param['TSTOP'] = 1.e5 +sim.param['ISTEP_OUT'] = 100 +sim.param['ISTEP_DUMP'] = 100 +sim.param['CHK_QMIN_COORD'] = "HELIO" +sim.param['CHK_QMIN'] = swiftest.RSun / swiftest.AU2M +sim.param['CHK_QMIN_RANGE'] = f"{swiftest.RSun / swiftest.AU2M} 1000.0" +sim.param['CHK_RMIN'] = swiftest.RSun / swiftest.AU2M +sim.param['CHK_RMAX'] = 1000.0 +sim.param['CHK_EJECT'] = 1000.0 +sim.param['OUT_STAT'] = "REPLACE" +sim.param['RHILL_PRESENT'] = "NO" +sim.param['GR'] = 'NO' +sim.param['CHK_CLOSE'] = "NO" + +bodyid = { + "Sun": 0, + "Neptune": 8, + "Pluto" : 9 +} + +for name, id in bodyid.items(): + sim.add(name, idval=id, date="2027-04-30") + +sim.save("param.in") + + diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index 89e1c8dee..932c6c960 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/python/swiftest/swiftest/init_cond.py @@ -40,6 +40,7 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date, ds): if plname in planetid: ispl = True + idval = planetid[plname] else: ispl = False print(f"\nMassive body {plname} not found or not yet supported") @@ -117,144 +118,129 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date, ds): solarrot = planetrot['Sun'] * param['TU2S'] rotcb = solarpole.cartesian * solarrot Ipsun = np.array([0.0, 0.0, planetIpz['Sun']]) - - cbid = np.array([0]) - cvec = np.vstack([GMcb, Rcb, J2RP2, J4RP4]) - if param['ROTATION'] == 'YES': - cvec = np.vstack([cvec, Ipsun[0], Ipsun[1], Ipsun[2], rotcb.x, rotcb.y, rotcb.z]) - - # Horizons date time internal variables - tstart = datetime.date.fromisoformat(ephemerides_start_date) - tstep = datetime.timedelta(days=1) - tend = tstart + tstep - ephemerides_end_date = tend.isoformat() - ephemerides_step = '1d' - + param_tmp = param param_tmp['OUT_FORM'] = 'XVEL' - clab, plab, tlab = swiftest.io.make_swiftest_labels(param) - dims = ['time', 'id', 'vec'] - infodims = ['id', 'vec'] - t = np.array([0.0]) - - key = plname - if ispl: - val = planetid[key] - else: - val = -1 - if key == "Sun" : # Create central body + if plname == "Sun" : # Create central body print("Creating the Sun as a central body") - cbframe = np.expand_dims(cvec.T, axis=0) - cbxr = xr.DataArray(cbframe, dims=dims, coords={'time': t, 'id': cbid, 'vec': clab}) - cbxr = cbxr.assign_coords(name=('id', [key])) - cbds = cbxr.to_dataset(dim='vec') - ds = xr.combine_by_coords([ds, cbds]) + v1 = None + v2 = None + v3 = None + v4 = None + v5 = None + v6 = None + rhill = None + GMpl = GMcb + Rpl = Rcb + J2 = J2RP2 + J4 = J4RP4 + if param['ROTATION'] == 'YES': + Ip1 = [Ipsun[0]] + Ip2 = [Ipsun[1]] + Ip3 = [Ipsun[2]] + rotx = [rotcb.x] + roty = [rotcb.y] + rotz = [rotcb.z] + else: + Ip1 = None + Ip2 = None + Ip3 = None + rotx = None + roty = None + rotz = None else: # Fetch solar system ephemerides from Horizons - print(f"Fetching ephemerides data for {key} from JPL/Horizons") + print(f"Fetching ephemerides data for {plname} from JPL/Horizons") - p1 = [] - p2 = [] - p3 = [] - p4 = [] - p5 = [] - p6 = [] - p7 = [] - p8 = [] - p9 = [] - p10 = [] - p11 = [] - p12 = [] - rhill = [] - Rpl = [] - GMpl = [] - Ip1 = [] - Ip2 = [] - Ip3 = [] - rotx = [] - roty = [] - rotz = [] + # Horizons date time internal variables + tstart = datetime.date.fromisoformat(ephemerides_start_date) + tstep = datetime.timedelta(days=1) + tend = tstart + tstep + ephemerides_end_date = tend.isoformat() + ephemerides_step = '1d' + + v1 = [] + v2 = [] + v3 = [] + v4 = [] + v5 = [] + v6 = [] + J2 = None + J4 = None pldata = {} - if ispl: - pldata[key] = Horizons(id=val, id_type='majorbody', location='@sun', - epochs={'start': ephemerides_start_date, 'stop': ephemerides_end_date, - 'step': ephemerides_step}) - else: - pldata[key] = Horizons(id=key, id_type='smallbody', location='@sun', - epochs={'start': ephemerides_start_date, 'stop': ephemerides_end_date, - 'step': ephemerides_step}) + pldata[plname] = Horizons(id=idval, location='@sun', + epochs={'start': ephemerides_start_date, 'stop': ephemerides_end_date, + 'step': ephemerides_step}) - if (param['OUT_FORM'] == 'XV' or param['OUT_FORM'] == 'XVEL'): - p1.append(pldata[key].vectors()['x'][0] * DCONV) - p2.append(pldata[key].vectors()['y'][0] * DCONV) - p3.append(pldata[key].vectors()['z'][0] * DCONV) - p4.append(pldata[key].vectors()['vx'][0] * VCONV) - p5.append(pldata[key].vectors()['vy'][0] * VCONV) - p6.append(pldata[key].vectors()['vz'][0] * VCONV) - p7.append(pldata[key].elements()['a'][0] * DCONV) - p8.append(pldata[key].elements()['e'][0]) - p9.append(pldata[key].elements()['incl'][0]) - p10.append(pldata[key].elements()['Omega'][0]) - p11.append(pldata[key].elements()['w'][0]) - p12.append(pldata[key].elements()['M'][0]) - elif param['OUT_FORM'] == 'EL': - p1.append(pldata[key].elements()['a'][0] * DCONV) - p2.append(pldata[key].elements()['e'][0]) - p3.append(pldata[key].elements()['incl'][0]) - p4.append(pldata[key].elements()['Omega'][0]) - p5.append(pldata[key].elements()['w'][0]) - p6.append(pldata[key].elements()['M'][0]) - p7.append(pldata[key].vectors()['x'][0] * DCONV) - p8.append(pldata[key].vectors()['y'][0] * DCONV) - p9.append(pldata[key].vectors()['z'][0] * DCONV) - p10.append(pldata[key].vectors()['vx'][0] * VCONV) - p11.append(pldata[key].vectors()['vy'][0] * VCONV) - p12.append(pldata[key].vectors()['vz'][0] * VCONV) - pvec = np.vstack([p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12]) + if param['IN_FORM'] == 'XV': + v1.append(pldata[plname].vectors()['x'][0] * DCONV) + v2.append(pldata[plname].vectors()['y'][0] * DCONV) + v3.append(pldata[plname].vectors()['z'][0] * DCONV) + v4.append(pldata[plname].vectors()['vx'][0] * VCONV) + v5.append(pldata[plname].vectors()['vy'][0] * VCONV) + v6.append(pldata[plname].vectors()['vz'][0] * VCONV) + elif param['IN_FORM'] == 'EL': + v1.append(pldata[plname].elements()['a'][0] * DCONV) + v2.append(pldata[plname].elements()['e'][0]) + v3.append(pldata[plname].elements()['incl'][0]) + v4.append(pldata[plname].elements()['Omega'][0]) + v5.append(pldata[plname].elements()['w'][0]) + v6.append(pldata[plname].elements()['M'][0]) if ispl: - Rpl.append(planetradius[key] * DCONV) - GMpl.append(GMcb[0] / MSun_over_Mpl[key]) - pvec = np.vstack([pvec, GMpl, Rpl]) + GMpl = [] + GMpl.append(GMcb[0] / MSun_over_Mpl[plname]) + if param['CHK_CLOSE'] == 'YES': + Rpl = [] + Rpl.append(planetradius[plname] * DCONV) + else: + Rpl = None # Generate planet value vectors if (param['RHILL_PRESENT'] == 'YES'): - rhill.append(pldata[key].elements()['a'][0] * DCONV * (3 * MSun_over_Mpl[key]) ** (-THIRDLONG)) - pvec = np.vstack([pvec, rhill]) + rhill = [] + rhill.append(pldata[plname].elements()['a'][0] * DCONV * (3 * MSun_over_Mpl[plname]) ** (-THIRDLONG)) + else: + rhill = None if (param['ROTATION'] == 'YES'): - RA = pldata[key].ephemerides()['NPole_RA'][0] - DEC = pldata[key].ephemerides()['NPole_DEC'][0] + Ip1 = [] + Ip2 = [] + Ip3 = [] + rotx = [] + roty = [] + rotz = [] + RA = pldata[plname].ephemerides()['NPole_RA'][0] + DEC = pldata[plname].ephemerides()['NPole_DEC'][0] rotpole = SkyCoord(ra=RA * u.degree, dec=DEC * u.degree) - rotrate = planetrot[key] * param['TU2S'] + rotrate = planetrot[plname] * param['TU2S'] rot = rotpole.cartesian * rotrate - Ip = np.array([0.0, 0.0, planetIpz[key]]) + Ip = np.array([0.0, 0.0, planetIpz[plname]]) Ip1.append(Ip[0]) Ip2.append(Ip[1]) Ip3.append(Ip[2]) rotx.append(rot.x) roty.append(rot.y) rotz.append(rot.z) - pvec = np.vstack([pvec, Ip1, Ip2, Ip3, rotx, roty, rotz]) - else: - plab = tlab.copy() - - if idval is None: - plid = np.array([planetid[key]], dtype=int) + else: + Ip1 = None + Ip2 = None + Ip3 = None + rotx = None + roty = None + rotz = None else: - plid = np.array([idval], dtype=int) + GMpl = None - # Prepare frames by adding an extra axis for the time coordinate - plframe = np.expand_dims(pvec.T, axis=0) - plxr = xr.DataArray(plframe, dims=dims, coords={'time': t, 'id': plid, 'vec': plab}) - plxr = plxr.assign_coords(name=('id', [plname])) - plds = plxr.to_dataset(dim='vec') - ds = xr.combine_by_coords([ds, plds]) + if idval is None: + plid = np.array([planetid[plname]], dtype=int) + else: + plid = np.array([idval], dtype=int) - return ds + return plid,[plname],v1,v2,v3,v4,v5,v6,GMpl,Rpl,rhill,Ip1,Ip2,Ip3,rotx,roty,rotz,J2,J4 -def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, rhill=None, Ip1=None, Ip2=None, Ip3=None, rotx=None, roty=None, rotz=None, t=0.0): +def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, rhill=None, Ip1=None, Ip2=None, Ip3=None, rotx=None, roty=None, rotz=None, J2=None, J4=None,t=0.0): """ Converts and stores the variables of all bodies in an xarray dataset. @@ -302,6 +288,11 @@ def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, ------- ds : xarray dataset """ + if v1 is None: # This is the central body + iscb = True + else: + iscb = False + if param['ROTATION'] == 'YES': if Ip1 is None: Ip1 = np.full_like(v1, 0.4) @@ -318,12 +309,12 @@ def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, dims = ['time', 'id', 'vec'] infodims = ['id', 'vec'] - if GMpl is not None: + if not iscb and GMpl is not None: ispl = True else: ispl = False - if ispl and Rpl is None: + if ispl and param['CHK_CLOSE'] == 'YES' and Rpl is None: print("Massive bodies need a radius value.") return None if ispl and rhill is None and param['RHILL_PRESENT'] == 'YES': @@ -335,76 +326,84 @@ def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, param['OUT_FORM'] = param['IN_FORM'] clab, plab, tlab, infolab_float, infolab_int, infolab_str = swiftest.io.make_swiftest_labels(param) param['OUT_FORM'] = old_out_form - vec_float = np.vstack([v1, v2, v3, v4, v5, v6]) vec_str = np.vstack([namevals]) label_str = ["name"] - if ispl: - label_float = plab.copy() - vec_float = np.vstack([vec_float, GMpl, Rpl]) - if param['RHILL_PRESENT'] == 'YES': - vec_float = np.vstack([vec_float, rhill]) + if iscb: + label_float = clab.copy() + vec_float = np.vstack([GMpl,Rpl,J2,J4]) if param['ROTATION'] == 'YES': vec_float = np.vstack([vec_float, Ip1, Ip2, Ip3, rotx, roty, rotz]) + particle_type = "Central Body" else: - label_float = tlab.copy() - if param['OUT_STAT'] == 'APPEND': # Include particle information in initial conditions when appending an existing run + vec_float = np.vstack([v1, v2, v3, v4, v5, v6]) if ispl: + label_float = plab.copy() + vec_float = np.vstack([vec_float, GMpl]) + if param['CHK_CLOSE'] == 'YES': + vec_float = np.vstack([vec_float, Rpl]) + if param['RHILL_PRESENT'] == 'YES': + vec_float = np.vstack([vec_float, rhill]) + if param['ROTATION'] == 'YES': + vec_float = np.vstack([vec_float, Ip1, Ip2, Ip3, rotx, roty, rotz]) particle_type = np.repeat("Massive Body",idvals.size) else: + label_float = tlab.copy() particle_type = np.repeat("Test Particle",idvals.size) - origin_type = np.repeat("User Added Body",idvals.size) - origin_time = np.full_like(v1,t) - collision_id = np.full_like(idvals,0) - origin_xhx = v1 - origin_xhy = v2 - origin_xhz = v3 - origin_vhx = v4 - origin_vhy = v5 - origin_vhz = v6 - discard_time = np.full_like(v1,-1.0) - status = np.repeat("ACTIVE",idvals.size) - discard_xhx = np.zeros_like(v1) - discard_xhy = np.zeros_like(v1) - discard_xhz = np.zeros_like(v1) - discard_vhx = np.zeros_like(v1) - discard_vhy = np.zeros_like(v1) - discard_vhz = np.zeros_like(v1) - discard_body_id = np.full_like(idvals,-1) - info_vec_float = np.vstack([ - origin_time, - origin_xhx, - origin_xhy, - origin_xhz, - origin_vhx, - origin_vhy, - origin_vhz, - discard_time, - discard_xhx, - discard_xhy, - discard_xhz, - discard_vhx, - discard_vhy, - discard_vhz]) - info_vec_int = np.vstack([collision_id, discard_body_id]) - info_vec_str = np.vstack([particle_type, origin_type, status]) - frame_float = info_vec_float.T - frame_int = info_vec_int.T - frame_str = info_vec_str.T - da_float = xr.DataArray(frame_float, dims=infodims, coords={'id': idvals, 'vec': infolab_float}) - da_int = xr.DataArray(frame_int, dims=infodims, coords={'id': idvals, 'vec': infolab_int}) - da_str = xr.DataArray(frame_str, dims=infodims, coords={'id': idvals, 'vec': infolab_str}) - ds_float = da_float.to_dataset(dim="vec") - ds_int = da_int.to_dataset(dim="vec") - ds_str = da_str.to_dataset(dim="vec") - info_ds = xr.combine_by_coords([ds_float, ds_int, ds_str]) + origin_type = np.repeat("User Added Body",idvals.size) + origin_time = np.full_like(v1,t) + collision_id = np.full_like(idvals,0) + origin_xhx = v1 + origin_xhy = v2 + origin_xhz = v3 + origin_vhx = v4 + origin_vhy = v5 + origin_vhz = v6 + discard_time = np.full_like(v1,-1.0) + status = np.repeat("ACTIVE",idvals.size) + discard_xhx = np.zeros_like(v1) + discard_xhy = np.zeros_like(v1) + discard_xhz = np.zeros_like(v1) + discard_vhx = np.zeros_like(v1) + discard_vhy = np.zeros_like(v1) + discard_vhz = np.zeros_like(v1) + discard_body_id = np.full_like(idvals,-1) + info_vec_float = np.vstack([ + origin_time, + origin_xhx, + origin_xhy, + origin_xhz, + origin_vhx, + origin_vhy, + origin_vhz, + discard_time, + discard_xhx, + discard_xhy, + discard_xhz, + discard_vhx, + discard_vhy, + discard_vhz]) + info_vec_int = np.vstack([collision_id, discard_body_id]) + info_vec_str = np.vstack([particle_type, origin_type, status]) + frame_float = info_vec_float.T + frame_int = info_vec_int.T + frame_str = info_vec_str.T + if param['IN_TYPE'] == 'NETCDF_FLOAT': + ftype=np.float32 + elif param['IN_TYPE'] == 'NETCDF_DOUBLE': + ftype=np.float64 + da_float = xr.DataArray(frame_float, dims=infodims, coords={'id': idvals, 'vec': infolab_float}).astype(ftype) + da_int = xr.DataArray(frame_int, dims=infodims, coords={'id': idvals, 'vec': infolab_int}) + da_str = xr.DataArray(frame_str, dims=infodims, coords={'id': idvals, 'vec': infolab_str}) + ds_float = da_float.to_dataset(dim="vec") + ds_int = da_int.to_dataset(dim="vec") + ds_str = da_str.to_dataset(dim="vec") + info_ds = xr.combine_by_coords([ds_float, ds_int, ds_str]) frame_float = np.expand_dims(vec_float.T, axis=0) frame_str = vec_str.T - da_float = xr.DataArray(frame_float, dims=dims, coords={'time': [t], 'id': idvals, 'vec': label_float}) + da_float = xr.DataArray(frame_float, dims=dims, coords={'time': [t], 'id': idvals, 'vec': label_float}).astype(ftype) da_str= xr.DataArray(frame_str, dims=infodims, coords={'id': idvals, 'vec': label_str}) 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]) - if param['OUT_STAT'] == 'APPEND': - ds = xr.combine_by_coords([ds,info_ds]) + 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 916c9a5f0..8d4d56167 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -643,7 +643,7 @@ def swiftest_stream(f, param): t11 = f.read_reals(np.float64) t12 = f.read_reals(np.float64) - clab, plab, tlab = make_swiftest_labels(param) + clab, plab, tlab, infolab_float, infolab_int, infolab_str = make_swiftest_labels(param) if npl > 0: pvec = np.vstack([p1, p2, p3, p4, p5, p6]) @@ -967,7 +967,7 @@ def select_active_from_frame(ds, param, framenum=-1): iactive = iframe['id'].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['xhx'])), drop=True).id else: iactive = iframe['id'].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['a'])), drop = True).id - frame = frame.isel(id=iactive.values) + frame = frame.sel(id=iactive.values) return frame @@ -995,6 +995,7 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram # Note: xarray will call the character array dimension string32. The Fortran code # will rename this after reading frame = unclean_string_values(frame) + print(f"Writing initial conditions to file {infile_name}") frame.to_netcdf(path=infile_name) return frame diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index a62723bbf..218ca553e 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -19,11 +19,9 @@ def __init__(self, codename="Swiftest", param_file="param.in", readbin=True, ver 'T0': "0.0", 'TSTOP': "0.0", 'DT': "0.0", - 'PL_IN': "pl.in", - 'TP_IN': "tp.in", - 'CB_IN': "cb.in", - 'IN_TYPE': "ASCII", - 'IN_FORM': "EL", + 'IN_FORM': "XV", + 'IN_TYPE': "NETCDF_DOUBLE", + 'NC_IN' : "init_cond.nc", 'ISTEP_OUT': "1", 'ISTEP_DUMP': "1", 'BIN_OUT': "bin.nc", @@ -51,8 +49,8 @@ def __init__(self, codename="Swiftest", param_file="param.in", readbin=True, ver 'TIDES': "NO", 'ENERGY': "NO", 'GR': "YES", - 'INTERACTION_LOOPS': "ADAPTIVE", - 'ENCOUNTER_CHECK': "ADAPTIVE" + 'INTERACTION_LOOPS': "TRIANGULAR", + 'ENCOUNTER_CHECK': "TRIANGULAR" } self.codename = codename self.verbose = verbose @@ -82,11 +80,12 @@ def add(self, plname, date=date.today().isoformat(), idval=None): ------- self.ds : xarray dataset """ - self.ds = init_cond.solar_system_horizons(plname, idval, self.param, date, self.ds) + #self.ds = init_cond.solar_system_horizons(plname, idval, self.param, date, self.ds) + self.addp(*init_cond.solar_system_horizons(plname, idval, self.param, date, self.ds)) return - def addp(self, idvals, namevals, t1, t2, t3, t4, t5, t6, GMpl=None, Rpl=None, rhill=None, Ip1=None, Ip2=None, Ip3=None, rotx=None, roty=None, rotz=None, t=None): + def addp(self, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None, rhill=None, Ip1=None, Ip2=None, Ip3=None, rotx=None, roty=None, rotz=None, J2=None,J4=None,t=None): """ Adds a body (test particle or massive body) to the internal DataSet given a set up 6 vectors (orbital elements or cartesian state vectors, depending on the value of self.param). Input all angles in degress @@ -95,17 +94,17 @@ def addp(self, idvals, namevals, t1, t2, t3, t4, t5, t6, GMpl=None, Rpl=None, rh ---------- idvals : integer Array of body index values. - t1 : float + v1 : float xh for param['IN_FORM'] == "XV"; a for param['IN_FORM'] == "EL" - t2 : float + v2 : float yh for param['IN_FORM'] == "XV"; e for param['IN_FORM'] == "EL" - t3 : float + v3 : float zh for param['IN_FORM'] == "XV"; inc for param['IN_FORM'] == "EL" - t4 : float + v4 : float vhxh for param['IN_FORM'] == "XV"; capom for param['IN_FORM'] == "EL" - t5 : float + v5 : float vhyh for param['IN_FORM'] == "XV"; omega for param['IN_FORM'] == "EL" - t6 : float + v6 : float vhzh for param['IN_FORM'] == "XV"; capm for param['IN_FORM'] == "EL" Gmass : float Optional: Array of G*mass values if these are massive bodies @@ -126,7 +125,7 @@ def addp(self, idvals, namevals, t1, t2, t3, t4, t5, t6, GMpl=None, Rpl=None, rh if t is None: t = self.param['T0'] - dsnew = init_cond.vec2xr(self.param, idvals, namevals, t1, t2, t3, t4, t5, t6, GMpl, Rpl, rhill, Ip1, Ip2, Ip3, rotx, roty, rotz, t) + dsnew = init_cond.vec2xr(self.param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl, Rpl, rhill, Ip1, Ip2, Ip3, rotx, roty, rotz, J2, J4, t) if dsnew is not None: self.ds = xr.combine_by_coords([self.ds, dsnew]) self.ds['ntp'] = self.ds['id'].where(np.isnan(self.ds['Gmass'])).count(dim="id") @@ -325,7 +324,7 @@ def save(self, param_file, framenum=-1, codename="Swiftest"): """ if codename == "Swiftest": - io.swiftest_xr2infile(self.ds, self.param, framenum) + io.swiftest_xr2infile(ds=self.ds, param=self.param, framenum=framenum,infile_name=self.param['NC_IN']) self.write_param(param_file) elif codename == "Swifter": if self.codename == "Swiftest": diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index ca2d97347..a86dc2312 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -1595,6 +1595,40 @@ module function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1) end function util_solve_rkf45 end interface + interface + module recursive pure subroutine qsort_DP(arr, ind) + implicit none + ! Arguments + real(DP), dimension(:), intent(inout) :: arr + integer(I4B),dimension(:),intent(out), optional :: ind + end subroutine qsort_DP + + module recursive pure subroutine qsort_I4B(arr, ind) + implicit none + integer(I4B), dimension(:), intent(inout) :: arr + integer(I4B), dimension(:), intent(out), optional :: ind + end subroutine qsort_I4B + + module recursive pure subroutine qsort_I4B_I8Bind(arr, ind) + implicit none + integer(I4B), dimension(:), intent(inout) :: arr + integer(I8B), dimension(:), intent(out), optional :: ind + end subroutine qsort_I4B_I8Bind + + module recursive pure subroutine qsort_I8B_I8Bind(arr, ind) + implicit none + integer(I8B), dimension(:), intent(inout) :: arr + integer(I8B), dimension(:), intent(out), optional :: ind + end subroutine qsort_I8B_I8Bind + + module recursive pure subroutine qsort_SP(arr, ind) + implicit none + real(SP), dimension(:), intent(inout) :: arr + integer(I4B),dimension(:),intent(out), optional :: ind + end subroutine qsort_SP + + end interface + interface util_sort module pure subroutine util_sort_i4b(arr) implicit none @@ -1739,8 +1773,11 @@ module subroutine util_sort_tp(self, sortby, ascending) character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order end subroutine util_sort_tp + end interface + + interface util_spill module subroutine util_spill_arr_char_string(keeps, discards, lspill_list, ldestructive) implicit none From 659cb10db5b0bd1d1caaa4681188e34fed8b6b2d Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 18 Oct 2022 13:38:01 -0400 Subject: [PATCH 2/5] Fixed names of J2 and J4 variables --- python/swiftest/swiftest/io.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 8d4d56167..0263f4e0a 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -477,7 +477,7 @@ def make_swiftest_labels(param): plab.append('radius') if param['RHILL_PRESENT'] == 'YES': plab.append('rhill') - clab = ['Gmass', 'radius', 'J_2', 'J_4'] + clab = ['Gmass', 'radius', 'j2rp2', 'j4rp4'] if param['ROTATION'] == 'YES': clab.append('Ip1') clab.append('Ip2') @@ -1002,16 +1002,16 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram # All other file types need seperate files for each of the inputs cb = frame.where(frame.id == 0, drop=True) pl = frame.where(frame.id > 0, drop=True) - pl = pl.where(np.invert(np.isnan(pl['Gmass'])), drop=True).drop_vars(['J_2', 'J_4'],errors="ignore") - tp = frame.where(np.isnan(frame['Gmass']), drop=True).drop_vars(['Gmass', 'radius', 'J_2', 'J_4'],errors="ignore") + pl = pl.where(np.invert(np.isnan(pl['Gmass'])), drop=True).drop_vars(['j2rp2', 'j2rp2'],errors="ignore") + tp = frame.where(np.isnan(frame['Gmass']), drop=True).drop_vars(['Gmass', 'radius', 'j2rp2', 'j4rp4'],errors="ignore") GMSun = np.double(cb['Gmass']) if param['CHK_CLOSE'] == 'YES': RSun = np.double(cb['radius']) else: RSun = param['CHK_RMIN'] - J2 = np.double(cb['J_2']) - J4 = np.double(cb['J_4']) + J2 = np.double(cb['j2rp2']) + J4 = np.double(cb['j4rp4']) cbname = cb['name'].values[0] if param['ROTATION'] == 'YES': Ip1cb = np.double(cb['Ip1']) @@ -1187,16 +1187,16 @@ def swifter_xr2infile(ds, param, framenum=-1): cb = frame.where(frame.id == 0, drop=True) pl = frame.where(frame.id > 0, drop=True) - pl = pl.where(np.invert(np.isnan(pl['Gmass'])), drop=True).drop_vars(['J_2', 'J_4']) - tp = frame.where(np.isnan(frame['Gmass']), drop=True).drop_vars(['Gmass', 'radius', 'J_2', 'J_4']) + pl = pl.where(np.invert(np.isnan(pl['Gmass'])), drop=True).drop_vars(['j2rp2', 'j4rp4']) + tp = frame.where(np.isnan(frame['Gmass']), drop=True).drop_vars(['Gmass', 'radius', 'j2rp2', 'j4rp4']) GMSun = np.double(cb['Gmass']) if param['CHK_CLOSE'] == 'YES': RSun = np.double(cb['radius']) else: RSun = param['CHK_RMIN'] - param['J2'] = np.double(cb['J_2']) - param['J4'] = np.double(cb['J_4']) + param['J2'] = np.double(cb['j2rp2']) + param['J4'] = np.double(cb['j4rp4']) if param['IN_TYPE'] == 'ASCII': # Swiftest Central body file From eb770cd6664b383e5d9b44021321f50d75a5c298 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 19 Oct 2022 14:59:41 -0400 Subject: [PATCH 3/5] Rolled back putting the qsort functions as module functions and added a new check that turns off energy if this is a new run with a NetCDF input file --- src/io/io.f90 | 4 ++++ src/modules/swiftest_classes.f90 | 34 -------------------------------- 2 files changed, 4 insertions(+), 34 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 2ecac13b0..27068aeb9 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -1401,6 +1401,10 @@ module subroutine io_read_in_system(self, param) allocate(tmp_param, source=param) tmp_param%outfile = param%in_netcdf tmp_param%out_form = param%in_form + if (.not. param%lrestart) then + ! Turn off energy computation so we don't have to feed it into the initial conditions + tmp_param%lenergy = .false. + end if ierr = self%read_frame(tmp_param%nciu, tmp_param) deallocate(tmp_param) if (ierr /=0) call util_exit(FAILURE) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index a86dc2312..1ef589770 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -1595,40 +1595,6 @@ module function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1) end function util_solve_rkf45 end interface - interface - module recursive pure subroutine qsort_DP(arr, ind) - implicit none - ! Arguments - real(DP), dimension(:), intent(inout) :: arr - integer(I4B),dimension(:),intent(out), optional :: ind - end subroutine qsort_DP - - module recursive pure subroutine qsort_I4B(arr, ind) - implicit none - integer(I4B), dimension(:), intent(inout) :: arr - integer(I4B), dimension(:), intent(out), optional :: ind - end subroutine qsort_I4B - - module recursive pure subroutine qsort_I4B_I8Bind(arr, ind) - implicit none - integer(I4B), dimension(:), intent(inout) :: arr - integer(I8B), dimension(:), intent(out), optional :: ind - end subroutine qsort_I4B_I8Bind - - module recursive pure subroutine qsort_I8B_I8Bind(arr, ind) - implicit none - integer(I8B), dimension(:), intent(inout) :: arr - integer(I8B), dimension(:), intent(out), optional :: ind - end subroutine qsort_I8B_I8Bind - - module recursive pure subroutine qsort_SP(arr, ind) - implicit none - real(SP), dimension(:), intent(inout) :: arr - integer(I4B),dimension(:),intent(out), optional :: ind - end subroutine qsort_SP - - end interface - interface util_sort module pure subroutine util_sort_i4b(arr) implicit none From 36a69b3f024c9db1f83961c57a7877794ff59369 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 21 Oct 2022 09:44:43 -0400 Subject: [PATCH 4/5] Updated initial conditions to compute energy --- examples/helio_neptune_pluto/init_cond.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/helio_neptune_pluto/init_cond.py b/examples/helio_neptune_pluto/init_cond.py index c3c6c96c7..975003c3a 100755 --- a/examples/helio_neptune_pluto/init_cond.py +++ b/examples/helio_neptune_pluto/init_cond.py @@ -1,7 +1,7 @@ #!/usr/bin/env python3 import swiftest -sim = swiftest.Simulation() +sim = swiftest.Simulation(readbin=False) sim.param['IN_TYPE'] = "NETCDF_DOUBLE" sim.param['NC_IN'] = "init_cond.nc" sim.param['IN_FORM'] = "XV" @@ -15,8 +15,8 @@ sim.param['T0'] = 0.0 sim.param['DT'] = 5.0 sim.param['TSTOP'] = 1.e5 -sim.param['ISTEP_OUT'] = 100 -sim.param['ISTEP_DUMP'] = 100 +sim.param['ISTEP_OUT'] = 20 +sim.param['ISTEP_DUMP'] = 20 sim.param['CHK_QMIN_COORD'] = "HELIO" sim.param['CHK_QMIN'] = swiftest.RSun / swiftest.AU2M sim.param['CHK_QMIN_RANGE'] = f"{swiftest.RSun / swiftest.AU2M} 1000.0" @@ -25,6 +25,7 @@ sim.param['CHK_EJECT'] = 1000.0 sim.param['OUT_STAT'] = "REPLACE" sim.param['RHILL_PRESENT'] = "NO" +sim.param['ENERGY'] = "YES" sim.param['GR'] = 'NO' sim.param['CHK_CLOSE'] = "NO" @@ -36,7 +37,7 @@ for name, id in bodyid.items(): sim.add(name, idval=id, date="2027-04-30") - + sim.save("param.in") From 40232d24554e2782f64e398747d7f67c477bc9b7 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 21 Oct 2022 09:45:40 -0400 Subject: [PATCH 5/5] Fixed NaN issue in potential energy computation. Inconsistent array bounds on mask (this may be the core reason why fraggle fails sometimes) --- src/util/util_get_energy_momentum.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index ae59c158a..65e9b90d0 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -196,10 +196,12 @@ subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, x !$omp firstprivate(npl) & !$omp reduction(+:pe) do i = 1, npl - do concurrent(j = i+1:npl, lmask(i) .and. lmask(j)) - pepl(j) = - (Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j)) - end do - pe = pe + sum(pepl(i+1:npl), lmask(i)) + if (lmask(i)) then + do concurrent(j = i+1:npl, lmask(i) .and. lmask(j)) + pepl(j) = - (Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j)) + end do + pe = pe + sum(pepl(i+1:npl), lmask(i+1:npl)) + end if end do !$omp end parallel do pe = pe + sum(pecb(1:npl), lmask(1:npl))