From 4fd9671628495ccf23922556893546ea65036c65 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 4 Feb 2023 14:37:22 -0500 Subject: [PATCH 01/10] Fixed issu that caused solar collisions to not conserve angular momentum when restarted. cb%dL was not being computed on restarted runs. --- src/fraggle/fraggle_module.f90 | 2 +- src/swiftest/swiftest_io.f90 | 8 +++++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index 01412154d..fda3ebe1b 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -44,7 +44,7 @@ end subroutine fraggle_generate_disrupt module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) implicit none - class(collision_fraggle), intent(inout) :: self !! Collision system object + class(collision_fraggle), intent(inout) :: self !! Collision system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions real(DP), intent(in) :: t !! Time of collision diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index af434b2c5..79daa0220 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -579,7 +579,7 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param) integer(I4B) :: itmax, idmax, tslot real(DP), dimension(:), allocatable :: vals real(DP), dimension(1) :: rtemp - real(DP), dimension(NDIM) :: rot0, Ip0 + real(DP), dimension(NDIM) :: rot0, Ip0, L real(DP) :: mass0 associate (nc => self%system_history%nc, cb => self%cb) @@ -628,6 +628,8 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param) rot0(:) = rot0(:) * DEG2RAD call netcdf_io_check( nf90_get_var(nc%id, nc%Ip_varid, Ip0, start=[1,1,tslot], count=[NDIM,1,1]), "netcdf_io_get_t0_values_system Ip_varid" ) cb%L0(:) = Ip0(3) * mass0 * cb%R0**2 * rot0(:) + L(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:) + cb%dL(:) = L(:) - cb%L0 end if ! Retrieve the current bookkeeping variables @@ -2400,8 +2402,8 @@ module subroutine swiftest_io_param_writer(self, unit, iotype, v_list, iostat, i end if call io_param_writer_one("FIRSTKICK",param%lfirstkick, unit) - if (param%GMTINY > 0.0_DP) call io_param_writer_one("GMTINY",param%GMTINY, unit) - if (param%min_GMfrag > 0.0_DP) call io_param_writer_one("MIN_GMFRAG",param%min_GMfrag, unit) + if (param%GMTINY >= 0.0_DP) call io_param_writer_one("GMTINY",param%GMTINY, unit) + if (param%min_GMfrag >= 0.0_DP) call io_param_writer_one("MIN_GMFRAG",param%min_GMfrag, unit) call io_param_writer_one("COLLISION_MODEL",param%collision_model, unit) if (param%collision_model == "FRAGGLE" ) then nseeds = size(param%seed) From 42fc8445347184aeab82a4a0eb7ca1455f1b8cdc Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 4 Feb 2023 20:19:54 -0500 Subject: [PATCH 02/10] Tweaks to angular momentum conservation --- src/fraggle/fraggle_generate.f90 | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index a5151d2d3..21ea232e6 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -469,11 +469,11 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu logical, intent(out) :: lfailure !! Did the velocity computation fail? ! Internals real(DP), parameter :: ENERGY_SUCCESS_METRIC = 1.0e-4_DP !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount) - real(DP) :: MOMENTUM_SUCCESS_METRIC = 10*epsilon(1.0_DP) !! Relative angular momentum error to accept as a success (should be *much* stricter than energy) + real(DP) :: MOMENTUM_SUCCESS_METRIC = 2*epsilon(1.0_DP) !! Relative angular momentum error to accept as a success (should be *much* stricter than energy) integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps, nsteps_best logical :: lhitandrun, lsupercat - real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, dL, drot, rot_new - real(DP) :: vimp, vmag, vesc, dE, E_residual, E_residual_best, E_residual_last, ke_min, ke_avail, ke_remove, dE_best, fscale, dE_metric, mfrag, dL_metric, dL_best, rn + real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, dL, drot, rot_new, dL_metric, dL_best + real(DP) :: vimp, vmag, vesc, dE, E_residual, E_residual_best, E_residual_last, ke_min, ke_avail, ke_remove, dE_best, fscale, dE_metric, mfrag, rn integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale real(DP), parameter :: L_ROT_VEL_RATIO = 0.2_DP ! Ratio of angular momentum to put into rotation relative to velocity shear of fragments @@ -526,7 +526,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu lfailure = .false. dE_metric = huge(1.0_DP) dE_best = huge(1.0_DP) - dL_best = huge(1.0_DP) + dL_best(:) = huge(1.0_DP) nsteps_best = 0 nsteps = 0 outer: do try = 1, MAXOUTER @@ -575,11 +575,11 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu angmtm: do j = 1, MAXANGMTM call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) - dL_metric = .mag.L_residual(:) / .mag.(collider_local%L_total(:,1)) + dL_metric(:) = abs(L_residual(:)) / .mag.(collider_local%L_total(:,1)) / MOMENTUM_SUCCESS_METRIC - if (dL_metric / MOMENTUM_SUCCESS_METRIC <= 1.0_DP) exit angmtm + if (all(dL_metric(:) <= 1.0_DP)) exit angmtm L_residual_unit(:) = .unit. L_residual(:) - do i = 1, fragments%nbody + do i = fragments%nbody,1,-1 mfrag = sum(fragments%mass(i:fragments%nbody)) drot(:) = -L_ROT_VEL_RATIO * L_residual(:) / (fragments%mtot * fragments%Ip(3,i) * fragments%radius(i)**2) rot_new(:) = fragments%rot(:,i) + drot(:) @@ -614,7 +614,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu call fragments%set_coordinate_system() call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") ke_avail = 0.0_DP - do i = 1, fragments%nbody + do i = fragments%nbody, 1, -1 ke_avail = ke_avail + 0.5_DP * fragments%mass(i) * max(fragments%vmag(i) - vesc,0.0_DP)**2 end do @@ -623,14 +623,14 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu E_residual = dE + impactors%Qloss L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) - dL_metric = .mag.L_residual(:) / .mag.collider_local%L_total(:,1) + dL_metric(:) = abs(L_residual(:)) / .mag.collider_local%L_total(:,1) / MOMENTUM_SUCCESS_METRIC ! Check if we've converged on our constraints - if (dL_metric < MOMENTUM_SUCCESS_METRIC) then + if (all(dL_metric(:) <= 1.0_DP)) then if ((abs(E_residual) < abs(E_residual_best)) .or. ((dE < 0.0_DP) .and. (dE_best >= 0.0_DP))) then ! This is our best case so far. Save it for posterity E_residual_best = E_residual dE_best = dE - dL_best = dL_metric + dL_best(:) = dL_metric(:) nsteps_best = nsteps do concurrent(i = 1:fragments%nbody) @@ -672,12 +672,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu vmin_guess = vmin_guess + delta_v vmax_guess = vmax_guess - delta_v end do outer - lfailure = (dE_best > 0.0_DP) .or. (dL_best > MOMENTUM_SUCCESS_METRIC) + lfailure = (dE_best > 0.0_DP) .or. (any(dL_best(:) > 1.0_DP)) write(message, *) nsteps if (lfailure) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle velocity calculation failed to converge after " // trim(adjustl(message)) // " steps. The best solution found had:") - write(message,*) "|dL|/|L0| = ",dL_best + write(message,*) "|dL|/|L0| = ",.mag.(dL_best) * MOMENTUM_SUCCESS_METRIC call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) write(message,*) " dE = ",dE_best call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) From cb756a470c030d3abd2d4983f42f077ed96208c0 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 4 Feb 2023 23:10:19 -0500 Subject: [PATCH 03/10] Relaxed the L criterion again --- src/fraggle/fraggle_generate.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 21ea232e6..d0afb7f7d 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -469,7 +469,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu logical, intent(out) :: lfailure !! Did the velocity computation fail? ! Internals real(DP), parameter :: ENERGY_SUCCESS_METRIC = 1.0e-4_DP !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount) - real(DP) :: MOMENTUM_SUCCESS_METRIC = 2*epsilon(1.0_DP) !! Relative angular momentum error to accept as a success (should be *much* stricter than energy) + real(DP) :: MOMENTUM_SUCCESS_METRIC = 10*epsilon(1.0_DP) !! Relative angular momentum error to accept as a success (should be *much* stricter than energy) integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps, nsteps_best logical :: lhitandrun, lsupercat real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, dL, drot, rot_new, dL_metric, dL_best From feabec73212c85cb41ad2ef20923460ee95af98b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 5 Feb 2023 19:27:29 -0500 Subject: [PATCH 04/10] Fixed potential issue wehre append subroutines might fail if the appended array has 0 size. --- src/swiftest/swiftest_util.f90 | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 05f9383ee..74a795856 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -31,6 +31,8 @@ module subroutine swiftest_util_append_arr_char_string(arr, source, nold, nsrc, if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(nold+nnew)) else @@ -59,6 +61,8 @@ module subroutine swiftest_util_append_arr_DP(arr, source, nold, nsrc, lsource_m if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(nold+nnew)) else @@ -87,6 +91,8 @@ module subroutine swiftest_util_append_arr_DPvec(arr, source, nold, nsrc, lsourc if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(NDIM,nold+nnew)) else @@ -117,6 +123,8 @@ module subroutine swiftest_util_append_arr_I4B(arr, source, nold, nsrc, lsource_ if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(nold+nnew)) else @@ -146,6 +154,8 @@ module subroutine swiftest_util_append_arr_info(arr, source, nold, nsrc, lsource if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(nold+nnew)) else @@ -178,6 +188,8 @@ module subroutine swiftest_util_append_arr_kin(arr, source, nold, nsrc, lsource_ if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(nold+nnew)) else @@ -206,6 +218,8 @@ module subroutine swiftest_util_append_arr_logical(arr, source, nold, nsrc, lsou if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(nold+nnew)) else From 4ae28d92610ab34a90adb00d5f10a4f46ad42b2d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 5 Feb 2023 19:28:17 -0500 Subject: [PATCH 05/10] New angular momentum constraint will add residual L to the center of mass velocity of the fragments. --- src/fraggle/fraggle_generate.f90 | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index d0afb7f7d..315c8d8b5 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -472,7 +472,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP) :: MOMENTUM_SUCCESS_METRIC = 10*epsilon(1.0_DP) !! Relative angular momentum error to accept as a success (should be *much* stricter than energy) integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps, nsteps_best logical :: lhitandrun, lsupercat - real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, dL, drot, rot_new, dL_metric, dL_best + real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, L_residual_best, dL, drot, rot_new, dL_metric real(DP) :: vimp, vmag, vesc, dE, E_residual, E_residual_best, E_residual_last, ke_min, ke_avail, ke_remove, dE_best, fscale, dE_metric, mfrag, rn integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale @@ -526,7 +526,6 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu lfailure = .false. dE_metric = huge(1.0_DP) dE_best = huge(1.0_DP) - dL_best(:) = huge(1.0_DP) nsteps_best = 0 nsteps = 0 outer: do try = 1, MAXOUTER @@ -613,6 +612,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) call fragments%set_coordinate_system() call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") + ke_avail = 0.0_DP do i = fragments%nbody, 1, -1 ke_avail = ke_avail + 0.5_DP * fragments%mass(i) * max(fragments%vmag(i) - vesc,0.0_DP)**2 @@ -629,14 +629,10 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu if (all(dL_metric(:) <= 1.0_DP)) then if ((abs(E_residual) < abs(E_residual_best)) .or. ((dE < 0.0_DP) .and. (dE_best >= 0.0_DP))) then ! This is our best case so far. Save it for posterity E_residual_best = E_residual + L_residual_best(:) = L_residual(:) dE_best = dE - dL_best(:) = dL_metric(:) nsteps_best = nsteps - do concurrent(i = 1:fragments%nbody) - fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) - end do - if (allocated(collider%fragments)) deallocate(collider%fragments) allocate(collider%fragments, source=fragments) dE_metric = abs(E_residual) / impactors%Qloss @@ -672,18 +668,24 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu vmin_guess = vmin_guess + delta_v vmax_guess = vmax_guess - delta_v end do outer - lfailure = (dE_best > 0.0_DP) .or. (any(dL_best(:) > 1.0_DP)) + lfailure = (dE_best > 0.0_DP) + + call fraggle_generate_velocity_torque(-L_residual_best(:), fragments%mtot, impactors%rbcom, impactors%vbcom) + + do concurrent(i = 1:collider%fragments%nbody) + collider%fragments%vb(:,i) = collider%fragments%vc(:,i) + impactors%vbcom(:) + end do write(message, *) nsteps if (lfailure) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle velocity calculation failed to converge after " // trim(adjustl(message)) // " steps. The best solution found had:") - write(message,*) "|dL|/|L0| = ",.mag.(dL_best) * MOMENTUM_SUCCESS_METRIC - call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) - write(message,*) " dE = ",dE_best - call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) else call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Fraggle velocity calculation converged after " // trim(adjustl(message)) // " steps.") end if + write(message,*) "dL/|L0| = ",(L_residual_best(:))/.mag.collider_local%L_total(:,1) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + write(message,*) "dE/Qloss = ",-dE_best / impactors%Qloss + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) write(message,*) nsteps_best call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Best solution came after " // trim(adjustl(message)) // " steps.") From 00cfd5909ee28996b0b97019628d27a2234266be Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 6 Feb 2023 07:10:46 -0500 Subject: [PATCH 06/10] Moved the velocity torque subroutine to the general collision_util submodule. Added velocity torque to merge case to improve angular momentum conservation --- src/collision/collision_generate.f90 | 5 ++++- src/collision/collision_module.f90 | 7 ++++++ src/collision/collision_util.f90 | 28 +++++++++++++++++++++++ src/fraggle/fraggle_generate.f90 | 33 ++-------------------------- 4 files changed, 41 insertions(+), 32 deletions(-) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 4ef6ec215..f0fc6d9e7 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -172,7 +172,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) real(DP), intent(in) :: t !! The time of the collision ! Internals integer(I4B) :: i, j, k, ibiggest - real(DP), dimension(NDIM) :: L_spin_new + real(DP), dimension(NDIM) :: L_spin_new, L_residual real(DP) :: volume character(len=STRMAX) :: message @@ -228,6 +228,9 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) ! Get the energy of the system after the collision call self%get_energy_and_momentum(nbody_system, param, phase="after") + L_residual(:) = (self%L_total(:,2) - self%L_total(:,1)) + call collision_util_velocity_torque(-L_residual(:), fragments%mass(1), fragments%rb(:,1), fragments%vb(:,1)) + ! Update any encounter lists that have the removed bodies in them so that they instead point to the new body do k = 1, nbody_system%plpl_encounter%nenc do j = 1, impactors%ncoll diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index d6323f522..2957dffaf 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -521,6 +521,13 @@ module subroutine collision_util_setup_fragments(self, n) integer(I4B), intent(in) :: n !! Number of particles to allocate space for end subroutine collision_util_setup_fragments + module subroutine collision_util_velocity_torque(dL, mass, r, v) + implicit none + real(DP), dimension(:), intent(in) :: dL !! Change in angular momentum to apply + real(DP), intent(in) :: mass !! Mass of body + real(DP), dimension(:), intent(in) :: r !! Position of body wrt system center of mass + real(DP), dimension(:), intent(inout) :: v !! Velocity of body wrt system center of mass + end subroutine collision_util_velocity_torque end interface contains diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 63d623de4..7adbe2ec5 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -1028,4 +1028,32 @@ module subroutine collision_util_set_original_scale_factors(self) end subroutine collision_util_set_original_scale_factors + module subroutine collision_util_velocity_torque(dL, mass, r, v) + !! author: David A. Minton + !! + !! Applies a torque to a body's center of mass velocity given a change in angular momentum + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: dL !! Change in angular momentum to apply + real(DP), intent(in) :: mass !! Mass of body + real(DP), dimension(:), intent(in) :: r !! Position of body wrt system center of mass + real(DP), dimension(:), intent(inout) :: v !! Velocity of body wrt system center of mass + ! Internals + real(DP), dimension(NDIM) :: dL_unit, r_unit, r_lever, vapply + real(DP) :: rmag, r_lever_mag + + dL_unit(:) = .unit. dL + r_unit(:) = .unit.r(:) + rmag = .mag.r(:) + ! Project the position vector onto the plane defined by the angular momentum vector and the origin to get the "lever arm" distance + r_lever(:) = dL_unit(:) .cross. (r(:) .cross. dL_unit(:)) + r_lever_mag = .mag.r_lever(:) + if ((r_lever_mag > epsilon(1.0_DP)) .and. (rmag > epsilon(1.0_DP))) then + vapply(:) = (dL(:) .cross. r(:)) / (mass * rmag**2) + v(:) = v(:) + vapply(:) + end if + + return + end subroutine collision_util_velocity_torque + end submodule s_collision_util \ No newline at end of file diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 315c8d8b5..8b805cacb 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -602,7 +602,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end if dL(:) = -L_residual(:) * fragments%mass(i) / fragments%mtot - drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 - call fraggle_generate_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i)) + call collision_util_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i)) call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) fragments%vmag(i) = .mag.fragments%vc(:,i) @@ -670,7 +670,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end do outer lfailure = (dE_best > 0.0_DP) - call fraggle_generate_velocity_torque(-L_residual_best(:), fragments%mtot, impactors%rbcom, impactors%vbcom) + call collision_util_velocity_torque(-L_residual_best(:), fragments%mtot, impactors%rbcom, impactors%vbcom) do concurrent(i = 1:collider%fragments%nbody) collider%fragments%vb(:,i) = collider%fragments%vc(:,i) + impactors%vbcom(:) @@ -695,33 +695,4 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end subroutine fraggle_generate_vel_vec - subroutine fraggle_generate_velocity_torque(dL, mass, r, v) - !! author: David A. Minton - !! - !! Applies a torque to a body's center of mass velocity given a change in angular momentum - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: dL !! Change in angular momentum to apply - real(DP), intent(in) :: mass !! Mass of body - real(DP), dimension(:), intent(in) :: r !! Position of body wrt system center of mass - real(DP), dimension(:), intent(inout) :: v !! Velocity of body wrt system center of mass - ! Internals - real(DP), dimension(NDIM) :: dL_unit, r_unit, r_lever, vapply - real(DP) :: rmag, r_lever_mag - - dL_unit(:) = .unit. dL - r_unit(:) = .unit.r(:) - rmag = .mag.r(:) - ! Project the position vector onto the plane defined by the angular momentum vector and the origin to get the "lever arm" distance - r_lever(:) = dL_unit(:) .cross. (r(:) .cross. dL_unit(:)) - r_lever_mag = .mag.r_lever(:) - if ((r_lever_mag > epsilon(1.0_DP)) .and. (rmag > epsilon(1.0_DP))) then - vapply(:) = (dL(:) .cross. r(:)) / (mass * rmag**2) - v(:) = v(:) + vapply(:) - end if - - return - end subroutine fraggle_generate_velocity_torque - - end submodule s_fraggle_generate From e0c8bac07016b7ea32d48d172fa94ddd49983f6d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 6 Feb 2023 13:06:23 -0500 Subject: [PATCH 07/10] More improvements to Fraggle aimed at getting angular momentum and success rates better --- src/fraggle/fraggle_generate.f90 | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 8b805cacb..35a323f6e 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -26,6 +26,7 @@ module subroutine fraggle_generate(self, nbody_system, param, t) real(DP), intent(in) :: t !! Time of collision ! Internals integer(I4B) :: i, ibiggest, nfrag + real(DP), dimension(NDIM) :: L_residual, vbcom_orig, dvb character(len=STRMAX) :: message logical :: lfailure @@ -37,6 +38,9 @@ module subroutine fraggle_generate(self, nbody_system, param, t) ! Set the coordinate system of the impactors call impactors%set_coordinate_system() + + vbcom_orig(:) = impactors%vbcom(:) + select case (impactors%regime) case (COLLRESOLVE_REGIME_HIT_AND_RUN) call self%hitandrun(nbody_system, param, t) @@ -84,9 +88,21 @@ module subroutine fraggle_generate(self, nbody_system, param, t) end if end if + ! Get the energy and momentum of the system before and after the collision + call self%get_energy_and_momentum(nbody_system, param, phase="before") + call self%get_energy_and_momentum(nbody_system, param, phase="after") + L_residual(:) = (self%L_total(:,2) - self%L_total(:,1)) + associate (fragments => self%fragments) - ! Populate the list of new bodies nfrag = fragments%nbody + + ! Put any residual angular momentum into orbital velocity + call collision_util_velocity_torque(-L_residual(:), fragments%mtot, impactors%rbcom(:), impactors%vbcom(:)) + dvb(:) = impactors%vbcom(:) - vbcom_orig(:) + do concurrent(i = 1:nfrag) + fragments%vb(:,i) = fragments%vb(:,i) + dvb(:) + end do + select case(impactors%regime) case(COLLRESOLVE_REGIME_DISRUPTION) status = DISRUPTED @@ -481,7 +497,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP) :: vmin_guess = 1.01_DP real(DP) :: vmax_guess real(DP) :: delta_v, GC - integer(I4B), parameter :: MAXINNER = 50 + integer(I4B), parameter :: MAXINNER = 100 integer(I4B), parameter :: MAXOUTER = 10 integer(I4B), parameter :: MAXANGMTM = 10000 class(collision_fraggle), allocatable :: collider_local From 7caf2f1f347df79ed5a121643a2f2ae4f8317213 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 6 Feb 2023 17:44:34 -0500 Subject: [PATCH 08/10] Fixed issues with failed runs not having enough variables to do the dL computation correctly --- src/fraggle/fraggle_generate.f90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 35a323f6e..9d3dac19c 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -58,7 +58,6 @@ module subroutine fraggle_generate(self, nbody_system, param, t) end select call self%set_mass_dist(param) call self%disrupt(nbody_system, param, t, lfailure) - if (lfailure) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find an energy-losing solution. Simplifying the collisional model.") @@ -68,7 +67,6 @@ module subroutine fraggle_generate(self, nbody_system, param, t) impactors%regime = COLLRESOLVE_REGIME_DISRUPTION call self%set_mass_dist(param) call self%disrupt(nbody_system, param, t, lfailure) - if (lfailure) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find an energy-losing solution. Treating this as a bounce.") call collision_util_bounce_one(impactors%rb(:,1),impactors%vb(:,1),impactors%rbcom(:),impactors%vbcom(:),impactors%radius(1)) @@ -81,10 +79,14 @@ module subroutine fraggle_generate(self, nbody_system, param, t) fragments%radius(1:2) = impactors%radius(1:2) fragments%rb(:,1:2) = impactors%rb(:,1:2) fragments%vb(:,1:2) = impactors%vb(:,1:2) + do concurrent(i = 1:2) + fragments%rc(:,i) = fragments%rb(:,i) - impactors%rbcom(:) + fragments%vc(:,i) = fragments%vb(:,i) - impactors%vbcom(:) + end do fragments%Ip(:,1:2) = impactors%Ip(:,1:2) fragments%rot(:,1:2) = impactors%rot(:,1:2) + fragments%mtot = sum(fragments%mass(1:2)) end associate - end if end if From 63a20eebe3da31ad08ce25e0acd568ad3433c396 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 7 Feb 2023 09:31:22 -0500 Subject: [PATCH 09/10] Improved the handling of reading a run where the initial condition and param file exists, but not the data file. --- python/swiftest/swiftest/simulation_class.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 6e80909c4..3eaa65ee4 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -342,7 +342,7 @@ def __init__(self,read_param: bool = False, read_data: bool = False, simdir: os. # If the file doesn't exist, flag it for now so we know to create it param_file_found = False if read_param or read_data: - if self.read_param(read_init_cond = read_data): + if self.read_param(read_init_cond = True): # We will add the parameter file to the kwarg list. This will keep the set_parameter method from # overriding everything with defaults when there are no arguments passed to Simulation() kwargs['param_file'] = self.param_file @@ -394,9 +394,9 @@ def _type_scrub(output_data): iloop = int((self.param['TSTART'] - self.param['T0']) / self.param['DT']) twidth = int(np.ceil(np.log10(self.param['TSTOP']/(self.param['DT'] * self.param['ISTEP_OUT'])))) pre_message = f"Time: {self.param['TSTART']:.{twidth}e} / {self.param['TSTOP']:.{twidth}e} {self.TU_name} " - post_message = f"npl: {self.data['npl'].values[0]} ntp: {self.data['ntp'].values[0]}" - if "nplm" in self.data: - post_message += f" nplm: {self.data['nplm'].values[0]}" + post_message = f"npl: {self.init_cond['npl'].values[0]} ntp: {self.init_cond['ntp'].values[0]}" + if "nplm" in self.init_cond: + post_message += f" nplm: {self.init_cond['nplm'].values[0]}" if self.param['ENERGY']: post_message += f" dL/L0: {0.0:.5e} dE/|E0|: {0.0:+.5e}" post_message += f" Wall time / step: {0.0:.5e} s" From 91df16f2b135bffa9c607bc133a65c8c99cf9c60 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 7 Feb 2023 11:24:56 -0500 Subject: [PATCH 10/10] Fixed bugs in the velocity guess values --- src/fraggle/fraggle_generate.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 9d3dac19c..c91676621 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -60,7 +60,6 @@ module subroutine fraggle_generate(self, nbody_system, param, t) call self%disrupt(nbody_system, param, t, lfailure) if (lfailure) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find an energy-losing solution. Simplifying the collisional model.") - impactors%mass_dist(1) = impactors%mass(1) impactors%mass_dist(2) = max(0.5_DP * impactors%mass(2), self%min_mfrag) impactors%mass_dist(3) = impactors%mass(2) - impactors%mass_dist(2) @@ -496,7 +495,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP), dimension(collider%fragments%nbody) :: vscale real(DP), parameter :: L_ROT_VEL_RATIO = 0.2_DP ! Ratio of angular momentum to put into rotation relative to velocity shear of fragments ! 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.01_DP + real(DP) :: vmin_guess real(DP) :: vmax_guess real(DP) :: delta_v, GC integer(I4B), parameter :: MAXINNER = 100 @@ -539,6 +538,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu vimp = .mag. (impactors%vc(:,2) - impactors%vc(:,1)) vmax_guess = 1.1_DP * vimp + vmin_guess = 1.001_DP * vesc E_residual_best = huge(1.0_DP) lfailure = .false. @@ -682,7 +682,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu collider_local%fail_scale = collider_local%fail_scale * 1.001_DP ! Bring the minimum and maximum velocities closer together - delta_v = 0.125_DP * (vmax_guess - vmin_guess) + delta_v = (vmax_guess - vmin_guess) / 16.0_DP vmin_guess = vmin_guess + delta_v vmax_guess = vmax_guess - delta_v end do outer