From e31cda21633de4ce8aea0a539882264571d9fdfd Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 20 Aug 2021 10:49:12 -0400 Subject: [PATCH] Cleaned up fragmentation code and fixed issues related to adding sorting to the index method in the calculation of the system energy --- src/fragmentation/fragmentation.f90 | 72 ++++++++++++++++------------- 1 file changed, 40 insertions(+), 32 deletions(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 501ed050a..f020fa36f 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -100,7 +100,6 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, call reset_fragments() call define_coordinate_system() - call construct_temporary_system() ! Calculate the initial energy of the system without the collisional family call calculate_system_energy(linclude_fragments=.false.) @@ -150,25 +149,25 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, end do call restore_scale_factors() - write(*, "(' -------------------------------------------------------------------------------------')") - write(*, "(' Final diagnostic')") - write(*, "(' -------------------------------------------------------------------------------------')") - call calculate_system_energy(linclude_fragments=.true.) - if (lfailure) then - write(*,*) "symba_frag_pos failed after: ",try," tries" - do ii = 1, nfrag - vb_frag(:, ii) = vcom(:) - end do - else - write(*,*) "symba_frag_pos succeeded after: ",try," tries" - write(*, "(' dL_tot should be very small' )") - write(*,fmtlabel) ' dL_tot |', dLmag / Lmag_before - write(*, "(' dE_tot should be negative and equal to Qloss' )") - write(*,fmtlabel) ' dE_tot |', dEtot / abs(Etot_before) - write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before) - write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before) - end if - write(*, "(' -------------------------------------------------------------------------------------')") + ! write(*, "(' -------------------------------------------------------------------------------------')") + ! write(*, "(' Final diagnostic')") + ! write(*, "(' -------------------------------------------------------------------------------------')") + ! call calculate_system_energy(linclude_fragments=.true.) + ! if (lfailure) then + ! write(*,*) "symba_frag_pos failed after: ",try," tries" + ! do ii = 1, nfrag + ! vb_frag(:, ii) = vcom(:) + ! end do + ! else + ! write(*,*) "symba_frag_pos succeeded after: ",try," tries" + ! write(*, "(' dL_tot should be very small' )") + ! write(*,fmtlabel) ' dL_tot |', dLmag / Lmag_before + ! write(*, "(' dE_tot should be negative and equal to Qloss' )") + ! write(*,fmtlabel) ' dE_tot |', dEtot / abs(Etot_before) + ! write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before) + ! write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before) + ! end if + ! write(*, "(' -------------------------------------------------------------------------------------')") call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily @@ -355,7 +354,7 @@ subroutine construct_temporary_system() associate(pl => system%pl, npl => system%pl%nbody, cb => system%cb) if (size(lexclude) /= npl + nfrag) then - allocate(lexclude_tmp(npl_new)) + allocate(lexclude_tmp(npl + nfrag)) lexclude_tmp(1:npl) = lexclude(1:npl) call move_alloc(lexclude_tmp, lexclude) end if @@ -364,8 +363,9 @@ subroutine construct_temporary_system() elsewhere lexclude(1:npl) = .false. end where - lexclude(npl+1:npl_new) = .true. + lexclude(npl+1:(npl + nfrag)) = .true. if (allocated(tmpparam)) deallocate(tmpparam) + if (allocated(tmpsys)) deallocate(tmpsys) allocate(tmpparam, source=param) call setup_construct_system(tmpsys, param) call tmpsys%tp%setup(0, param) @@ -388,8 +388,11 @@ subroutine add_fragments_to_tmpsys() implicit none ! Internals integer(I4B) :: i + class(swiftest_pl), allocatable :: pl_discards + logical, dimension(:), allocatable :: lexclude_tmp associate(pl => system%pl, npl => system%pl%nbody) + npl_new = npl + nfrag tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag) tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * tmpparam%GU @@ -407,11 +410,20 @@ subroutine add_fragments_to_tmpsys() ! Disable the collisional family for subsequent energy calculations and coordinate shifts lexclude(family(:)) = .true. lexclude(npl+1:npl_new) = .false. - where(lexclude(:)) - tmpsys%pl%status(:) = INACTIVE + where(lexclude(1:npl_new)) + tmpsys%pl%status(1:npl_new) = INACTIVE elsewhere - tmpsys%pl%status(:) = ACTIVE + tmpsys%pl%status(1:npl_new) = ACTIVE end where + allocate(pl_discards, mold=tmpsys%pl) + call tmpsys%pl%spill(pl_discards, lspill_list=lexclude(1:npl_new), ldestructive=.true.) + npl_new = count(.not.lexclude(:)) + + if (size(lexclude) /= npl_new) then + allocate(lexclude_tmp(npl_new)) + call move_alloc(lexclude_tmp, lexclude) + end if + lexclude(1:npl_new) = .false. end associate @@ -438,13 +450,6 @@ subroutine calculate_system_energy(linclude_fragments) ! is in the main program. This will temporarily expand the planet list in a temporary system object called tmpsys to feed it into symba_energy associate(pl => system%pl, npl => system%pl%nbody, cb => system%cb) - where (lexclude(1:npl_new)) - tmpsys%pl%status(1:npl_new) = INACTIVE - elsewhere - tmpsys%pl%status(1:npl_new) = ACTIVE - end where - - ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by ! allocating two of these ginormous arrays simulteouously. This is not particularly efficient, but as this @@ -452,6 +457,9 @@ subroutine calculate_system_energy(linclude_fragments) lk_plpl = allocated(pl%k_plpl) if (lk_plpl) deallocate(pl%k_plpl) + call construct_temporary_system() + if (linclude_fragments) call add_fragments_to_tmpsys() + call tmpsys%pl%index(param) call tmpsys%get_energy_and_momentum(param)