diff --git a/swiftest/simulation_class.py b/swiftest/simulation_class.py index d6627c325..a4a2b12b8 100644 --- a/swiftest/simulation_class.py +++ b/swiftest/simulation_class.py @@ -413,7 +413,7 @@ def run(self,dask: bool = False, **kwargs): self.set_parameter(**kwargs) # Write out the current parameter set before executing run - self.write_param() + self.save(verbose=False) if self.codename != "Swiftest": warnings.warn(f"Running an integration is not yet supported for {self.codename}",stacklevel=2) diff --git a/tests/test_swiftest.py b/tests/test_swiftest.py index e6d1cf4c9..28f4e4a83 100755 --- a/tests/test_swiftest.py +++ b/tests/test_swiftest.py @@ -259,123 +259,4 @@ def test05_gr(self): if __name__ == '__main__': - unittest.main() - if os.path.exists("simdir"): - shutil.rmtree("simdir") - - # Add the modern planets and the Sun using the JPL Horizons Database. - # Add the modern planets and the Sun using the JPL Horizons Database. - sim.add_solar_system_body(major_bodies) - - # Display the run configuration parameters. - param = sim.get_parameter(verbose=False) - sim.save() - - for f in file_list: - self.assertTrue(os.path.exists(f)) - - print("\ntest_read_ic: Test whether we can read back initial conditions files created by test_gen_ic") - sim = swiftest.Simulation(read_param=True) - # Check if all names in Dataset read in from file match the expected list of names - self.assertTrue((major_bodies == sim.init_cond['name']).all(), msg="Name mismatch in Dataset") - - # Check to see if all parameter values read in from file match the expected parameters saved when generating the file - self.assertTrue(all([v == param[k] for k,v in sim.param.items() if k in param])) - - def test02_integrators(self): - """ - Tests that Swiftest is able to integrate a collection of massive bodies and test particles with all available integrators - """ - print("\ntest_integrators: Tests that Swiftest is able to integrate a collection of massive bodies and test particles with all available integrators") - sim = swiftest.Simulation() - sim.clean() - - # Add the modern planets and the Sun using the JPL Horizons Database. - # Add the modern planets and the Sun using the JPL Horizons Database. - sim.add_solar_system_body(major_bodies) - - # Display the run configuration parameters. - param = sim.get_parameter(verbose=False) - sim.save() - - # Add 10 user-defined test particles. - ntp = 10 - - name_tp = [f"TestParticle_{i:02}" for i in range(1,ntp+1)] - a_tp = rng.uniform(0.3, 1.5, ntp) - e_tp = rng.uniform(0.0, 0.2, ntp) - inc_tp = rng.uniform(0.0, 10, ntp) - capom_tp = rng.uniform(0.0, 360.0, ntp) - omega_tp = rng.uniform(0.0, 360.0, ntp) - capm_tp = rng.uniform(0.0, 360.0, ntp) - - integrators= ["helio","whm","rmvs","symba"] - sim.add_body(name=name_tp, a=a_tp, e=e_tp, inc=inc_tp, capom=capom_tp, omega=omega_tp, capm=capm_tp) - sim.set_parameter(tstart=0.0, tstop=0.10, dt=0.01, istep_out=1, dump_cadence=0) - for i in integrators: - sim.run(integrator=i) - - - def test03_conservation(self): - """ - Tests that Swiftest conserves mass, energy, and momentum to within acceptable tolerances. - """ - print("\ntest_conservation: Tests that Swiftest conserves mass, energy, and momentum to within acceptable tolerances.") - - # Error limits - L_slope_limit = 1e-10 - E_slope_limit = 1e-8 - GM_limit = 1e-14 - - sim = swiftest.Simulation() - sim.clean() - - sim.add_solar_system_body(major_bodies) - - dt = 0.01 - nout = 1000 - tstop = 1e4 - tstep_out = tstop / nout - - sim.run(tstart=0.0, tstop=tstop, dt=dt, tstep_out=tstep_out, dump_cadence=0, compute_conservation_values=True, integrator="symba") - - def fit_func(x,slope,b): - """ - Linear function for curve fitting - """ - return slope * x + b - - # Calculate the angular momentum error - sim.data['L_tot'] = sim.data['L_orbit'] + sim.data['L_spin'] + sim.data['L_escape'] - sim.data['DL'] = sim.data['L_tot'] - sim.data['L_tot'].isel(time=0) - L_error = swiftest.tool.magnitude(sim.data,'DL') / swiftest.tool.magnitude(sim.data.isel(time=0), 'L_tot') - - # Calculate the energy error - E_error = (sim.data['TE'] - sim.data['TE'].isel(time=0)) / sim.data['TE'].isel(time=0) - - # Calculate the mass error - sim.data['GMtot'] = sim.data['Gmass'].sum(dim='name',skipna=True) + sim.data['GMescape'] - GM_error = (sim.data['GMtot'] - sim.data['GMtot'].isel(time=0)) / sim.data['GMtot'].isel(time=0) - GM_final = GM_error.isel(time=-1).values - - # Compute the slope of the error vs time fit - E_fit_result = E_error.curvefit("time",fit_func) - L_fit_result = L_error.curvefit("time",fit_func) - - E_slope = E_fit_result['curvefit_coefficients'].sel(param='slope').values - L_slope = L_fit_result['curvefit_coefficients'].sel(param='slope').values - - # Print the final errors - print("\n") - print(f"Angular momentum error slope: {L_slope:.2e}/{sim.TU_name}") - print(f"Energy error slope: {E_slope:.2e}/{sim.TU_name}") - print(f"Final Mass Error: {GM_final:.2e}") - - self.assertLess(np.abs(L_slope),L_slope_limit, msg=f"Angular Momentum Error of {L_slope:.2e}/{sim.TU_name} higher than threshold value of {L_slope_limit:.2e}/{sim.TU_name}") - self.assertLess(np.abs(E_slope),E_slope_limit, msg=f"Energy Error of {E_slope:.2e}/{sim.TU_name} higher than threshold value of {E_slope_limit:.2e}/{sim.TU_name}") - self.assertLess(np.abs(GM_final),GM_limit, msg=f"Mass Error of {GM_final:.2e} higher than threshold value of {GM_limit:.2e}") - -if __name__ == '__main__': - unittest.main() - if os.path.exists("simdir"): - shutil.rmtree("simdir") + unittest.main() \ No newline at end of file