Merge remote-tracking branch 'origin/debug' into debug

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***"
+)
+

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
 ```

examples/Basic_Simulation/.gitignore

@@ -2,4 +2,5 @@
 !.gitignore
 !initial_conditions.py
 !output_reader.py
+!errors.py
 !README.txt

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.
+see the Swiftest User Manual.

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.")

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,31 +48,31 @@
 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)
 
 # Add 10 user-defined test particles.
 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)

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")