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(:)