From e28c01e3148d1c38452dfdd802c91e34f58f8944 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 23 Dec 2022 19:20:23 -0500 Subject: [PATCH] Fixed bugs that were preventing disruptions from working in the simple model --- src/collision/collision_generate.f90 | 42 +++++++++++++++------------- src/collision/collision_resolve.f90 | 37 ++++++++++++------------ 2 files changed, 41 insertions(+), 38 deletions(-) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 5898311a3..4e55b7f51 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -209,7 +209,7 @@ module subroutine collision_generate_simple(self, nbody_system, param, t) class is (swiftest_nbody_system) select type(pl => nbody_system%pl) class is (swiftest_pl) - associate(impactors => self%impactors, fragments => self%fragments, status => self%status) + associate(impactors => self%impactors, status => self%status) select case (impactors%regime) case (COLLRESOLVE_REGIME_HIT_AND_RUN) @@ -226,30 +226,34 @@ module subroutine collision_generate_simple(self, nbody_system, param, t) write(*,*) "Error in swiftest_collision, unrecognized collision regime" call util_exit(FAILURE) end select + call self%set_mass_dist(param) + call self%disrupt(nbody_system, param, t) dpe = self%pe(2) - self%pe(1) nbody_system%Ecollisions = nbody_system%Ecollisions - dpe nbody_system%Euntracked = nbody_system%Euntracked + dpe - ! Populate the list of new bodies - nfrag = fragments%nbody - write(message, *) nfrag - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") - select case(impactors%regime) - case(COLLRESOLVE_REGIME_DISRUPTION) - status = DISRUPTED - ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1)) - fragments%id(1) = pl%id(ibiggest) - fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] - param%maxid = fragments%id(nfrag) - case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - status = SUPERCATASTROPHIC - fragments%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] - param%maxid = fragments%id(nfrag) - end select + associate (fragments => self%fragments) + ! Populate the list of new bodies + nfrag = fragments%nbody + write(message, *) nfrag + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") + select case(impactors%regime) + case(COLLRESOLVE_REGIME_DISRUPTION) + status = DISRUPTED + ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1)) + fragments%id(1) = pl%id(ibiggest) + fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] + param%maxid = fragments%id(nfrag) + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + status = SUPERCATASTROPHIC + fragments%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] + param%maxid = fragments%id(nfrag) + end select - call collision_resolve_mergeaddsub(nbody_system, param, t, status) - end associate + call collision_resolve_mergeaddsub(nbody_system, param, t, status) + end associate + end associate end select end select return diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 9e91b1ecd..f29e6abce 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -507,22 +507,21 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) associate(plpl_encounter => self, plpl_collision => nbody_system%plpl_collision, & collision_history => nbody_system%collision_history, pl => nbody_system%pl, cb => nbody_system%cb, & collider => nbody_system%collider, fragments => nbody_system%collider%fragments, impactors => nbody_system%collider%impactors) - associate( idx1 => plpl_collision%index1, idx2 => plpl_collision%index2) + if (plpl_collision%nenc == 0) return ! No collisions to resolve - if (plpl_collision%nenc == 0) return ! No collisions to resolve + ! Make sure that the heliocentric and barycentric coordinates are consistent with each other + call pl%vb2vh(nbody_system%cb) + call pl%rh2rb(nbody_system%cb) - ! Make sure that the heliocentric and barycentric coordinates are consistent with each other - call pl%vb2vh(nbody_system%cb) - call pl%rh2rb(nbody_system%cb) - - ! Get the energy before the collision is resolved - if (param%lenergy) then - call nbody_system%get_energy_and_momentum(param) - Eorbit_before = nbody_system%te - end if + ! Get the energy before the collision is resolved + if (param%lenergy) then + call nbody_system%get_energy_and_momentum(param) + Eorbit_before = nbody_system%te + end if - do loop = 1, MAXCASCADE + do loop = 1, MAXCASCADE + associate( idx1 => plpl_collision%index1, idx2 => plpl_collision%index2) ncollisions = plpl_collision%nenc write(timestr,*) t call swiftest_io_log_one_message(COLLISION_LOG_OUT, "") @@ -572,15 +571,15 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) write(*,*) "Consider reducing the step size or changing the parameters in the collisional model to reduce the number of fragments." call util_exit(FAILURE) end if - end do + end associate + end do - if (param%lenergy) then - call nbody_system%get_energy_and_momentum(param) - Eorbit_after = nbody_system%te - nbody_system%Ecollisions = nbody_system%Ecollisions + (Eorbit_after - Eorbit_before) - end if + if (param%lenergy) then + call nbody_system%get_energy_and_momentum(param) + Eorbit_after = nbody_system%te + nbody_system%Ecollisions = nbody_system%Ecollisions + (Eorbit_after - Eorbit_before) + end if - end associate end associate end select end select