Merge branch 'debug' into RevampedRestart

python/swiftest/swiftest/init_cond.py

@@ -286,7 +286,7 @@ def vec2xr(param: Dict, **kwargs: Any):
         if "rot" not in kwargs and "Gmass" in kwargs:
             kwargs['rot'] = np.zeros((len(kwargs['Gmass']),3))
         if "Ip" not in kwargs and "Gmass" in kwargs:
-            kwargs['Ip'] = np.full_like(kwargs['Gmass'], 0.4)
+            kwargs['Ip'] = np.full((len(kwargs['Gmass']),3), 0.4)
 
     if "time" not in kwargs:
         kwargs["time"] = np.array([0.0])

python/swiftest/swiftest/simulation_class.py

@@ -2792,7 +2792,7 @@ def _preprocess(ds, param):
             return io.process_netcdf_input(ds,param)
         partial_func = partial(_preprocess, param=self.param)
 
-        self.collisions = xr.open_mfdataset(col_files,parallel=True, coords=["collision"], join="inner", preprocess=partial_func,mask_and_scale=True)
+        self.collisions = xr.open_mfdataset(col_files,parallel=True,  combine="nested", concat_dim="collision", preprocess=partial_func,mask_and_scale=True)
         self.collisions = io.process_netcdf_input(self.collisions, self.param)
 
         return

src/CMakeLists.txt

@@ -19,15 +19,18 @@ SET(STRICT_MATH_FILES
             ${SRC}/symba/symba_kick.f90
             ${SRC}/whm/whm_kick.f90
             ${SRC}/swiftest/swiftest_user.f90
+            ${SRC}/fraggle/fraggle_generate.f90
+            ${SRC}/fraggle/fraggle_util.f90
+            ${SRC}/fraggle/fraggle_module.f90
+            ${SRC}/misc/lambda_function_module.f90
+            ${SRC}/misc/solver_module.f90
 )
 
 SET(FAST_MATH_FILES
             ${SRC}/globals/globals_module.f90
             ${SRC}/base/base_module.f90
             ${SRC}/netcdf_io/netcdf_io_module.f90
-            ${SRC}/misc/lambda_function_module.f90
             ${SRC}/misc/io_progress_bar_module.f90
-            ${SRC}/misc/solver_module.f90
             ${SRC}/encounter/encounter_module.f90
             ${SRC}/collision/collision_module.f90
             ${SRC}/operator/operator_module.f90
@@ -37,7 +40,6 @@ SET(FAST_MATH_FILES
             ${SRC}/rmvs/rmvs_module.f90
             ${SRC}/helio/helio_module.f90
             ${SRC}/symba/symba_module.f90
-            ${SRC}/fraggle/fraggle_module.f90
             ${SRC}/collision/collision_check.f90  
             ${SRC}/collision/collision_generate.f90
             ${SRC}/collision/collision_io.f90     
@@ -47,8 +49,6 @@ SET(FAST_MATH_FILES
             ${SRC}/encounter/encounter_check.f90
             ${SRC}/encounter/encounter_io.f90
             ${SRC}/encounter/encounter_util.f90
-            ${SRC}/fraggle/fraggle_generate.f90
-            ${SRC}/fraggle/fraggle_util.f90
             ${SRC}/helio/helio_drift.f90
             ${SRC}/helio/helio_gr.f90
             ${SRC}/helio/helio_step.f90
@@ -65,7 +65,6 @@ SET(FAST_MATH_FILES
             ${SRC}/swiftest/swiftest_drift.f90    
             ${SRC}/swiftest/swiftest_gr.f90       
             ${SRC}/swiftest/swiftest_io.f90         
-            ${SRC}/swiftest/swiftest_kick.f90       
             ${SRC}/swiftest/swiftest_obl.f90    
             ${SRC}/swiftest/swiftest_orbel.f90
             ${SRC}/swiftest/swiftest_util.f90
@@ -85,7 +84,6 @@ SET(FAST_MATH_FILES
             ${SRC}/swiftest/swiftest_driver.f90
 )
 
-
 set(SWIFTEST_src ${FAST_MATH_FILES} ${STRICT_MATH_FILES})
 
 # Define the executable in terms of the source files

src/base/base_module.f90

@@ -20,7 +20,7 @@ module base
    !> User defined parameters that are read in from the parameters input file. 
    !>    Each paramter is initialized to a default values. 
    type, abstract :: base_parameters
-      character(len=:), allocatable           :: integrator                             !! Symbolic name of the nbody integrator  used
+      character(len=:), allocatable           :: integrator                             !! Name of the nbody integrator used
       character(len=:), allocatable           :: param_file_name                        !! The name of the parameter file
       integer(I4B)                            :: maxid                = -1              !! The current maximum particle id number 
       integer(I4B)                            :: maxid_collision      = 0               !! The current maximum collision id number

src/collision/collision_io.f90

@@ -303,7 +303,7 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param)
 
             select type(before =>self%collider%before)
             class is (swiftest_nbody_system)
-            select type(after =>self%collider%before)
+            select type(after =>self%collider%after)
             class is (swiftest_nbody_system)
                do stage = 1,2
                   if (allocated(pl)) deallocate(pl)

src/collision/collision_resolve.f90

@@ -316,18 +316,19 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
       real(DP),                  intent(in)    :: t      !! Time of collision
       integer(I4B),              intent(in)    :: status !! Status flag to assign to adds
       ! Internals
-      integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, nimpactors, nfrag
+      integer(I4B) :: i, ibiggest, ismallest, iother, nimpactors, nfrag
       logical, dimension(:), allocatable    :: lmask
       class(swiftest_pl), allocatable           :: plnew, plsub
       character(*), parameter :: FRAGFMT = '("Newbody",I0.7)'
       character(len=NAMELEN) :: newname, origin_type
+      real(DP) :: volume
 
       select type(nbody_system)
       class is (swiftest_nbody_system)
       select type(param)
       class is (swiftest_parameters)
          associate(pl => nbody_system%pl, pl_discards => nbody_system%pl_discards, info => nbody_system%pl%info, pl_adds => nbody_system%pl_adds, cb => nbody_system%cb, npl => pl%nbody, &
-            collision_basic => nbody_system%collider, impactors => nbody_system%collider%impactors,fragments => nbody_system%collider%fragments)
+            collider => nbody_system%collider, impactors => nbody_system%collider%impactors,fragments => nbody_system%collider%fragments)
 
             ! Add the impactors%id bodies to the subtraction list
             nimpactors = impactors%ncoll
@@ -351,7 +352,10 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
             plnew%mass(1:nfrag) = fragments%mass(1:nfrag)
             plnew%Gmass(1:nfrag) = param%GU * fragments%mass(1:nfrag)
             plnew%radius(1:nfrag) = fragments%radius(1:nfrag)
-            plnew%density(1:nfrag) = fragments%mass(1:nfrag) / fragments%radius(1:nfrag)
+            do concurrent(i = 1:nfrag)
+               volume = 4.0_DP/3.0_DP * PI * plnew%radius(i)**3
+               plnew%density(i) = fragments%mass(i) / volume
+            end do
             call plnew%set_rhill(cb)
 
             select case(status)
@@ -377,7 +381,7 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
                if (status == DISRUPTED) then
                   write(origin_type,*) "Disruption"
                else if (status == HIT_AND_RUN_DISRUPT) then
-                  write(origin_type,*) "Hit and run fragmention"
+                  write(origin_type,*) "Hit and run fragmentation"
                end if
                call plnew%info(1)%copy(pl%info(ibiggest))
                plnew%status(1) = OLD_PARTICLE
@@ -437,14 +441,12 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
             end where
 
             ! Log the properties of the new bodies
-            select type(after => collision_basic%after)
+            select type(after => collider%after)
             class is (swiftest_nbody_system)
                allocate(after%pl, source=plnew)
             end select
 
             ! Append the new merged body to the list 
-            nstart = pl_adds%nbody + 1
-            nend = pl_adds%nbody + nfrag
             call pl_adds%append(plnew, lsource_mask=[(.true., i=1, nfrag)])
 
             ! Add the discarded bodies to the discard list
@@ -461,8 +463,6 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
             allocate(plsub, mold=pl)
             call pl%spill(plsub, lmask, ldestructive=.false.)
 
-            nstart = pl_discards%nbody + 1
-            nend = pl_discards%nbody + nimpactors
             call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nimpactors)])
 
             call plsub%setup(0, param)

