diff --git a/CMakeLists.txt b/CMakeLists.txt index 961683f3e..4d36d9612 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -73,3 +73,17 @@ ADD_SUBDIRECTORY(${SRC} ${BIN}) ADD_CUSTOM_TARGET(distclean COMMAND ${CMAKE_COMMAND} -P ${CMAKE_SOURCE_DIR}/distclean.cmake ) + +# 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 "***Test Complete***" +) + diff --git a/README.md b/README.md index 3671832cb..f45a812c3 100644 --- a/README.md +++ b/README.md @@ -46,7 +46,7 @@ Parallelization in Swiftest is done with OpenMP. Version 3.1.4 or higher is nece The easiest way to get Swiftest on your machine is to clone the GitHub repository. To do so, open a terminal window and type the following: ``` -$ git clone https://github.itap.purdue.edu/MintonGroup/swiftest.git +$ git clone https://github.com/carlislewishard/swiftest.git ``` If your cloned version is not already set to the master branch: @@ -85,7 +85,7 @@ CMake allows the user to specify a set of compiler flags to use during compilati As a general rule, the release flags are fully optimized and best used when running Swiftest with the goal of generating results. This is the default set of flags. When making changes to the Swiftest source code, it is best to compile Swiftest using the debug set of flags. You may also define your own set of compiler flags. -To build Swiftest with the release flags (default), type the following: +To build Swiftest with the release flags (default) using the Intel fortran compiler (ifort), type the following: ``` $ cmake .. ``` @@ -97,8 +97,16 @@ To build with another set of flags, simply replace ```DEBUG``` in the above line Add ```-CMAKE_PREFIX_PATH=/path/to/netcdf/``` to these commands as needed. -After building Swiftest, make the executable using: +If using the GCC fortran compiler (gfortran), add the following flags: +``` +-DCMAKE_Fortran_FLAGS="-I/usr/lib64/gfortran/modules/ -ffree-line-length-512" +``` +You can manually specify the compiler you wish to use with the following flag: +``` +c-DCMAKE_Fortran_COMPILER=$(which ifort) +``` +After building Swiftest, make the executable using: ``` $ make ``` diff --git a/examples/Basic_Simulation/.gitignore b/examples/Basic_Simulation/.gitignore index 728fe8873..9e4d02aad 100644 --- a/examples/Basic_Simulation/.gitignore +++ b/examples/Basic_Simulation/.gitignore @@ -2,4 +2,5 @@ !.gitignore !initial_conditions.py !output_reader.py +!errors.py !README.txt diff --git a/examples/Basic_Simulation/README.txt b/examples/Basic_Simulation/README.txt index 9db2d5f54..ae626a9c5 100644 --- a/examples/Basic_Simulation/README.txt +++ b/examples/Basic_Simulation/README.txt @@ -11,13 +11,14 @@ README.txt Swiftest Example : Basic_Simulation Author : Carlisle Wishard and David Minton -Date : December 6, 2022 +Date : June 27, 2023 Included in the Basic_Simulation example directory are the following files: - README.txt : This file - initial_conditions.py : A Python Script that generates and runs a set of initial conditions. - output_reader.py : A Python Script that processes out.nc and generates output.eps + - errors.py : A Python Script that processes out.nc and reports the simulation errors to the terminal This example is intended to be run with Swiftest SyMBA. For details on how to generate, run, and analyze this example, -see the Swiftest User Manual. \ No newline at end of file +see the Swiftest User Manual. diff --git a/examples/Basic_Simulation/errors.py b/examples/Basic_Simulation/errors.py new file mode 100644 index 000000000..5cbdfaa55 --- /dev/null +++ b/examples/Basic_Simulation/errors.py @@ -0,0 +1,69 @@ +#!/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. +""" + +""" +Reads and processes a Swiftest output file. + +Input +------ +data.nc : A NetCDF file containing the simulation output. + +Output +------ +None +""" + +import swiftest +import xarray as xr +import matplotlib.pyplot as plt + +# Read in the simulation output and store it as an Xarray dataset. +ds = swiftest.Simulation(read_data=True) + +# Calculate the angular momentum error +ds.data['L_tot'] = ds.data['L_orbit'] + ds.data['L_spin'] + ds.data['L_escape'] +ds.data['DL'] = ds.data['L_tot'] - ds.data['L_tot'].isel(time=0) +ds.data['L_error'] = swiftest.tool.magnitude(ds.data,'DL') / swiftest.tool.magnitude(ds.data.isel(time=0), 'L_tot') +L_final = ds.data['L_error'][-1].values + +# Calculate the energy error +ds.data['E_error'] = (ds.data['TE'] - ds.data['TE'].isel(time=0)) / ds.data['TE'].isel(time=0) +E_final = ds.data['E_error'][-1].values + +# Calculate the mass error +ds.data['GMtot'] = ds.data['Gmass'].sum(dim='name',skipna=True) + ds.data['GMescape'] +ds.data['GM_error'] = (ds.data['GMtot'] - ds.data['GMtot'].isel(time=0)) / ds.data['GMtot'].isel(time=0) +GM_final = ds.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 + +lerror = 0 +if (L_final > L_limit): + lerror = 1 + print("Angular Momentum Error of", L_final, " higher than threshold value of", L_limit, ". Test failed.") +if (E_final > E_limit): + lerror = 1 + print("Energy Error of", E_final, " higher than threshold value of", E_limit, ". Test failed.") +if (GM_final > GM_limit): + lerror = 1 + print("Mass Error of", GM_final, " higher than threshold value of", GM_limit, ". Test failed.") +if (lerror == 0): + print("Test successful.") diff --git a/examples/Basic_Simulation/initial_conditions.py b/examples/Basic_Simulation/initial_conditions.py index 813bd8dfb..83614fb07 100755 --- a/examples/Basic_Simulation/initial_conditions.py +++ b/examples/Basic_Simulation/initial_conditions.py @@ -39,6 +39,7 @@ #sim = swiftest.Simulation(fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8) sim = swiftest.Simulation() sim.clean() +rng = default_rng(seed=123) # Add the modern planets and the Sun using the JPL Horizons Database. sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune","Pluto"]) @@ -47,18 +48,18 @@ npl = 5 density_pl = 3000.0 / (sim.param['MU2KG'] / sim.param['DU2M'] ** 3) -name_pl = ["MassiveBody_01", "MassiveBody_02", "MassiveBody_03", "MassiveBody_04", "MassiveBody_05"] -a_pl = default_rng().uniform(0.3, 1.5, npl) -e_pl = default_rng().uniform(0.0, 0.3, npl) -inc_pl = default_rng().uniform(0.0, 90, npl) -capom_pl = default_rng().uniform(0.0, 360.0, npl) -omega_pl = default_rng().uniform(0.0, 360.0, npl) -capm_pl = default_rng().uniform(0.0, 360.0, npl) -M_pl = np.array([6e23, 8e23, 1e24, 3e24, 5e24]) * sim.KG2MU +name_pl = ["SemiBody_01", "SemiBody_02", "SemiBody_03", "SemiBody_04", "SemiBody_05"] +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.01 * np.max(M_pl) +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) @@ -66,12 +67,12 @@ ntp = 10 name_tp = ["TestParticle_01", "TestParticle_02", "TestParticle_03", "TestParticle_04", "TestParticle_05", "TestParticle_06", "TestParticle_07", "TestParticle_08", "TestParticle_09", "TestParticle_10"] -a_tp = default_rng().uniform(0.3, 1.5, ntp) -e_tp = default_rng().uniform(0.0, 0.3, ntp) -inc_tp = default_rng().uniform(0.0, 90, ntp) -capom_tp = default_rng().uniform(0.0, 360.0, ntp) -omega_tp = default_rng().uniform(0.0, 360.0, ntp) -capm_tp = default_rng().uniform(0.0, 360.0, ntp) +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) 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=1.0e3, dt=0.01, istep_out=100, dump_cadence=0, compute_conservation_values=True, mtiny=mtiny) diff --git a/examples/Basic_Simulation/output_reader.py b/examples/Basic_Simulation/output_reader.py index 851384bf3..a46b6cefc 100644 --- a/examples/Basic_Simulation/output_reader.py +++ b/examples/Basic_Simulation/output_reader.py @@ -39,8 +39,7 @@ ax.set(xlabel="Semimajor Axis (AU)", ylabel="Eccentricity", title="Simulation Initial Conditions (t=0)") ax.scatter(sim.data['a'].isel(time=0), sim.data['e'].isel(time=0), c=colors, s=sizes, edgecolor='black') ax.set_xlim(0, 2.0) -ax.set_ylim(0, 0.4) -ax.text(1.5, 0.35, f"t = {sim.data['time'].isel(time=0).values} years", size=10, ha="left") +ax.set_ylim(0, 0.25) +ax.text(1.5, 0.2, f"t = {sim.data['time'].isel(time=0).values} years", size=10, ha="left") plt.tight_layout() -plt.show() -fig.savefig("output.eps", dpi=300, facecolor='white', transparent=False, bbox_inches="tight") +fig.savefig("../examples/Basic_Simulation/output.eps", dpi=300, facecolor='white', transparent=False, bbox_inches="tight")