diff --git a/buildscripts/build_dependencies.sh b/buildscripts/build_dependencies.sh index c922824d8..a09b99922 100755 --- a/buildscripts/build_dependencies.sh +++ b/buildscripts/build_dependencies.sh @@ -41,6 +41,11 @@ if ! command -v ninja &> /dev/null; then fi fi +# Get the OpenMP Libraries +if [ $OS = "MacOSX" ]; then + ${SCRIPT_DIR}/get_lomp.sh ${ARGS} +fi + ${SCRIPT_DIR}/build_zlib.sh ${ARGS} ${SCRIPT_DIR}/build_hdf5.sh ${ARGS} ${SCRIPT_DIR}/build_netcdf-c.sh ${ARGS} diff --git a/buildscripts/build_hdf5.sh b/buildscripts/build_hdf5.sh index 9863f33ed..3f93d036b 100755 --- a/buildscripts/build_hdf5.sh +++ b/buildscripts/build_hdf5.sh @@ -12,7 +12,6 @@ # If not, see: https://www.gnu.org/licenses. HDF5_VER="1_14_2" -PLUGIN_VER="1.14.0" ZLIB_VER="1.3" SCRIPT_DIR=$(realpath $(dirname $0)) @@ -20,10 +19,7 @@ set -a ARGS=$@ . ${SCRIPT_DIR}/_build_getopts.sh ${ARGS} . ${SCRIPT_DIR}/set_compilers.sh -# Get the OpenMP Libraries -if [ $OS = "MacOSX" ]; then - ${SCRIPT_DIR}/get_lomp.sh ${ARGS} -fi + NPROC=$(nproc) @@ -54,7 +50,6 @@ fi if [ ! -d ${DEPENDENCY_DIR}/hdfsrc ]; then curl -s -L https://github.com/HDFGroup/hdf5/releases/download/hdf5-${HDF5_VER}/hdf5-${HDF5_VER}.tar.gz | tar xvz -C ${DEPENDENCY_DIR} fi - printf "\n" printf "*********************************************************\n" printf "* BUILDING HDF5 LIBRARY *\n" @@ -68,9 +63,6 @@ printf "LDFLAGS: ${LDFLAGS}\n" printf "*********************************************************\n" cd ${DEPENDENCY_DIR}/hdfsrc -ZLIB_TGZ_NAME="zlib-${ZLIB_VER}.tar.gz" -ZLIB_TGZ_ORIGPATH="https://github.com/madler/zlib/releases/download/v${ZLIB_VER}/" -curl -L "https://github.com/HDFGroup/hdf5_plugins/archive/refs/tags/${PLUGIN_VER}.tar.gz" -o hdf5_plugins.tar.gz HDF5_ROOT=${PREFIX} ZLIB_ROOT=${PREFIX} @@ -108,18 +100,6 @@ else sudo cmake --install build fi -# tar xvzf hdf5_plugins.tar.gz -# PLUGIN_SOURCE=hdf5_plugins-${PLUGIN_VER} - -# BUILD_OPTIONS="-DTGZPATH:PATH=${PLUGIN_SOURCE}/libs -DH5PL_ALLOW_EXTERNAL_SUPPORT:STRING=\"TGZ\"" -# cmake -B ${PLUGIN_SOURCE}/build -C ${PLUGIN_SOURCE}/config/cmake/cacheinit.cmake -DCMAKE_BUILD_TYPE:STRING=Release ${BUILD_OPTIONS} -G Ninja ${PLUGIN_SOURCE} -# cmake --build ${PLUGIN_SOURCE}/build -j${NPROC} --config Release -# if [ -w ${PREFIX} ]; then -# cmake --install ${PLUGIN_SOURCE}/build -# else -# sudo cmake --install ${PLUGIN_SOURCE}/build -# fi - if [ $? -ne 0 ]; then printf "hdf5 could not be compiled.\n" exit 1 diff --git a/buildscripts/build_netcdf-c.sh b/buildscripts/build_netcdf-c.sh index 6295eeda7..ee7415cbf 100755 --- a/buildscripts/build_netcdf-c.sh +++ b/buildscripts/build_netcdf-c.sh @@ -15,10 +15,6 @@ set -a ARGS=$@ . ${SCRIPT_DIR}/_build_getopts.sh ${ARGS} . ${SCRIPT_DIR}/set_compilers.sh -# Get the OpenMP Libraries -if [ $OS = "MacOSX" ]; then - ${SCRIPT_DIR}/get_lomp.sh ${ARGS} -fi NPROC=$(nproc) diff --git a/buildscripts/build_netcdf-fortran.sh b/buildscripts/build_netcdf-fortran.sh index 0457e74c7..c82f24573 100755 --- a/buildscripts/build_netcdf-fortran.sh +++ b/buildscripts/build_netcdf-fortran.sh @@ -15,10 +15,6 @@ set -a ARGS=$@ . ${SCRIPT_DIR}/_build_getopts.sh ${ARGS} . ${SCRIPT_DIR}/set_compilers.sh -# Get the OpenMP Libraries -if [ $OS = "MacOSX" ]; then - ${SCRIPT_DIR}/get_lomp.sh ${ARGS} -fi NPROC=$(nproc) diff --git a/buildscripts/build_zlib.sh b/buildscripts/build_zlib.sh index c3ef95318..337e6a839 100755 --- a/buildscripts/build_zlib.sh +++ b/buildscripts/build_zlib.sh @@ -15,10 +15,6 @@ set -a ARGS=$@ . ${SCRIPT_DIR}/_build_getopts.sh ${ARGS} . ${SCRIPT_DIR}/set_compilers.sh -# Get the OpenMP Libraries -if [ $OS = "MacOSX" ]; then - ${SCRIPT_DIR}/get_lomp.sh ${ARGS} -fi NPROC=$(nproc) diff --git a/pyproject.toml b/pyproject.toml index 106f976f9..7de757846 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "swiftest" -version = "2023.9.4" +version = "2023.10.0" authors=[ {name = 'David A. Minton', email='daminton@purdue.edu'}, {name = 'Carlisle Wishard'}, @@ -30,8 +30,9 @@ dependencies = [ 'dask>=2022.1', 'distributed>=2022.1', 'bottleneck>=1.3.5', - 'h5netcdf>=1.0.2', - 'netcdf4>=1.6.2', + 'h5netcdf', + 'h5py', + 'netcdf4', 'matplotlib>=3.7.1', 'astropy>=5.1', 'astroquery>=0.4.6', @@ -103,7 +104,7 @@ before-all = [ ] [tool.cibuildwheel.linux] -skip = "cp312-*" +skip = "cp312-* pp* -manylinux_i686* *-musllinux*" before-all = [ "yum install doxygen libxml2-devel libcurl-devel fftw-devel blas-devel lapack-devel -y", "buildscripts/build_dependencies.sh -p /usr/local" diff --git a/src/globals/globals_module.f90 b/src/globals/globals_module.f90 index 93ed3ac1e..d1a078329 100644 --- a/src/globals/globals_module.f90 +++ b/src/globals/globals_module.f90 @@ -48,7 +48,7 @@ module globals integer(I4B), parameter :: UPPERCASE_OFFSET = iachar('A') - iachar('a') !! ASCII character set parameter for lower to upper !! conversion - offset between upper and lower - character(*), parameter :: VERSION = "2023.9.4" !! Swiftest version + character(*), parameter :: VERSION = "2023.10.0" !! Swiftest version !> Symbolic name for integrator types character(*), parameter :: UNKNOWN_INTEGRATOR = "UKNOWN INTEGRATOR" diff --git a/swiftest/simulation_class.py b/swiftest/simulation_class.py index baa59c04d..b71f01a67 100644 --- a/swiftest/simulation_class.py +++ b/swiftest/simulation_class.py @@ -412,15 +412,16 @@ def run(self,dask: bool = False, **kwargs): if len(kwargs) > 0: self.set_parameter(**kwargs) - # Write out the current parameter set before executing run - self.write_param() - if self.codename != "Swiftest": warnings.warn(f"Running an integration is not yet supported for {self.codename}",stacklevel=2) return + # Save initial conditions if not self.restart: self.clean() + + # Write out the current parameter set before executing run + self.write_param(verbose=False,**kwargs) print(f"Running a {self.codename} {self.integrator} run from tstart={self.param['TSTART']} {self.TU_name} to tstop={self.param['TSTOP']} {self.TU_name}") @@ -2176,10 +2177,8 @@ def add_solar_system_body(self, dsnew = init_cond.vec2xr(self.param,**kwargs) dsnew = self._combine_and_fix_dsnew(dsnew) - if dsnew['npl'] > 0 or dsnew['ntp'] > 0: + if dsnew['id'].max(dim='name') > 0 and dsnew['name'].size > 0: self.save(verbose=False) - - self.init_cond = self.data.copy(deep=True) diff --git a/tests/test_swiftest.py b/tests/test_swiftest.py index e6d1cf4c9..63f07829c 100755 --- a/tests/test_swiftest.py +++ b/tests/test_swiftest.py @@ -1,5 +1,4 @@ #!/usr/bin/env python3 - """ Copyright 2023 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh This file is part of Swiftest. @@ -25,22 +24,19 @@ import swiftest import unittest import os -import shutil import numpy as np from numpy.random import default_rng from astroquery.jplhorizons import Horizons import datetime - 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): + def test_gen_ic(self): """ Tests that Swiftest is able to successfully generate a set of initial conditions in a file without any exceptions being raised """ @@ -51,99 +47,56 @@ def test01_gen_ic(self): sim = swiftest.Simulation() sim.clean() -#!/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 shutil -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)) - + return + - def test02_read_ic(self): + def test_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) + sim = swiftest.Simulation() + sim.clean() + + # Add the modern planets and the Sun using the JPL Horizons Database. + sim.add_solar_system_body(major_bodies) + sim.save() # 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])) - + return + - def test03_integrators(self): + def test_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) + sim = swiftest.Simulation() + + # Add the modern planets and the Sun using the JPL Horizons Database. + sim.add_solar_system_body(major_bodies) + sim.clean() # 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) + 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) @@ -154,9 +107,9 @@ def test03_integrators(self): sim.run(integrator=i) except: self.fail(f"Failed with integrator {i}") + return - - def test04_conservation(self): + def test_conservation(self): """ Tests that Swiftest conserves mass, energy, and momentum to within acceptable tolerances. """ @@ -214,9 +167,9 @@ def fit_func(x,slope,b): 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}") - + return - def test05_gr(self): + def test_gr(self): """ Tests that GR is working correctly by computing the precession of Mercury's longitude of periapsis and comparing it to values obtained from the JPL/Horizons ephemeris service @@ -253,129 +206,12 @@ def test05_gr(self): varpi_sim = sim.data['varpi'].sel(name="Mercury") dvarpi_gr = np.diff(varpi_sim) / tstep_out dvarpi_err = np.mean(dvarpi_obs - dvarpi_gr) / dvarpi_obs_mean - self.assertLess(np.abs(dvarpi_err),dvarpi_limit,msg=f'Mercury precession rate error of of {dvarpi_err:.2e} "/{sim.TU_name} higher than threshold value of {dvarpi_limit:.2e} "/{sim.TU_name}') + print(f'{i}: Mercury precession rate error {dvarpi_err:.2e} "/{sim.TU_name}') + self.assertLess(np.abs(dvarpi_err),dvarpi_limit,msg=f'{dvarpi_err:.2e} /{sim.TU_name} is higher than threshold value of {dvarpi_limit:.2e} "/{sim.TU_name}') return 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") + os.environ["HDF5_USE_FILE_LOCKING"]="FALSE" + unittest.main() \ No newline at end of file diff --git a/version.txt b/version.txt index 6895b410d..3bae60812 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -2023.9.4 \ No newline at end of file +2023.10.0 \ No newline at end of file