From 9a3723fd0902a69133daf05728da8569d462ff3d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 23 Dec 2022 09:00:44 -0500 Subject: [PATCH] cleaned up some missing bits from the append utility --- src/collision/collision_generate.f90 | 5 +- src/collision/collision_regime.f90 | 13 --- src/collision/collision_util.f90 | 10 ++- src/fraggle/fraggle_generate.f90 | 124 +++++++++++++++++---------- src/misc/minimizer_module.f90 | 3 +- src/swiftest/swiftest_util.f90 | 15 +++- 6 files changed, 108 insertions(+), 62 deletions(-) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 2d4a006a8..2be2408bb 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -149,12 +149,15 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) select type(pl => nbody_system%pl) class is (swiftest_pl) + ! Get coordinate system + call self%set_coordinate_system() + ! Generate the merged body as a single fragment call self%setup_fragments(1) ! Calculate the initial energy of the nbody_system without the collisional family call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) - + ! The new body's metadata will be taken from the largest of the two impactor bodies, so we need ! its index in the main pl structure ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1)) diff --git a/src/collision/collision_regime.f90 b/src/collision/collision_regime.f90 index 197a84f42..a1ae47fe0 100644 --- a/src/collision/collision_regime.f90 +++ b/src/collision/collision_regime.f90 @@ -32,10 +32,6 @@ module subroutine collision_regime_impactors(self, nbody_system, param) select type(param) class is (swiftest_parameters) - mtot = sum(impactors%mass(:)) - impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / mtot - impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / mtot - select case(param%collision_model) case("MERGE") impactors%regime = COLLRESOLVE_REGIME_MERGE @@ -105,15 +101,6 @@ subroutine collision_regime_LS12(impactors, nbody_system, param) impactors%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) impactors%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) - ! Find the center of mass of the collisional system - mtot = sum(impactors%mass(:)) - impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / mtot - impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / mtot - - ! Find the point of impact between the two bodies - runit(:) = impactors%rb(:,2) - impactors%rb(:,1) - runit(:) = runit(:) / (.mag. runit(:)) - impactors%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * runit(:) ! Convert quantities back to the system units and save them into the fragment system impactors%mass_dist(:) = (impactors%mass_dist(:) / param%MU2KG) diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 784802dba..35d4e10f2 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -387,7 +387,7 @@ module subroutine collision_util_set_coordinate_system(self) ! Internals integer(I4B) :: i real(DP), dimension(NDIM) :: delta_r, delta_v, Ltot - real(DP) :: L_mag + real(DP) :: L_mag, mtot real(DP), dimension(NDIM, self%fragments%nbody) :: L_sigma associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody) @@ -413,6 +413,14 @@ module subroutine collision_util_set_coordinate_system(self) impactors%x_unit(:) = impactors%y_unit(:) .cross. impactors%z_unit(:) impactors%v_unit(:) = .unit.delta_v(:) + ! Find the center of mass of the collisional system + mtot = sum(impactors%mass(:)) + impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / mtot + impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / mtot + + ! Find the point of impact between the two bodies + impactors%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * impactors%y_unit(:) + ! The "bounce" unit vector is the projection of impactors%vbimp(:) = dot_product(delta_v(:),impactors%x_unit(:)) * impactors%x_unit(:) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 8e52da0c5..3cfb3ee82 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -118,7 +118,7 @@ module subroutine fraggle_generate_hitandrun(collider, nbody_system, param, t) class is (swiftest_nbody_system) select type(after => collider%after) class is (swiftest_nbody_system) - associate(impactors => collider%impactors, fragments => collider%fragments, pl => nbody_system%pl) + associate(impactors => collider%impactors,pl => nbody_system%pl) message = "Hit and run between" call collision_io_collider_message(nbody_system%pl, impactors%id, message) call swiftest_io_log_one_message(COLLISION_LOG_OUT, trim(adjustl(message))) @@ -150,7 +150,7 @@ module subroutine fraggle_generate_hitandrun(collider, nbody_system, param, t) call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Should have been a pure hit and run instead") nfrag = 0 else - nfrag = fragments%nbody + nfrag = collider%fragments%nbody write(message, *) nfrag call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") end if @@ -163,9 +163,9 @@ module subroutine fraggle_generate_hitandrun(collider, nbody_system, param, t) allocate(after%pl, source=before%pl) ! Be sure to save the pl so that snapshots still work else 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) + collider%fragments%id(1) = pl%id(ibiggest) + collider%fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] + param%maxid = collider%fragments%id(nfrag) status = HIT_AND_RUN_DISRUPT call collision_resolve_mergeaddsub(nbody_system, param, t, status) end if @@ -293,44 +293,46 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.) call collider%set_budgets() - call fraggle_generate_tan_vel(collider, lfailure) - if (lfailure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find tangential velocities") - cycle - end if - - call fraggle_generate_rad_vel(collider, lfailure) - if (lfailure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find radial velocities") - cycle - end if - - call fraggle_generate_spins(collider, f_spin, lfailure) - if (lfailure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find spins") - cycle - end if - - call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.) - dEtot = collider%Etot(2) - collider%Etot(1) - dLmag = .mag. (collider%Ltot(:,2) - collider%Ltot(:,1)) - exit - - lfailure = ((abs(dEtot + impactors%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) - if (lfailure) then - write(message, *) dEtot, abs(dEtot + impactors%Qloss) / FRAGGLE_ETOL - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high energy error: " // & - trim(adjustl(message))) - cycle - end if - - lfailure = ((abs(dLmag) / (.mag.collider%Ltot(:,1))) > FRAGGLE_LTOL) - if (lfailure) then - write(message,*) dLmag / (.mag.collider%Ltot(:,1)) - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high angular momentum error: " // & - trim(adjustl(message))) - cycle - end if + call fraggle_generate_vel_vec_guess(collider) + + ! call fraggle_generate_tan_vel(collider, lfailure) + ! if (lfailure) then + ! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find tangential velocities") + ! cycle + ! end if + + ! call fraggle_generate_rad_vel(collider, lfailure) + ! if (lfailure) then + ! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find radial velocities") + ! cycle + ! end if + + ! call fraggle_generate_spins(collider, f_spin, lfailure) + ! if (lfailure) then + ! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find spins") + ! cycle + ! end if + + ! call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.) + ! dEtot = collider%Etot(2) - collider%Etot(1) + ! dLmag = .mag. (collider%Ltot(:,2) - collider%Ltot(:,1)) + ! exit + + ! lfailure = ((abs(dEtot + impactors%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) + ! if (lfailure) then + ! write(message, *) dEtot, abs(dEtot + impactors%Qloss) / FRAGGLE_ETOL + ! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high energy error: " // & + ! trim(adjustl(message))) + ! cycle + ! end if + + ! lfailure = ((abs(dLmag) / (.mag.collider%Ltot(:,1))) > FRAGGLE_LTOL) + ! if (lfailure) then + ! write(message,*) dLmag / (.mag.collider%Ltot(:,1)) + ! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high angular momentum error: " // & + ! trim(adjustl(message))) + ! cycle + ! end if ! Check if any of the usual floating point exceptions happened, and fail the try if so call ieee_get_flag(ieee_usual, fpe_flag) @@ -442,6 +444,42 @@ subroutine fraggle_generate_pos_vec(collider, r_max_start) end subroutine fraggle_generate_pos_vec + subroutine fraggle_generate_vel_vec_guess(collider) + !! Author: David A. Minton + !! + !! Generates an initial "guess" for the velocitity vectors + implicit none + ! Arguments + class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object + ! Internals + integer(I4B) :: i, j + logical :: lcat + real(DP), dimension(NDIM) :: vimp_unit + real(DP), dimension(NDIM,collider%fragments%nbody) :: vnoise + real(DP), parameter :: VNOISE_MAG = 0.10_DP + real(DP) :: vmag + + associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) + lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + + ! "Bounce" the first two bodies + do i = 1,2 + fragments%vc(:,i) = impactors%vb(:,i) - impactors%vbcom(:) - 2 * impactors%vbimp(:) + end do + + vmag = .mag.impactors%vbcom(:) + vimp_unit(:) = .unit. impactors%vbimp(:) + call random_number(vnoise) + vnoise = (2 * vnoise - 1.0_DP) * vmag + do i = 3, nfrag + vimp_unit(:) = .unit. (fragments%rc(:,i) + impactors%rbcom(:) - impactors%rbimp(:)) + fragments%vc(:,i) = vmag * vimp_unit(:) + vnoise(:,i) + end do + end associate + return + end subroutine fraggle_generate_vel_vec_guess + + subroutine fraggle_generate_spins(collider, f_spin, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! diff --git a/src/misc/minimizer_module.f90 b/src/misc/minimizer_module.f90 index 5a189ad03..772db418a 100644 --- a/src/misc/minimizer_module.f90 +++ b/src/misc/minimizer_module.f90 @@ -140,9 +140,10 @@ module subroutine minimize_bfgs(f, N, x0, eps, maxloop, lerr, x1) call ieee_get_status(original_fpe_status) ! Save the original floating point exception status call ieee_set_flag(ieee_all, .false.) ! Set all flags to quiet - allocate(fpe_flag(size(ieee_usual))) + if (.not.allocated(fpe_flag)) allocate(fpe_flag(size(ieee_usual))) lerr = .false. + if (allocated(x1)) deallocate(x1) allocate(x1, source=x0) ! Initialize approximate Hessian with the identity matrix (i.e. begin with method of steepest descent) ! Get initial gradient and initialize arrays for updated values of gradient and x diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 3cd3793e1..3d6cd2047 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -235,11 +235,13 @@ module subroutine swiftest_util_append_body(self, source, lsource_mask) nsrc = source%nbody nnew = count(lsource_mask(1:nsrc)) - call swiftest_util_append(self%info, source%info, nold, nsrc, lsource_mask) call swiftest_util_append(self%id, source%id, nold, nsrc, lsource_mask) + call swiftest_util_append(self%info, source%info, nold, nsrc, lsource_mask) + call swiftest_util_append(self%lmask, source%lmask, nold, nsrc, lsource_mask) call swiftest_util_append(self%status, source%status, nold, nsrc, lsource_mask) call swiftest_util_append(self%ldiscard, source%ldiscard, nold, nsrc, lsource_mask) - call swiftest_util_append(self%lmask, source%lmask, nold, nsrc, lsource_mask) + call swiftest_util_append(self%lencounter, source%lencounter, nold, nsrc, lsource_mask) + call swiftest_util_append(self%lcollision, source%lcollision, nold, nsrc, lsource_mask) call swiftest_util_append(self%mu, source%mu, nold, nsrc, lsource_mask) call swiftest_util_append(self%rh, source%rh, nold, nsrc, lsource_mask) call swiftest_util_append(self%vh, source%vh, nold, nsrc, lsource_mask) @@ -250,6 +252,9 @@ module subroutine swiftest_util_append_body(self, source, lsource_mask) call swiftest_util_append(self%atide, source%atide, nold, nsrc, lsource_mask) call swiftest_util_append(self%agr, source%agr, nold, nsrc, lsource_mask) call swiftest_util_append(self%ir3h, source%ir3h, nold, nsrc, lsource_mask) + call swiftest_util_append(self%isperi, source%isperi, nold, nsrc, lsource_mask) + call swiftest_util_append(self%peri, source%peri, nold, nsrc, lsource_mask) + call swiftest_util_append(self%atp, source%atp, nold, nsrc, lsource_mask) call swiftest_util_append(self%a, source%a, nold, nsrc, lsource_mask) call swiftest_util_append(self%e, source%e, nold, nsrc, lsource_mask) call swiftest_util_append(self%inc, source%inc, nold, nsrc, lsource_mask) @@ -282,15 +287,19 @@ module subroutine swiftest_util_append_pl(self, source, lsource_mask) call swiftest_util_append(self%rhill, source%rhill, nold, nsrc, lsource_mask) call swiftest_util_append(self%renc, source%renc, nold, nsrc, lsource_mask) call swiftest_util_append(self%radius, source%radius, nold, nsrc, lsource_mask) + call swiftest_util_append(self%density, source%density, nold, nsrc, lsource_mask) call swiftest_util_append(self%rbeg, source%rbeg, nold, nsrc, lsource_mask) call swiftest_util_append(self%rend, source%rend, nold, nsrc, lsource_mask) call swiftest_util_append(self%vbeg, source%vbeg, nold, nsrc, lsource_mask) - call swiftest_util_append(self%density, source%density, nold, nsrc, lsource_mask) call swiftest_util_append(self%Ip, source%Ip, nold, nsrc, lsource_mask) call swiftest_util_append(self%rot, source%rot, nold, nsrc, lsource_mask) call swiftest_util_append(self%k2, source%k2, nold, nsrc, lsource_mask) call swiftest_util_append(self%Q, source%Q, nold, nsrc, lsource_mask) call swiftest_util_append(self%tlag, source%tlag, nold, nsrc, lsource_mask) + call swiftest_util_append(self%kin, source%kin, nold, nsrc, lsource_mask) + call swiftest_util_append(self%lmtiny, source%lmtiny, nold, nsrc, lsource_mask) + call swiftest_util_append(self%nplenc, source%nplenc, nold, nsrc, lsource_mask) + call swiftest_util_append(self%ntpenc, source%ntpenc, nold, nsrc, lsource_mask) if (allocated(self%k_plpl)) deallocate(self%k_plpl)