src/fraggle/fraggle_generate.f90

@@ -51,7 +51,6 @@ module subroutine fraggle_generate(self, nbody_system, param, t)
             call self%set_mass_dist(param) 
             call self%disrupt(nbody_system, param, t)
 
-
             associate (fragments => self%fragments)
                ! Populate the list of new bodies
                nfrag = fragments%nbody
@@ -96,6 +95,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
       logical                              :: lk_plpl, lfailure_local
       logical, dimension(size(IEEE_ALL))   :: fpe_halting_modes, fpe_quiet_modes
       real(DP)                             :: dE, dL
+      integer(I4B)                         :: i
       character(len=STRMAX)                :: message
       real(DP), parameter                  :: fail_scale_initial = 1.001_DP
 
@@ -135,6 +135,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
          call fraggle_generate_vel_vec(self,lfailure_local)
          call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
          call self%set_original_scale()
+
          dE = self%Etot(2) - self%Etot(1)
          dL = .mag.(self%Ltot(:,2) - self%Ltot(:,1))
 
@@ -514,14 +515,14 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure)
          lfailure = E_residual < 0.0_DP
 
          do concurrent(i = 1:nfrag)
-            fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:)
+            collider%fragments%vb(:,i) = collider%fragments%vc(:,i) + impactors%vbcom(:)
          end do
 
          impactors%vbcom(:) = 0.0_DP
          do concurrent(i = 1:nfrag)
