diff --git a/.gitignore b/.gitignore index 5056a889a..13ec87692 100644 --- a/.gitignore +++ b/.gitignore @@ -23,6 +23,7 @@ dump* !examples/** *ipynb_checkpoints **/.DS_Store +!test/test_suite.py #Documentation !docs/* diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 585aedb3e..047ec6a1c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -170,14 +170,9 @@ ENDIF (USE_OPENMP) # Add an automatic test + ADD_CUSTOM_TARGET(Run_Test ALL - COMMAND echo "***Running Python script examples/Basic_Simulation/initial_conditions.py***" - COMMAND python ${CMAKE_SOURCE_DIR}/examples/Basic_Simulation/initial_conditions.py - COMMAND echo "***Saving data to directory examples/Basic_Simulation/simdata***" - COMMAND echo "***Running Python script examples/Basic_Simulation/output_reader.py***" - COMMAND python ${CMAKE_SOURCE_DIR}/examples/Basic_Simulation/output_reader.py - COMMAND echo "***Plotting results as examples/Basic_Simulation/output.eps***" - COMMAND echo "***Calculating errors with Python script examples/Basic_Simulation/errors.py***" - COMMAND python ${CMAKE_SOURCE_DIR}/examples/Basic_Simulation/errors.py + COMMAND echo "***Running Python script tests/test_suite.py***" + COMMAND python ${CMAKE_SOURCE_DIR}/test/test_suite.py COMMAND echo "***Test Complete***" ) diff --git a/test/test_suite.py b/test/test_suite.py new file mode 100755 index 000000000..36588314e --- /dev/null +++ b/test/test_suite.py @@ -0,0 +1,161 @@ +#!/usr/bin/env python3 + +""" + Copyright 2023 - 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. +""" + +""" +Tests that energy and momentum errors are within tolerances in a Swiftest simulation + +Input +------ + +Output +------ +None +""" + +import swiftest +import unittest +import os +import numpy as np +from numpy.random import default_rng + +rng = default_rng(seed=123) + +major_bodies = ["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"] +param = {} + + +class TestSwiftest(unittest.TestCase): + + def test01_gen_ic(self): + """ + Tests that Swiftest is able to successfully generate a set of initial conditions in a file without any exceptions being raised + """ + print("\ntest_gen_ic: Test whether we can generate simulation initial conditions test") + # Files that are expected to be generated: + simdir = "simdata" + file_list = [simdir, os.path.join(simdir,"param.in"), os.path.join(simdir,"init_cond.nc")] + + 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() + + for f in file_list: + self.assertTrue(os.path.exists(f)) + + def test02_read_ic(self): + """ + Tests that Swiftest is able to read a set of pre-existing initial conditions files and that they contain the correct data + """ + 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 test03_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(read_param=True) + + # 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= ["whm","helio","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.02, dt=0.01, istep_out=1, dump_cadence=0) + integrators= ["whm","helio","rmvs","symba"] + for i in integrators: + try: + sim.run(integrator=i) + except: + self.fail(f"Failed with integrator {i}") + + + def test04_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.") + sim = swiftest.Simulation(read_param=True) + + # Add 5 user-defined semi-interacting massive bodies. + npl = 5 + density_pl = 3000.0 / (sim.param['MU2KG'] / sim.param['DU2M'] ** 3) + + name_pl = [f"SemiBody_{i:02}" for i in range(1,npl+1)] + a_pl = rng.uniform(0.3, 1.5, npl) + e_pl = rng.uniform(0.0, 0.2, npl) + inc_pl = rng.uniform(0.0, 10, npl) + capom_pl = rng.uniform(0.0, 360.0, npl) + omega_pl = rng.uniform(0.0, 360.0, npl) + capm_pl = rng.uniform(0.0, 360.0, npl) + M_pl = np.array([6e20, 8e20, 1e21, 3e21, 5e21]) * sim.KG2MU + R_pl = np.full(npl, (3 * M_pl/ (4 * np.pi * density_pl)) ** (1.0 / 3.0)) + Ip_pl = np.full((npl,3),0.4,) + rot_pl = np.zeros((npl,3)) + mtiny = 1.1 * np.max(M_pl) + + sim.add_body(name=name_pl, a=a_pl, e=e_pl, inc=inc_pl, capom=capom_pl, omega=omega_pl, capm=capm_pl, mass=M_pl, radius=R_pl, Ip=Ip_pl, rot=rot_pl) + + sim.run(start=0.0, tstop=1.0e3, dt=0.01, istep_out=100, dump_cadence=0, compute_conservation_values=True, mtiny=mtiny, integrator="symba") + # 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) + sim.data['L_error'] = swiftest.tool.magnitude(sim.data,'DL') / swiftest.tool.magnitude(sim.data.isel(time=0), 'L_tot') + L_final = sim.data['L_error'][-1].values + + # Calculate the energy error + sim.data['E_error'] = (sim.data['TE'] - sim.data['TE'].isel(time=0)) / sim.data['TE'].isel(time=0) + E_final = sim.data['E_error'][-1].values + + # Calculate the mass error + sim.data['GMtot'] = sim.data['Gmass'].sum(dim='name',skipna=True) + sim.data['GMescape'] + sim.data['GM_error'] = (sim.data['GMtot'] - sim.data['GMtot'].isel(time=0)) / sim.data['GMtot'].isel(time=0) + GM_final = sim.data['GM_error'][-1].values + + # Print the final errors + print("Final Angular Momentum Error: ", L_final) + print("Final Energy Error: ", E_final) + print("Final Mass Error: ", GM_final) + + # Determine if the errors are within bounds + L_limit = 1e-10 + E_limit = 1e-5 + GM_limit = 1e-14 + + self.assertLess(L_final,L_limit, msg=f"Angular Momentum Error of {L_final} higher than threshold value of {L_limit}") + self.assertLess(E_final,E_limit, msg=f"Energy Error of {E_final} higher than threshold value of {E_limit}") + self.assertLess(GM_final,GM_limit, msg=f"Mass Error of {GM_final} higher than threshold value of {GM_limit}") + +if __name__ == '__main__': + unittest.main() +