From 88610dc314fe555ae9f3bee991ae13366c830708 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 21 May 2021 15:57:56 -0400 Subject: [PATCH] Improved energy balancing due to collisions with central body by allowing radius to change and also to keep track of energy terms that can no longer be accounted for (such as potential energy of a body that has collided with the central body) --- src/io/io_conservation_report.f90 | 8 +----- src/main/swiftest_symba.f90 | 3 ++- src/modules/module_interfaces.f90 | 7 +++--- src/modules/swiftest_data_structures.f90 | 8 +++--- src/symba/symba_discard_conserve_mtm.f90 | 32 +++++++++++++++++------- src/symba/symba_rearray.f90 | 4 +-- 6 files changed, 37 insertions(+), 25 deletions(-) diff --git a/src/io/io_conservation_report.f90 b/src/io/io_conservation_report.f90 index 6e245e6f9..0463cb193 100644 --- a/src/io/io_conservation_report.f90 +++ b/src/io/io_conservation_report.f90 @@ -64,18 +64,12 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, write(*, egytermfmt) Lerror, Ecoll_error, Etotal_error, Merror if (Ecoll_error > 0.0_DP) then write(*,*) 'Something has gone wrong! Collisional energy should not be positive!' - write(*,*) 'Comparisons with last time: ' - write(*,*) 'ke_orb : ',ke_orbit - ke_orb_last, (ke_orbit - ke_orb_last) / abs(Eorbit_last) - write(*,*) 'ke_spin: ',ke_spin - ke_spin_last, (ke_spin - ke_spin_last) / abs(Eorbit_last) - write(*,*) 'pe : ',pe - pe_last, (pe - pe_last) / abs(Eorbit_last) - write(*,*) 'Etot : ',Eorbit - Eorbit_last, (Eorbit - Eorbit_last) / abs(Eorbit_last) - write(*,*) end if + end if ke_orb_last = ke_orbit ke_spin_last = ke_spin pe_last = pe - Ltot_last(:) = Ltot_now(:) Eorbit_last = Eorbit end associate return diff --git a/src/main/swiftest_symba.f90 b/src/main/swiftest_symba.f90 index 123b456a1..77c1f5609 100644 --- a/src/main/swiftest_symba.f90 +++ b/src/main/swiftest_symba.f90 @@ -146,6 +146,7 @@ program swiftest_symba ! Save initial mass and angular momentum of the central body symba_plA%helio%swiftest%Mcb_initial = symba_plA%helio%swiftest%mass(1) + symba_plA%helio%swiftest%Rcb_initial = symba_plA%helio%swiftest%radius(1) symba_plA%helio%swiftest%Lcb_initial(:) = symba_plA%helio%swiftest%Ip(3,1) * symba_plA%helio%swiftest%mass(1) * & symba_plA%helio%swiftest%radius(1)**2 * symba_plA%helio%swiftest%rot(:,1) @@ -189,7 +190,7 @@ program swiftest_symba if (param%lenergy) then call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Eorbit_after, Ltot) - Ecollision = Eorbit_before - Eorbit_after ! Energy change resulting in this collisional event Total running energy offset from collision in this step + Ecollision = Eorbit_before - Eorbit_after ! Energy change resulting in this collisional event Total running energy offset from collision in this step symba_plA%helio%swiftest%Ecollisions = symba_plA%helio%swiftest%Ecollisions + Ecollision end if !if (ntp > 0) call util_dist_index_pltp(nplm, ntp, symba_plA, symba_tpA) diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index 8450dd115..a5918240f 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -877,12 +877,13 @@ END SUBROUTINE symba_discard_sun_pl END INTERFACE INTERFACE - subroutine symba_discard_conserve_mtm(swiftest_plA, ipl, lescape) + subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape) use swiftest_globals use swiftest_data_structures implicit none - integer(I4B), intent(in) :: ipl + type(user_input_parameters), intent(inout) :: param type(swiftest_pl), intent(inout) :: swiftest_plA + integer(I4B), intent(in) :: ipl logical, intent(in) :: lescape end subroutine END INTERFACE @@ -1093,7 +1094,7 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, type(symba_merger), intent(inout) :: mergeadd_list logical(LGT), intent(in) :: ldiscard_pl, ldiscard_tp real(DP), intent(in) :: mtiny - type(user_input_parameters), intent(in) :: param + type(user_input_parameters), intent(inout) :: param logical, dimension(:), allocatable, intent(inout) :: discard_l_pl integer(I4B), dimension(:), allocatable, intent(inout) :: discard_stat_list end subroutine symba_rearray diff --git a/src/modules/swiftest_data_structures.f90 b/src/modules/swiftest_data_structures.f90 index 3625b5a04..46b30ae6a 100644 --- a/src/modules/swiftest_data_structures.f90 +++ b/src/modules/swiftest_data_structures.f90 @@ -44,11 +44,13 @@ module swiftest_data_structures real(DP), dimension(:,:), allocatable :: Ip !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). Principal axis rotation assumed. real(DP), dimension(:,:), allocatable :: rot !! Body rotation vector in inertial coordinate frame real(DP), dimension(NDIM) :: Lcb_initial !! Initial angular momentum of the central body - real(DP), dimension(NDIM) :: dLcb = (/0.0_DP, 0.0_DP, 0.0_DP/) !! Change in angular momentum of the central body - real(DP) :: Mcb_initial !! Initial mass and change in mass of the central body + real(DP), dimension(NDIM) :: dLcb = [0.0_DP, 0.0_DP, 0.0_DP] !! Change in angular momentum of the central body + real(DP), dimension(NDIM) :: Lescape = [0.0_DP, 0.0_DP, 0.0_DP] !! Angular momentum of bodies that escaped the system (used for bookeeping) + real(DP) :: Mcb_initial !! Initial mass of the central body real(DP) :: dMcb = 0.0_DP !! Change in mass of the central body - real(DP), dimension(NDIM) :: Lescape = (/0.0_DP, 0.0_DP, 0.0_DP/)!! Angular momentum of bodies that escaped the system (used for bookeeping) real(DP) :: Mescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping) + real(DP) :: Rcb_initial !! Initial radius of the central body + real(DP) :: dRcb = 0.0_DP!! Change in the radius of the central body real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions integer(I4B), dimension(:,:), allocatable :: k_plpl integer(I8B) :: num_plpl_comparisons diff --git a/src/symba/symba_discard_conserve_mtm.f90 b/src/symba/symba_discard_conserve_mtm.f90 index 68abb0f2d..90e6aa7c5 100644 --- a/src/symba/symba_discard_conserve_mtm.f90 +++ b/src/symba/symba_discard_conserve_mtm.f90 @@ -1,45 +1,59 @@ -subroutine symba_discard_conserve_mtm(swiftest_plA, ipl, lescape) +subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape) !! author: David A. Minton !! !! Conserves system momentum when a body is lost from the system or collides with central body use swiftest - use module_interfaces, EXCEPT_THIS_ONE => symba_discard_conserve_mtm + use module_interfaces, except_this_one => symba_discard_conserve_mtm implicit none - - integer(I4B), intent(in) :: ipl + ! Arguments type(swiftest_pl), intent(inout) :: swiftest_plA + type(user_input_parameters), intent(inout) :: param + integer(I4B), intent(in) :: ipl logical, intent(in) :: lescape - + ! Internals real(DP), dimension(NDIM) :: Lpl, Lcb, xcom, vcom + reaL(DP) :: pe, ke associate(npl => swiftest_plA%nbody, xb => swiftest_plA%xb, vb => swiftest_plA%vb, & - rot => swiftest_plA%rot, Ip => swiftest_plA%Ip, rad => swiftest_plA%radius, mass => swiftest_plA%mass, & + rot => swiftest_plA%rot, Ip => swiftest_plA%Ip, radius => swiftest_plA%radius, mass => swiftest_plA%mass, & statipl => swiftest_plA%status(ipl), xbpl => swiftest_plA%xb(:,ipl), xbcb => swiftest_plA%xb(:,1)) xcom(:) = (mass(ipl) * xb(:, ipl) + mass(1) * xb(:,1)) / (mass(1) + mass(ipl)) vcom(:) = (mass(ipl) * vb(:, ipl) + mass(1) * vb(:,1)) / (mass(1) + mass(ipl)) - call util_crossproduct(xb(:,ipl) - xcom(:),vb(:,ipl) - vcom(:), Lpl) - Lpl(:) = mass(ipl) * (Lpl(:) + rad(ipl)**2 * Ip(3,ipl) * rot(:, ipl)) + call util_crossproduct(xb(:,ipl) - xcom(:), vb(:,ipl) - vcom(:), Lpl) + Lpl(:) = mass(ipl) * (Lpl(:) + radius(ipl)**2 * Ip(3,ipl) * rot(:, ipl)) call util_crossproduct(xb(:, 1) - xcom(:), vb(:, 1) - vcom(:), Lcb) Lcb(:) = mass(1) * Lcb(:) ! Add planet mass to central body accumulator if (lescape) then + ! Add the potential and kinetic energy of the lost body to the records + pe = -mass(1) * mass(ipl) / norm2(xb(:, ipl) - xb(:, 1)) + ke = 0.5_DP * mass(ipl) * dot_product(vb(:, ipl), vb(:, ipl)) + swiftest_plA%Elost_bodies = swiftest_plA%Elost_bodies + pe + ke swiftest_plA%Mescape = swiftest_plA%Mescape + mass(ipl) else + ! Add the potential energy of the lost body to the records + pe = -mass(1) * mass(ipl) / norm2(xb(:, ipl) - xb(:, 1)) + ke = 0.0_DP + swiftest_plA%Elost_bodies = swiftest_plA%Elost_bodies + pe swiftest_plA%dMcb = swiftest_plA%dMcb + mass(ipl) + swiftest_plA%dRcb = swiftest_plA%dRcb + 1.0_DP / 3.0_DP * (radius(ipl) / radius(1))**3 - 2.0_DP / 9.0_DP * (radius(ipl) / radius(1))**6 ! Update mass of central body to be consistent with its total mass mass(1) = swiftest_plA%Mcb_initial + swiftest_plA%dMcb + radius(1) = swiftest_plA%Rcb_initial + swiftest_plA%dRcb + param%rmin = radius(1) end if + swiftest_plA%Ecollisions = swiftest_plA%Ecollisions - (ke + pe) ! Add planet angular momentum to central body accumulator swiftest_plA%dLcb(:) = Lpl(:) + Lcb(:) + swiftest_plA%dLcb(:) ! Update rotation of central body to by consistent with its angular momentum - swiftest_plA%rot(:,1) = (swiftest_plA%Lcb_initial(:) + swiftest_plA%dLcb(:)) / (Ip(3, 1) * mass(1) * rad(1)**2) + swiftest_plA%rot(:,1) = (swiftest_plA%Lcb_initial(:) + swiftest_plA%dLcb(:)) / (Ip(3, 1) * mass(1) * radius(1)**2) ! Update position and velocity of central body xb(:, 1) = xcom(:) diff --git a/src/symba/symba_rearray.f90 b/src/symba/symba_rearray.f90 index 47c31a35c..53ea421ab 100644 --- a/src/symba/symba_rearray.f90 +++ b/src/symba/symba_rearray.f90 @@ -22,7 +22,7 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, type(symba_merger), intent(inout) :: mergeadd_list logical, intent(in) :: ldiscard_pl, ldiscard_tp real(DP), intent(in) :: mtiny - type(user_input_parameters), intent(in) :: param + type(user_input_parameters), intent(inout) :: param logical, dimension(:), allocatable, intent(inout) :: discard_l_pl integer(I4B), dimension(:), allocatable, intent(inout) :: discard_stat_list @@ -50,7 +50,7 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, cycle end if ! Resolve the discard - call symba_discard_conserve_mtm(symba_plA%helio%swiftest, discard_index_list(i), lescape) + call symba_discard_conserve_mtm(param, symba_plA%helio%swiftest, discard_index_list(i), lescape) ! Flip the main status flag to the discard state end do symba_plA%helio%swiftest%status(discard_index_list(:)) = discard_stat_list(:)