From 8d406e11fe704477859d9fb0bf7923490b18da9b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 17 Dec 2022 10:38:33 -0500 Subject: [PATCH] Too soon. Not quite compiled yet. More cleanup --- src/collision/collision_util.f90 | 61 +++++++++- src/encounter/encounter_util.f90 | 195 ++++++++++++++++++++++++++++++ src/fraggle/fraggle_setup.f90 | 2 +- src/fraggle/fraggle_util.f90 | 2 +- src/modules/collision_classes.f90 | 8 +- 5 files changed, 258 insertions(+), 10 deletions(-) diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 819f9df76..3a1aea3fa 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -386,7 +386,6 @@ module subroutine collision_util_set_coordinate_system(self) return end subroutine collision_util_set_coordinate_system - subroutine collision_util_save_snapshot(collision_history, snapshot) !! author: David A. Minton !! @@ -432,4 +431,64 @@ subroutine collision_util_save_snapshot(collision_history, snapshot) end subroutine collision_util_save_snapshot + module subroutine collision_util_snapshot(self, param, system, t, arg) + !! author: David A. Minton + !! + !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories + !! can be played back through the encounter + implicit none + ! Internals + class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store + real(DP), intent(in), optional :: t !! Time of snapshot if different from system time + character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision. + ! Arguments + class(collision_snapshot), allocatable :: snapshot + type(symba_pl) :: pl + character(len=:), allocatable :: stage + + if (present(arg)) then + stage = arg + else + stage = "" + end if + + select type (system) + class is (symba_nbody_system) + + select case(stage) + case("before") + ! Saves the states of the bodies involved in the collision before the collision is resolved + associate (idx => system%collision_system%impactors%idx, ncoll => system%collision_system%impactors%ncoll) + call pl%setup(ncoll, param) + pl%id(:) = system%pl%id(idx(:)) + pl%Gmass(:) = system%pl%Gmass(idx(:)) + pl%radius(:) = system%pl%radius(idx(:)) + pl%rot(:,:) = system%pl%rot(:,idx(:)) + pl%Ip(:,:) = system%pl%Ip(:,idx(:)) + pl%rh(:,:) = system%pl%rh(:,idx(:)) + pl%vh(:,:) = system%pl%vh(:,idx(:)) + pl%info(:) = system%pl%info(idx(:)) + !end select + allocate(system%collision_system%before%pl, source=pl) + end associate + case("after") + allocate(collision_snapshot :: snapshot) + allocate(snapshot%collision_system, source=system%collision_system) + snapshot%t = t + select type(param) + class is (symba_parameters) + call collision_util_save_snapshot(param%collision_history,snapshot) + end select + case default + write(*,*) "collision_util_snapshot requies either 'before' or 'after' passed to 'arg'" + end select + + end select + + return + end subroutine collision_util_snapshot + + end submodule s_collision_util \ No newline at end of file diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index c4acc1df9..0297c769c 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -410,4 +410,199 @@ subroutine encounter_util_save_snapshot(encounter_history, snapshot) end subroutine encounter_util_save_snapshot + module subroutine encounter_util_snapshot(self, param, system, t, arg) + !! author: David A. Minton + !! + !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories + !! can be played back through the encounter + implicit none + ! Internals + class(encounter_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store + real(DP), intent(in), optional :: t !! Time of snapshot if different from system time + character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) + ! Arguments + class(encounter_snapshot), allocatable :: snapshot + integer(I4B) :: i, pi, pj, k, npl_snap, ntp_snap, iflag + real(DP), dimension(NDIM) :: rrel, vrel, rcom, vcom + real(DP) :: Gmtot, a, q, capm, tperi + real(DP), dimension(NDIM,2) :: rb,vb + + if (.not.present(t)) then + write(*,*) "encounter_util_snapshot_encounter requires `t` to be passed" + return + end if + + if (.not.present(arg)) then + write(*,*) "encounter_util_snapshot_encounter requires `arg` to be passed" + return + end if + + select type(param) + class is (symba_parameters) + select type (system) + class is (symba_nbody_system) + select type(pl => system%pl) + class is (symba_pl) + select type (tp => system%tp) + class is (symba_tp) + associate(npl => pl%nbody, ntp => tp%nbody) + if (npl + ntp == 0) return + allocate(encounter_snapshot :: snapshot) + allocate(snapshot%pl, mold=pl) + allocate(snapshot%tp, mold=tp) + snapshot%iloop = param%iloop + + select type(pl_snap => snapshot%pl) + class is (symba_pl) + select type(tp_snap => snapshot%tp) + class is (symba_tp) + + select case(arg) + case("trajectory") + snapshot%t = t + + npl_snap = npl + ntp_snap = ntp + + if (npl > 0) then + pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE .and. pl%levelg(1:npl) == system%irec + npl_snap = count(pl%lmask(1:npl)) + end if + if (ntp > 0) then + tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE .and. tp%levelg(1:ntp) == system%irec + ntp_snap = count(tp%lmask(1:ntp)) + end if + + if (npl_snap + ntp_snap == 0) return ! Nothing to snapshot + + pl_snap%nbody = npl_snap + + ! Take snapshot of the currently encountering massive bodies + if (npl_snap > 0) then + call pl_snap%setup(npl_snap, param) + pl_snap%levelg(:) = pack(pl%levelg(1:npl), pl%lmask(1:npl)) + pl_snap%id(:) = pack(pl%id(1:npl), pl%lmask(1:npl)) + pl_snap%info(:) = pack(pl%info(1:npl), pl%lmask(1:npl)) + pl_snap%Gmass(:) = pack(pl%Gmass(1:npl), pl%lmask(1:npl)) + do i = 1, NDIM + pl_snap%rh(i,:) = pack(pl%rh(i,1:npl), pl%lmask(1:npl)) + pl_snap%vh(i,:) = pack(pl%vb(i,1:npl), pl%lmask(1:npl)) + end do + if (param%lclose) then + pl_snap%radius(:) = pack(pl%radius(1:npl), pl%lmask(1:npl)) + end if + + if (param%lrotation) then + do i = 1, NDIM + pl_snap%Ip(i,:) = pack(pl%Ip(i,1:npl), pl%lmask(1:npl)) + pl_snap%rot(i,:) = pack(pl%rot(i,1:npl), pl%lmask(1:npl)) + end do + end if + call pl_snap%sort("id", ascending=.true.) + end if + + ! Take snapshot of the currently encountering test particles + tp_snap%nbody = ntp_snap + if (ntp_snap > 0) then + call tp_snap%setup(ntp_snap, param) + tp_snap%id(:) = pack(tp%id(1:ntp), tp%lmask(1:ntp)) + tp_snap%info(:) = pack(tp%info(1:ntp), tp%lmask(1:ntp)) + do i = 1, NDIM + tp_snap%rh(i,:) = pack(tp%rh(i,1:ntp), tp%lmask(1:ntp)) + tp_snap%vh(i,:) = pack(tp%vh(i,1:ntp), tp%lmask(1:ntp)) + end do + end if + + ! Save the snapshot + param%encounter_history%nid = param%encounter_history%nid + ntp_snap + npl_snap + call encounter_util_save_encounter(param%encounter_history,snapshot) + case("closest") + associate(plplenc_list => system%plplenc_list, pltpenc_list => system%pltpenc_list) + if (any(plplenc_list%lclosest(:))) then + call pl_snap%setup(2, param) + do k = 1, plplenc_list%nenc + if (plplenc_list%lclosest(k)) then + pi = plplenc_list%index1(k) + pj = plplenc_list%index2(k) + pl_snap%levelg(:) = pl%levelg([pi,pj]) + pl_snap%id(:) = pl%id([pi,pj]) + pl_snap%info(:) = pl%info([pi,pj]) + pl_snap%Gmass(:) = pl%Gmass([pi,pj]) + Gmtot = sum(pl_snap%Gmass(:)) + if (param%lclose) pl_snap%radius(:) = pl%radius([pi,pj]) + if (param%lrotation) then + do i = 1, NDIM + pl_snap%Ip(i,:) = pl%Ip(i,[pi,pj]) + pl_snap%rot(i,:) = pl%rot(i,[pi,pj]) + end do + end if + + ! Compute pericenter passage time to get the closest approach parameters + rrel(:) = plplenc_list%r2(:,k) - plplenc_list%r1(:,k) + vrel(:) = plplenc_list%v2(:,k) - plplenc_list%v1(:,k) + call orbel_xv2aqt(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi) + snapshot%t = t + tperi + if ((snapshot%t < maxval(pl_snap%info(:)%origin_time)) .or. & + (snapshot%t > minval(pl_snap%info(:)%discard_time))) cycle + + ! Computer the center mass of the pair + rcom(:) = (plplenc_list%r1(:,k) * pl_snap%Gmass(1) + plplenc_list%r2(:,k) * pl_snap%Gmass(2)) / Gmtot + vcom(:) = (plplenc_list%v1(:,k) * pl_snap%Gmass(1) + plplenc_list%v2(:,k) * pl_snap%Gmass(2)) / Gmtot + rb(:,1) = plplenc_list%r1(:,k) - rcom(:) + rb(:,2) = plplenc_list%r2(:,k) - rcom(:) + vb(:,1) = plplenc_list%v1(:,k) - vcom(:) + vb(:,2) = plplenc_list%v2(:,k) - vcom(:) + + ! Drift the relative orbit to get the new relative position and velocity + call drift_one(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), tperi, iflag) + if (iflag /= 0) write(*,*) "Danby error in encounter_util_snapshot_encounter. Closest approach positions and vectors may not be accurate." + + ! Get the new position and velocity vectors + rb(:,1) = -(pl_snap%Gmass(2) / Gmtot) * rrel(:) + rb(:,2) = (pl_snap%Gmass(1)) / Gmtot * rrel(:) + + vb(:,1) = -(pl_snap%Gmass(2) / Gmtot) * vrel(:) + vb(:,2) = (pl_snap%Gmass(1)) / Gmtot * vrel(:) + + ! Move the CoM assuming constant velocity over the time it takes to reach periapsis + rcom(:) = rcom(:) + vcom(:) * tperi + + ! Compute the heliocentric position and velocity vector at periapsis + pl_snap%rh(:,1) = rb(:,1) + rcom(:) + pl_snap%rh(:,2) = rb(:,2) + rcom(:) + pl_snap%vh(:,1) = vb(:,1) + vcom(:) + pl_snap%vh(:,2) = vb(:,2) + vcom(:) + + call pl_snap%sort("id", ascending=.true.) + call encounter_util_save_encounter(param%encounter_history,snapshot) + end if + end do + + plplenc_list%lclosest(:) = .false. + end if + + if (any(pltpenc_list%lclosest(:))) then + do k = 1, pltpenc_list%nenc + end do + pltpenc_list%lclosest(:) = .false. + end if + end associate + case default + write(*,*) "encounter_util_snapshot_encounter requires `arg` to be either `trajectory` or `closest`" + end select + end select + end select + end associate + end select + end select + end select + end select + + return + end subroutine encounter_util_snapshot + + + end submodule s_encounter_util \ No newline at end of file diff --git a/src/fraggle/fraggle_setup.f90 b/src/fraggle/fraggle_setup.f90 index 7b6551d88..5b9442ff6 100644 --- a/src/fraggle/fraggle_setup.f90 +++ b/src/fraggle/fraggle_setup.f90 @@ -45,7 +45,7 @@ module subroutine fraggle_setup_fragments(self, n, param) integer(I4B), intent(in) :: n class(swiftest_parameters), intent(in) :: param - call collision_util_setup_fragments(n, param) + call self%collision_fragments%setup(n, param) if (n < 0) return if (allocated(self%v_r_mag)) deallocate(self%v_r_mag) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index c1e2fccb5..ae9cb8854 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -139,7 +139,7 @@ module subroutine fraggle_util_dealloc_fragments(self) ! Arguments class(fraggle_fragments), intent(inout) :: self - call collision_util_deallocate_fragments(self) + call collision_util_dealloc_fragments(self) if (allocated(self%v_r_mag)) deallocate(self%v_r_mag) if (allocated(self%v_t_mag)) deallocate(self%v_t_mag) diff --git a/src/modules/collision_classes.f90 b/src/modules/collision_classes.f90 index 5c31ab592..635a3fb33 100644 --- a/src/modules/collision_classes.f90 +++ b/src/modules/collision_classes.f90 @@ -81,7 +81,6 @@ module collision_classes procedure :: accel => collision_util_placeholder_accel !! Placeholder subroutine to fulfill requirement for an accel method procedure :: kick => collision_util_placeholder_kick !! Placeholder subroutine to fulfill requirement for a kick method procedure :: step => collision_util_placeholder_step !! Placeholder subroutine to fulfill requirement for a step method - procedure :: set_coordinate_system => collision_set_coordinate_fragments !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. procedure :: setup => collision_setup_fragments !! Allocates arrays for n fragments in a Fraggle system. Passing n = 0 deallocates all arrays. procedure :: dealloc => collision_util_dealloc_fragments !! Deallocates all allocatable arrays end type collision_fragments @@ -152,7 +151,7 @@ end subroutine abstract_set_mass_dist end type collision_io_parameters type, extends(encounter_snapshot) :: collision_snapshot - logical :: lcollision !! Indicates that this snapshot contains at least one collision + logical :: lcollision !! Indicates that this snapshot contains at least one collision class(collision_system), allocatable :: collision_system !! impactors object at this snapshot contains procedure :: write_frame => collision_io_write_frame_snapshot !! Writes a frame of encounter data to file @@ -233,11 +232,6 @@ module subroutine collision_set_coordinate_impactors(self) class(collision_impactors), intent(inout) :: self !! Collider system object end subroutine collision_set_coordinate_impactors - module subroutine collision_set_coordinate_fragments(self) - implicit none - class(collision_fragments), intent(inout) :: self !! Fragment system object - end subroutine collision_set_coordinate_fragments - module subroutine collision_util_set_coordinate_system(self) implicit none class(collision_system), intent(inout) :: self !! Collisional system