-            impactors%vbcom(:) = impactors%vbcom(:) + fragments%mass(i) * fragments%vb(:,i) 
+            impactors%vbcom(:) = impactors%vbcom(:) + collider%fragments%mass(i) * collider%fragments%vb(:,i) 
          end do
-         impactors%vbcom(:) = impactors%vbcom(:) / fragments%mtot
+         impactors%vbcom(:) = impactors%vbcom(:) / collider%fragments%mtot
 
       end associate
       return

src/fraggle/fraggle_util.f90

@@ -188,7 +188,7 @@ module subroutine fraggle_util_set_mass_dist(self, param)
 
             fragments%density(istart:nfrag) = fragments%mtot / sum(volume(:))
             fragments%radius(istart:nfrag) = (3 * fragments%mass(istart:nfrag) / (4 * PI * fragments%density(istart:nfrag)))**(1.0_DP / 3.0_DP)
-            do i = istart, nfrag
+            do concurrent(i = istart:nfrag)
                fragments%Ip(:, i) = Ip_avg(:)
             end do
 
@@ -263,7 +263,7 @@ module subroutine fraggle_util_set_natural_scale_factors(self)
          impactors%Lspin(:,:)     = impactors%Lspin(:,:)     / collider%Lscale
          impactors%Lorbit(:,:)    = impactors%Lorbit(:,:)    / collider%Lscale
 
-         do i = 1, 2
+         do concurrent(i = 1:2)
             impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i))
          end do
 
@@ -311,7 +311,7 @@ module subroutine fraggle_util_set_original_scale_factors(self)
          impactors%vc        = impactors%vc        * collider%vscale
          impactors%Lspin     = impactors%Lspin     * collider%Lscale
          impactors%Lorbit    = impactors%Lorbit    * collider%Lscale
-         do i = 1, 2
+         do concurrent(i = 1:2)
             impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i))
          end do
 
@@ -322,11 +322,8 @@ module subroutine fraggle_util_set_original_scale_factors(self)
          fragments%rot(:,:)  = fragments%rot(:,:)    / collider%tscale
          fragments%rc(:,:)   = fragments%rc(:,:)     * collider%dscale
          fragments%vc(:,:)   = fragments%vc(:,:)     * collider%vscale
-
-         do i = 1, fragments%nbody
-            fragments%rb(:, i) = fragments%rc(:, i) + impactors%rbcom(:)
-            fragments%vb(:, i) = fragments%vc(:, i) + impactors%vbcom(:)
-         end do
+         fragments%rb(:,:)   = fragments%rb(:,:)     * collider%dscale
+         fragments%vb(:,:)   = fragments%vb(:,:)     * collider%vscale
 
          impactors%Qloss = impactors%Qloss * collider%Escale
 

src/swiftest/swiftest_discard.f90

@@ -25,7 +25,7 @@ module subroutine swiftest_discard_system(self, param)
       lpl_check = allocated(self%pl_discards)
       ltp_check = allocated(self%tp_discards)
 
-      associate(nbody_system => self, tp => self%tp, pl => self%pl, tp_discards => self%tp_discards, pl_discards => self%pl_discards)
+      associate(nbody_system => self, tp => self%tp, pl => self%pl, tp_discards => self%tp_discards, pl_discards => self%pl_discards, nc => param%system_history%nc)
          lpl_discards = .false.
          ltp_discards = .false.
          if (lpl_check) then
@@ -37,11 +37,20 @@ module subroutine swiftest_discard_system(self, param)
             call tp%discard(nbody_system, param)
             ltp_discards = (tp_discards%nbody > 0)
          end if
-         if (ltp_discards) call tp_discards%write_info(param%system_history%nc, param)
-         if (lpl_discards) call pl_discards%write_info(param%system_history%nc, param)
-         if (lpl_discards .and. param%lenergy) call self%conservation_report(param, lterminal=.false.)
-         if (lpl_check) call pl_discards%setup(0,param) 
-         if (ltp_check) call tp_discards%setup(0,param) 
+
+         if (ltp_discards.or.lpl_discards) then
+            call nc%open(param)
+            if (lpl_discards) then
+               call pl_discards%write_info(param%system_history%nc, param)
+               if (param%lenergy) call self%conservation_report(param, lterminal=.false.)
+               call pl_discards%setup(0,param) 
+            end if
+            if (ltp_discards) then
+               call tp_discards%write_info(param%system_history%nc, param)
+               call tp_discards%setup(0,param) 
+            end if
+            call nc%close()
+         end if
 
       end associate
 

src/swiftest/swiftest_driver.f90

@@ -46,7 +46,7 @@ program swiftest_driver
    !> Read in the user-defined parameters file and the initial conditions of the nbody_system
    allocate(swiftest_parameters :: param)
    param%integrator = trim(adjustl(integrator))
-   call param%set_display(display_style)
+   param%display_style = trim(adjustl(display_style))
    call param%read_in(param_file_name)
 
    associate(t0       => param%t0, &