Merge branch 'debug'

cmake/Modules/SetFortranFlags.cmake

@@ -330,4 +330,15 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_PROFILE "${CMAKE_Fortran_FLAGS_RELEASE}"
 
 SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_PROFILE "${CMAKE_Fortran_FLAGS_PROFILE}"
                  Fortran "-check bounds,pointers,uninit"  # Intel
-                )                
+                         "-fcheck=bounds,pointer,mem"
+                )
+
+# Sanitize
+SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}"
+                 Fortran "-fsanitize=address,undefined"  # Gnu 
+                )
+
+
+SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}"
+                 Fortran "-fstack-check" # GNU 
+                )     

python/swiftest/swiftest/simulation_class.py

@@ -28,6 +28,7 @@
 import subprocess
 import shlex
 import warnings
+import sys
 from tqdm.auto import tqdm
 from typing import (
     Literal,
@@ -453,9 +454,11 @@ def _type_scrub(output_data):
                 if p.returncode != 0:
                     for line in res[1]:
                         print(line, end='')
-                        warnings.warn("Failure in swiftest_driver", stacklevel=2)
+                    warnings.warn("Failure in swiftest_driver", stacklevel=2)
+                    sys.exit()
         except:
             warnings.warn(f"Error executing main swiftest_driver program", stacklevel=2)
+            sys.exit()
 
         pbar.close()
         return
@@ -485,11 +488,20 @@ def run(self,**kwargs):
             warnings.warn(f"Running an integration is not yet supported for {self.codename}",stacklevel=2)
             return
 
-        if self.driver_executable is None:
+        if self.driver_executable is None or self.binary_source is None:
             msg = "Path to swiftest_driver has not been set!"
             msg += f"\nMake sure swiftest_driver is compiled and the executable is in {str(self.binary_path)}"
             warnings.warn(msg,stacklevel=2)
             return
+
+        if not self.driver_executable.exists():
+            if not self.binary_source.exists():
+                msg = "Path to swiftest_driver has not been set!"
+                msg += f"\nMake sure swiftest_driver is compiled and the executable is in {str(self.binary_path)}"
+                warnings.warn(msg,stacklevel=2)
+                return
+            else:
+                shutil.copy(self.binary_source, self.driver_executable) 
 
         if not self.restart:
             self.clean()
@@ -902,11 +914,16 @@ def set_integrator(self,
             self.param['! VERSION'] = f"{self.codename} input file"
             update_list.append("codename")
             if self.codename == "Swiftest":
-                self.binary_path = Path(_pyfile).parent.parent.parent.parent / "bin"
+                self.binary_source = Path(_pyfile).parent.parent.parent.parent / "bin" / "swiftest_driver"
+                self.binary_path = self.simdir.resolve()
                 self.driver_executable = self.binary_path / "swiftest_driver"
-                if not self.driver_executable.exists():
+                if not self.binary_source.exists():
                     warnings.warn(f"Cannot find the Swiftest driver in {str(self.binary_path)}",stacklevel=2)
                     self.driver_executable = None
+                else:
+                    if self.binary_path.exists():
+                        shutil.copy(self.binary_source, self.driver_executable) 
+                        self.driver_executable.resolve()
             else:
                 self.binary_path = "NOT IMPLEMENTED FOR THIS CODE"
                 self.driver_executable = None
@@ -1027,6 +1044,7 @@ def set_feature(self,
                     encounter_check_loops: Literal["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] | None = None,
                     encounter_save: Literal["NONE", "TRAJECTORY", "CLOSEST", "BOTH"] | None = None,
                     verbose: bool | None = None,
+                    simdir: str | os.PathLike = None, 
                     **kwargs: Any
                     ):
         """
@@ -1105,6 +1123,9 @@ def set_feature(self,
             If true, will restart an old run. The file given by `output_file_name` must exist before the run can
             execute. If false, will start a new run. If the file given by `output_file_name` exists, it will be replaced
             when the run is executed.
+        simdir: PathLike, optional
+            Directory where simulation data will be stored, including the parameter file, initial conditions file, output file,
+            dump files, and log files. 
         verbose: bool, optional
             If passed, it will override the Simulation object's verbose flag
         **kwargs
@@ -1222,6 +1243,16 @@ def set_feature(self,
             else:
                 self.param["ENCOUNTER_SAVE"] = encounter_save
                 update_list.append("encounter_save")
+
+        if simdir is not None:
+            self.simdir = Path(simdir)
+            if self.simdir.exists():
+                if not self.simdir.is_dir():
+                    msg = f"Cannot create the {self.simdir.resolve()} directory: File exists."
+                    msg += "\nDelete the file or change the location of param_file"
+                    raise NotADirectoryError(msg)
+            self.binary_path = self.simdir.resolve()
+            self.driver_executable = self.binary_path / "swiftest_driver"
 
         self.param["TIDES"] = False
 
@@ -2868,10 +2899,22 @@ def save(self,
             param = self.param
 
         self.simdir.mkdir(parents=True, exist_ok=True)
+
+
         if codename == "Swiftest":
             infile_name = Path(self.simdir) / param['NC_IN']
             io.swiftest_xr2infile(ds=self.data, param=param, in_type=self.param['IN_TYPE'], infile_name=infile_name, framenum=framenum, verbose=verbose)
             self.write_param(param_file=param_file,**kwargs)
+
+            if not self.driver_executable.exists():
+                if not self.binary_source.exists():
+                    msg = "Path to swiftest_driver has not been set!"
+                    msg += f"\nMake sure swiftest_driver is compiled and the executable is in {str(self.binary_path)}"
+                    warnings.warn(msg,stacklevel=2)
+                    return
+                else:
+                    shutil.copy(self.binary_source, self.driver_executable)  
+
         elif codename == "Swifter":
             if codename == "Swiftest":
                 swifter_param = io.swiftest2swifter_param(param)

src/fraggle/fraggle_generate.f90

@@ -475,7 +475,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
       integer(I4B), dimension(collider%fragments%nbody) :: vsign
       real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove
       ! For the initial "guess" of fragment velocities, this is the minimum and maximum velocity relative to escape velocity that the fragments will have
-      real(DP)                :: vmin_guess = 1.5_DP 
+      real(DP)                :: vmin_guess = 1.1_DP 
       real(DP)                :: vmax_guess = 10.0_DP
       real(DP)                :: delta_v, volume
       integer(I4B), parameter :: MAXLOOP = 200

src/swiftest/swiftest_driver.f90

@@ -79,8 +79,8 @@ program swiftest_driver
             else
                call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum
             end if
+            call nbody_system%conservation_report(param, lterminal=.true.)
          end if
-         call nbody_system%conservation_report(param, lterminal=.true.)
          call system_history%take_snapshot(param,nbody_system)
          call nbody_system%dump(param)
 

src/swiftest/swiftest_io.f90

@@ -800,12 +800,14 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param)
          end do
 
          ! Set special fill mode for discard time so that we can make use of it for non-discarded bodies.
-         select case (vartype)
-         case(NF90_FLOAT)
-            call netcdf_io_check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_SP)), "netcdf_io_initialize_output nf90_def_var_fill discard_time NF90_FLOAT"  )
-         case(NF90_DOUBLE)
-            call netcdf_io_check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_DP)), "netcdf_io_initialize_output nf90_def_var_fill discard_time NF90_DOUBLE"  )
-         end select
+         if (param%lclose) then
+            select case (vartype)
+            case(NF90_FLOAT)
+               call netcdf_io_check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_SP)), "netcdf_io_initialize_output nf90_def_var_fill discard_time NF90_FLOAT"  )
+            case(NF90_DOUBLE)
+               call netcdf_io_check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_DP)), "netcdf_io_initialize_output nf90_def_var_fill discard_time NF90_DOUBLE"  )
+            end select
+         end if
 
          ! Take the file out of define mode
          call netcdf_io_check( nf90_enddef(nc%id), "netcdf_io_initialize_output nf90_enddef"  )
@@ -1112,13 +1114,12 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier
             call netcdf_io_check( nf90_get_var(nc%id, nc%radius_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_io_read_frame_system nf90_getvar radius_varid"  )
             cb%radius = rtemp(1)
 
-            ! Set initial central body radius for SyMBA bookkeeping
-            cb%R0 = cb%radius
             if (npl > 0) pl%radius(:) = pack(rtemp, plmask)
          else
             cb%radius = param%rmin
             if (npl > 0) pl%radius(:) = 0.0_DP
          end if
+         cb%R0 = cb%radius
 
          if (param%lrotation) then
             call netcdf_io_check( nf90_get_var(nc%id, nc%Ip_varid,  vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_io_read_frame_system nf90_getvar Ip_varid"  )

src/swiftest/swiftest_util.f90

@@ -1446,7 +1446,11 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param)
             nbody_system%L_orbit(j) = Lcborbit(j) + sum(Lplorbit(j,1:npl), pl%lmask(1:npl)) 
          end do
 
-         nbody_system%be = sum(-3*pl%Gmass(1:npl)*pl%mass(1:npl)/(5*pl%radius(1:npl)), pl%lmask(1:npl)) 
+         if ((param%lclose)) then
+            nbody_system%be = sum(-3*pl%Gmass(1:npl)*pl%mass(1:npl)/(5*pl%radius(1:npl)), pl%lmask(1:npl)) 
+         else
+            nbody_system%be = 0.0_DP
+         end if
          nbody_system%te = nbody_system%ke_orbit + nbody_system%ke_spin + nbody_system%pe + nbody_system%be 
          nbody_system%L_total(:) = nbody_system%L_orbit(:) + nbody_system%L_spin(:)
       end associate