From e85df57540bb649a15d9da77257adaee757c46fb Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 21 Feb 2024 19:39:53 -0500 Subject: [PATCH] Started work on having pl-tp collisions recorded in the collision.nc file. Still need to get tp snapshots saved --- src/collision/collision_io.f90 | 2 + src/collision/collision_resolve.f90 | 51 +++++++++++++ src/collision/collision_util.f90 | 109 +++++++++++++++------------- src/swiftest/swiftest_discard.f90 | 2 +- 4 files changed, 112 insertions(+), 52 deletions(-) diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index ea8675089..907eb90ec 100644 --- a/src/collision/collision_io.f90 +++ b/src/collision/collision_io.f90 @@ -389,8 +389,10 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) if (allocated(pl)) deallocate(pl) select case(stage) case(1) + if (.not. allocated(before%pl)) cycle allocate(pl, source=before%pl) case(2) + if (.not. allocated(after%pl)) cycle allocate(pl, source=after%pl) end select npl = pl%nbody diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 545b882fc..ed676b32b 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -707,6 +707,14 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec) real(DP), intent(in) :: t !! Current simulation tim real(DP), intent(in) :: dt !! Current simulation step size integer(I4B), intent(in) :: irec !! Current recursion level + ! Internals + class(swiftest_pl), allocatable :: plsub + logical :: lpltp_collision + character(len=STRMAX) :: timestr, idstr + integer(I4B) :: i, j, nnew, loop + integer(I8B) :: k, ncollisions + integer(I4B), dimension(:), allocatable :: idnew + logical, dimension(:), allocatable :: lmask ! Make sure coordinate systems are all synced up due to being inside the recursion at this point select type(nbody_system) @@ -726,6 +734,49 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec) ! Discard the collider call nbody_system%tp%discard(nbody_system, param) + + associate(idx1 => pltp_collision%index1, idx2 => pltp_collision%index2) + ncollisions = pltp_collision%nenc + write(timestr,*) t + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "") + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "***********************************************************" // & + "***********************************************************") + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Collision between test particle and massive body detected at time t = " // & + trim(adjustl(timestr))) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "***********************************************************" // & + "***********************************************************") + + do k = 1_I8B, ncollisions + ! Advance the collision id number and save it + collider%maxid_collision = max(collider%maxid_collision, maxval(nbody_system%pl%info(:)%collision_id)) + collider%maxid_collision = collider%maxid_collision + 1 + collider%collision_id = collider%maxid_collision + write(idstr,*) collider%collision_id + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "collision_id " // trim(adjustl(idstr))) + collider%impactors%regime = COLLRESOLVE_REGIME_MERGE + allocate(lmask, mold=pl%lmask) + lmask(:) = .false. + lmask(idx1(k)) = .true. + + allocate(plsub, mold=pl) + call pl%spill(plsub, lmask, ldestructive=.false.) + + ! Save the before snapshots + select type(before => collider%before) + class is (swiftest_nbody_system) + call move_alloc(plsub, before%pl) + end select + + call collision_history%take_snapshot(param,nbody_system, t, "particle") + + call impactors%dealloc() + end do + + ! Destroy the collision list now that the collisions are resolved + call pltp_collision%setup(0_I8B) + + end associate + end associate end select end select diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 9ce2be434..fab184281 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -756,8 +756,8 @@ end subroutine collision_util_shift_vector_to_origin module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) !! author: David A. Minton !! - !! Takes a minimal snapshot of the state of the nbody_system during an encounter so that the trajectories - !! can be played back through the encounter + !! Takes a minimal snapshot of the state of the nbody_system during a collision to record the before and after states of the + !! system through the collision. implicit none ! Internals class(collision_storage), intent(inout) :: self !! Swiftest storage object @@ -791,63 +791,70 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) snapshot%t = t case ("after") phase_val = 2 + case ("particle") + phase_val = -1 + allocate(collision_snapshot :: snapshot) + allocate(snapshot%collider, source=nbody_system%collider) case default - write(*,*) "collision_util_snapshot requies either 'before' or 'after' passed to 'arg'" + write(*,*) "collision_util_snapshot requies either 'before', 'after', or 'particle' passed to 'arg'" return end select - ! Get and record the energy of the system before the collision - call nbody_system%get_energy_and_momentum(param) - snapshot%collider%L_orbit(:,phase_val) = nbody_system%L_orbit(:) - snapshot%collider%L_spin(:,phase_val) = nbody_system%L_spin(:) - snapshot%collider%L_total(:,phase_val) = nbody_system%L_total(:) - snapshot%collider%ke_orbit(phase_val) = nbody_system%ke_orbit - snapshot%collider%ke_spin(phase_val) = nbody_system%ke_spin - snapshot%collider%pe(phase_val) = nbody_system%pe - snapshot%collider%be(phase_val) = nbody_system%be - snapshot%collider%te(phase_val) = nbody_system%te - - if (stage == "after") then - select type(before_snap => snapshot%collider%before ) - class is (swiftest_nbody_system) - select type(before_orig => nbody_system%collider%before) + if (stage /= "particle" ) then + ! Get and record the energy of the system before the collision + call nbody_system%get_energy_and_momentum(param) + snapshot%collider%L_orbit(:,phase_val) = nbody_system%L_orbit(:) + snapshot%collider%L_spin(:,phase_val) = nbody_system%L_spin(:) + snapshot%collider%L_total(:,phase_val) = nbody_system%L_total(:) + snapshot%collider%ke_orbit(phase_val) = nbody_system%ke_orbit + snapshot%collider%ke_spin(phase_val) = nbody_system%ke_spin + snapshot%collider%pe(phase_val) = nbody_system%pe + snapshot%collider%be(phase_val) = nbody_system%be + snapshot%collider%te(phase_val) = nbody_system%te + + if (stage == "after") then + select type(before_snap => snapshot%collider%before ) class is (swiftest_nbody_system) - select type(plsub => before_orig%pl) - class is (swiftest_pl) - ! Log the properties of the old and new bodies - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Removing bodies:") - do i = 1, plsub%nbody - write(message,*) trim(adjustl(plsub%info(i)%name)), " (", trim(adjustl(plsub%info(i)%particle_type)),")" - call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) - end do - - allocate(before_snap%pl, source=plsub) + select type(before_orig => nbody_system%collider%before) + class is (swiftest_nbody_system) + select type(plsub => before_orig%pl) + class is (swiftest_pl) + ! Log the properties of the old and new bodies + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Removing bodies:") + do i = 1, plsub%nbody + write(message,*) trim(adjustl(plsub%info(i)%name)), " (", trim(adjustl(plsub%info(i)%particle_type)),")" + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + end do + + allocate(before_snap%pl, source=plsub) + end select + deallocate(before_orig%pl) end select - deallocate(before_orig%pl) - end select - end select - - - select type(after_snap => snapshot%collider%after ) - class is (swiftest_nbody_system) - select type(after_orig => nbody_system%collider%after) - class is (swiftest_nbody_system) - select type(plnew => after_orig%pl) - class is (swiftest_pl) - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Adding bodies:") - do i = 1, plnew%nbody - write(message,*) trim(adjustl(plnew%info(i)%name)), " (", trim(adjustl(plnew%info(i)%particle_type)),")" - call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) - end do - call swiftest_io_log_one_message(COLLISION_LOG_OUT, & - "***********************************************************" // & - "***********************************************************") - allocate(after_snap%pl, source=plnew) end select - deallocate(after_orig%pl) - end select - end select + select type(after_snap => snapshot%collider%after ) + class is (swiftest_nbody_system) + select type(after_orig => nbody_system%collider%after) + class is (swiftest_nbody_system) + select type(plnew => after_orig%pl) + class is (swiftest_pl) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Adding bodies:") + do i = 1, plnew%nbody + write(message,*) trim(adjustl(plnew%info(i)%name)), " (", trim(adjustl(plnew%info(i)%particle_type)),")" + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + end do + call swiftest_io_log_one_message(COLLISION_LOG_OUT, & + "***********************************************************" // & + "***********************************************************") + allocate(after_snap%pl, source=plnew) + end select + deallocate(after_orig%pl) + end select + end select + end if + end if + + if ((stage == "after") .or. (stage == "particle")) then ! Save the snapshot for posterity call self%save(snapshot) deallocate(snapshot) diff --git a/src/swiftest/swiftest_discard.f90 b/src/swiftest/swiftest_discard.f90 index 2d7d2c7a2..311bf21f1 100644 --- a/src/swiftest/swiftest_discard.f90 +++ b/src/swiftest/swiftest_discard.f90 @@ -274,7 +274,7 @@ subroutine swiftest_discard_pl_tp(tp, nbody_system, param) write(idstri, *) tp%id(i) write(idstrj, *) pl%id(j) write(timestr, *) nbody_system%t - write(message, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & + write(message, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & // " too close to massive body " // trim(adjustl(pl%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & // " at t = " // trim(adjustl(timestr)) call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)