From ac3491807d395d618774df4ba68f0012e33f8896 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 6 Aug 2021 15:27:47 -0400 Subject: [PATCH] Unified the variable names between integrators by refactoring mergesub_list to pl_discards (and mergeadd_list to pl_adds for consistency). Fixed some issues related to identifying collisions. --- docs/src/setup.f90 | 4 +- docs/src/symba_classes.f90 | 4 +- docs/src/symba_fragmentation.f90 | 72 +- docs/src/symba_io.f90 | 32 +- docs/src/symba_step.f90 | 6 +- docs/src/symba_util.f90 | 6 +- .../param.disruption_off_axis.in | 3 +- src/modules/symba_classes.f90 | 3 +- src/setup/setup.f90 | 4 +- src/symba/symba_fragmentation.f90 | 827 +++++++++--------- src/symba/symba_io.f90 | 90 +- src/symba/symba_step.f90 | 6 +- src/symba/symba_util.f90 | 17 +- 13 files changed, 541 insertions(+), 533 deletions(-) diff --git a/docs/src/setup.f90 b/docs/src/setup.f90 index 6cba6d27b..9346e8c12 100644 --- a/docs/src/setup.f90 +++ b/docs/src/setup.f90 @@ -54,8 +54,8 @@ module subroutine setup_construct_system(system, param) allocate(symba_pl :: system%pl) allocate(symba_tp :: system%tp) allocate(symba_tp :: system%tp_discards) - allocate(symba_merger :: system%mergeadd_list) - allocate(symba_merger :: system%mergesub_list) + allocate(symba_merger :: system%pl_adds) + allocate(symba_merger :: system%pl_discards) allocate(symba_plplenc :: system%plplenc_list) allocate(symba_pltpenc :: system%pltpenc_list) end select diff --git a/docs/src/symba_classes.f90 b/docs/src/symba_classes.f90 index 0e66ebf7c..4cf69a3e1 100644 --- a/docs/src/symba_classes.f90 +++ b/docs/src/symba_classes.f90 @@ -160,8 +160,8 @@ module symba_classes ! symba_nbody_system class definitions and method interfaces !******************************************************************************************************************************** type, extends(helio_nbody_system) :: symba_nbody_system - class(symba_merger), allocatable :: mergeadd_list !! List of added bodies in mergers or collisions - class(symba_merger), allocatable :: mergesub_list !! List of subtracted bodies in mergers or collisions + class(symba_merger), allocatable :: pl_adds !! List of added bodies in mergers or collisions + class(symba_merger), allocatable :: pl_discards !! List of subtracted bodies in mergers or collisions class(symba_pltpenc), allocatable :: pltpenc_list !! List of massive body-test particle encounters in a single step class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step integer(I4B) :: irec !! System recursion level diff --git a/docs/src/symba_fragmentation.f90 b/docs/src/symba_fragmentation.f90 index efdd8c0d7..9c13170af 100644 --- a/docs/src/symba_fragmentation.f90 +++ b/docs/src/symba_fragmentation.f90 @@ -34,7 +34,7 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, select type(pl => system%pl) class is (symba_pl) - associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb) + associate(pl_adds => system%pl_adds, pl_discards => system%pl_discards, cb => system%cb) ! Collisional fragments will be uniformly distributed around the pre-impact barycenter nfrag = NFRAG_DISRUPT allocate(m_frag(nfrag)) @@ -91,11 +91,11 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, lmask(:) = .false. lmask(family(:)) = .true. pl%status(family(:)) = MERGED - nstart = mergesub_list%nbody + 1 - nend = mergesub_list%nbody + nfamily - call mergesub_list%append(pl, lmask) + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + nfamily + call pl_discards%append(pl, lmask) ! Record how many bodies were subtracted in this event - mergesub_list%ncomp(nstart:nend) = nfamily + pl_discards%ncomp(nstart:nend) = nfamily allocate(plnew, mold=pl) call plnew%setup(nfrag, param) @@ -133,10 +133,10 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, end if ! Append the new merged body to the list and record how many we made - nstart = mergeadd_list%nbody + 1 - nend = mergeadd_list%nbody + plnew%nbody - call mergeadd_list%append(plnew) - mergeadd_list%ncomp(nstart:nend) = plnew%nbody + nstart = pl_adds%nbody + 1 + nend = pl_adds%nbody + plnew%nbody + call pl_adds%append(plnew) + pl_adds%ncomp(nstart:nend) = plnew%nbody call plnew%setup(0, param) deallocate(plnew) @@ -179,7 +179,7 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m select type(pl => system%pl) class is (symba_pl) - associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb) + associate(pl_adds => system%pl_adds, pl_discards => system%pl_discards, cb => system%cb) mtot = sum(mass(:)) xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot @@ -247,11 +247,11 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m lmask(:) = .false. lmask(family(:)) = .true. pl%status(family(:)) = MERGED - nstart = mergesub_list%nbody + 1 - nend = mergesub_list%nbody + nfamily - call mergesub_list%append(pl, lmask) + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + nfamily + call pl_discards%append(pl, lmask) ! Record how many bodies were subtracted in this event - mergesub_list%ncomp(nstart:nend) = nfamily + pl_discards%ncomp(nstart:nend) = nfamily allocate(plnew, mold=pl) call plnew%setup(nfrag, param) @@ -289,10 +289,10 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m end if ! Append the new merged body to the list and record how many we made - nstart = mergeadd_list%nbody + 1 - nend = mergeadd_list%nbody + plnew%nbody - call mergeadd_list%append(plnew) - mergeadd_list%ncomp(nstart:nend) = plnew%nbody + nstart = pl_adds%nbody + 1 + nend = pl_adds%nbody + plnew%nbody + call pl_adds%append(plnew) + pl_adds%ncomp(nstart:nend) = plnew%nbody call plnew%setup(0, param) deallocate(plnew) @@ -334,7 +334,7 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass, select type(pl => system%pl) class is (symba_pl) - associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb) + associate(pl_adds => system%pl_adds, pl_discards => system%pl_discards, cb => system%cb) status = MERGED write(*, '("Merging bodies ",99(I8,",",:))') pl%id(family(:)) mass_new = sum(mass(:)) @@ -386,11 +386,11 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass, lmask(:) = .false. lmask(family(:)) = .true. pl%status(family(:)) = MERGED - nstart = mergesub_list%nbody + 1 - nend = mergesub_list%nbody + nfamily - call mergesub_list%append(pl, lmask) + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + nfamily + call pl_discards%append(pl, lmask) ! Record how many bodies were subtracted in this event - mergesub_list%ncomp(nstart:nend) = nfamily + pl_discards%ncomp(nstart:nend) = nfamily ! Create the new merged body allocate(plnew, mold=pl) @@ -422,10 +422,10 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass, end if ! Append the new merged body to the list and record how many we made - nstart = mergeadd_list%nbody + 1 - nend = mergeadd_list%nbody + plnew%nbody - call mergeadd_list%append(plnew) - mergeadd_list%ncomp(nstart:nend) = plnew%nbody + nstart = pl_adds%nbody + 1 + nend = pl_adds%nbody + plnew%nbody + call pl_adds%append(plnew) + pl_adds%ncomp(nstart:nend) = plnew%nbody call plnew%setup(0, param) deallocate(plnew) @@ -468,7 +468,7 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family, select type(pl => system%pl) class is (symba_pl) - associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb) + associate(pl_adds => system%pl_adds, pl_discards => system%pl_discards, cb => system%cb) ! Collisional fragments will be uniformly distributed around the pre-impact barycenter nfrag = NFRAG_SUPERCAT allocate(m_frag(nfrag)) @@ -521,11 +521,11 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family, lmask(:) = .false. lmask(family(:)) = .true. pl%status(family(:)) = MERGED - nstart = mergesub_list%nbody + 1 - nend = mergesub_list%nbody + nfamily - call mergesub_list%append(pl, lmask) + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + nfamily + call pl_discards%append(pl, lmask) ! Record how many bodies were subtracted in this event - mergesub_list%ncomp(nstart:nend) = nfamily + pl_discards%ncomp(nstart:nend) = nfamily allocate(plnew, mold=pl) call plnew%setup(nfrag, param) @@ -563,10 +563,10 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family, end if ! Append the new merged body to the list and record how many we made - nstart = mergeadd_list%nbody + 1 - nend = mergeadd_list%nbody + plnew%nbody - call mergeadd_list%append(plnew) - mergeadd_list%ncomp(nstart:nend) = plnew%nbody + nstart = pl_adds%nbody + 1 + nend = pl_adds%nbody + plnew%nbody + call pl_adds%append(plnew) + pl_adds%ncomp(nstart:nend) = plnew%nbody call plnew%setup(0, param) deallocate(plnew) diff --git a/docs/src/symba_io.f90 b/docs/src/symba_io.f90 index 2e568dd7e..8e6420f61 100644 --- a/docs/src/symba_io.f90 +++ b/docs/src/symba_io.f90 @@ -225,10 +225,10 @@ module subroutine symba_io_write_discard(self, param) character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))' class(swiftest_body), allocatable :: pltemp - associate(pl => self%pl, npl => self%pl%nbody, mergesub_list => self%mergesub_list, mergeadd_list => self%mergeadd_list) + associate(pl => self%pl, npl => self%pl%nbody, pl_discards => self%pl_discards, pl_adds => self%pl_adds) if (self%tp_discards%nbody > 0) call io_write_discard(self, param) - if (mergesub_list%nbody == 0) return + if (pl_discards%nbody == 0) return select case(param%out_stat) case('APPEND') open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', iostat = ierr) @@ -240,31 +240,31 @@ module subroutine symba_io_write_discard(self, param) end select lfirst = .false. if (param%lgr) then - call mergesub_list%pv2v(param) - call mergeadd_list%pv2v(param) + call pl_discards%pv2v(param) + call pl_adds%pv2v(param) end if - write(LUN, HDRFMT) param%t, mergesub_list%nbody, param%lbig_discard + write(LUN, HDRFMT) param%t, pl_discards%nbody, param%lbig_discard iadd = 1 isub = 1 - do while (iadd <= mergeadd_list%nbody) - nadd = mergeadd_list%ncomp(iadd) - nsub = mergesub_list%ncomp(isub) + do while (iadd <= pl_adds%nbody) + nadd = pl_adds%ncomp(iadd) + nsub = pl_discards%ncomp(isub) do j = 1, nadd - if (iadd <= mergeadd_list%nbody) then - write(LUN, NAMEFMT) ADD, mergesub_list%id(iadd), mergesub_list%status(iadd) - write(LUN, VECFMT) mergeadd_list%xh(1, iadd), mergeadd_list%xh(2, iadd), mergeadd_list%xh(3, iadd) - write(LUN, VECFMT) mergeadd_list%vh(1, iadd), mergeadd_list%vh(2, iadd), mergeadd_list%vh(3, iadd) + if (iadd <= pl_adds%nbody) then + write(LUN, NAMEFMT) ADD, pl_discards%id(iadd), pl_discards%status(iadd) + write(LUN, VECFMT) pl_adds%xh(1, iadd), pl_adds%xh(2, iadd), pl_adds%xh(3, iadd) + write(LUN, VECFMT) pl_adds%vh(1, iadd), pl_adds%vh(2, iadd), pl_adds%vh(3, iadd) else exit end if iadd = iadd + 1 end do do j = 1, nsub - if (isub <= mergesub_list%nbody) then - write(LUN, NAMEFMT) SUB, mergesub_list%id(isub), mergesub_list%status(isub) - write(LUN, VECFMT) mergesub_list%xh(1, isub), mergesub_list%xh(2, isub), mergesub_list%xh(3, isub) - write(LUN, VECFMT) mergesub_list%vh(1, isub), mergesub_list%vh(2, isub), mergesub_list%vh(3, isub) + if (isub <= pl_discards%nbody) then + write(LUN, NAMEFMT) SUB, pl_discards%id(isub), pl_discards%status(isub) + write(LUN, VECFMT) pl_discards%xh(1, isub), pl_discards%xh(2, isub), pl_discards%xh(3, isub) + write(LUN, VECFMT) pl_discards%vh(1, isub), pl_discards%vh(2, isub), pl_discards%vh(3, isub) else exit end if diff --git a/docs/src/symba_step.f90 b/docs/src/symba_step.f90 index 41e7a3a74..7065625b4 100644 --- a/docs/src/symba_step.f90 +++ b/docs/src/symba_step.f90 @@ -237,7 +237,7 @@ module subroutine symba_step_reset_system(self) ! Internals integer(I4B) :: i - associate(system => self, pltpenc_list => self%pltpenc_list, plplenc_list => self%plplenc_list, mergeadd_list => self%mergeadd_list, mergesub_list => self%mergesub_list) + associate(system => self, pltpenc_list => self%pltpenc_list, plplenc_list => self%plplenc_list, pl_adds => self%pl_adds, pl_discards => self%pl_discards) select type(pl => system%pl) class is (symba_pl) select type(tp => system%tp) @@ -265,8 +265,8 @@ module subroutine symba_step_reset_system(self) pltpenc_list%nenc = 0 end if - call mergeadd_list%resize(0) - call mergesub_list%resize(0) + call pl_adds%resize(0) + call pl_discards%resize(0) end select end select end associate diff --git a/docs/src/symba_util.f90 b/docs/src/symba_util.f90 index 98c8889d8..0517cff9a 100644 --- a/docs/src/symba_util.f90 +++ b/docs/src/symba_util.f90 @@ -381,13 +381,13 @@ module subroutine symba_util_rearray_pl(self, system, param) ! Internals class(symba_pl), allocatable :: pl_discards !! The discarded body list. - associate(pl => self, mergeadd_list => system%mergeadd_list) + associate(pl => self, pl_adds => system%pl_adds) allocate(pl_discards, mold=pl) ! Remove the discards call pl%spill(pl_discards, lspill_list=(pl%ldiscard(:) .or. pl%status(:) == INACTIVE), ldestructive=.true.) ! Add in any new bodies - call pl%append(mergeadd_list) + call pl%append(pl_adds) ! If there are still bodies in the system, sort by mass in descending order and re-index if (pl%nbody > 0) then @@ -397,7 +397,7 @@ module subroutine symba_util_rearray_pl(self, system, param) call pl%eucl_index() end if - ! Destroy the discarded body list, since we already have what we need in the mergesub_list + ! Destroy the discarded body list, since we already have what we need in the pl_discards call pl_discards%setup(0,param) deallocate(pl_discards) end associate diff --git a/examples/symba_energy_momentum/param.disruption_off_axis.in b/examples/symba_energy_momentum/param.disruption_off_axis.in index b6f29564b..bb0915050 100644 --- a/examples/symba_energy_momentum/param.disruption_off_axis.in +++ b/examples/symba_energy_momentum/param.disruption_off_axis.in @@ -17,9 +17,8 @@ CHK_RMIN 0.005 CHK_RMAX 1e6 CHK_EJECT -1.0 ! ignore this check CHK_QMIN -1.0 ! ignore this check -!CHK_QMIN_COORD HELIO ! commented out here -!CHK_QMIN_RANGE 1.0 1000.0 ! commented out here ENC_OUT enc.disruption_off_axis.dat +DISCARD_OUT discard.disruption_off_axis.out EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT yes ! Hill's sphere radii in input file diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 0e66ebf7c..ed495ecfa 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -160,8 +160,7 @@ module symba_classes ! symba_nbody_system class definitions and method interfaces !******************************************************************************************************************************** type, extends(helio_nbody_system) :: symba_nbody_system - class(symba_merger), allocatable :: mergeadd_list !! List of added bodies in mergers or collisions - class(symba_merger), allocatable :: mergesub_list !! List of subtracted bodies in mergers or collisions + class(symba_merger), allocatable :: pl_adds !! List of added bodies in mergers or collisions class(symba_pltpenc), allocatable :: pltpenc_list !! List of massive body-test particle encounters in a single step class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step integer(I4B) :: irec !! System recursion level diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 6cba6d27b..9346e8c12 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -54,8 +54,8 @@ module subroutine setup_construct_system(system, param) allocate(symba_pl :: system%pl) allocate(symba_tp :: system%tp) allocate(symba_tp :: system%tp_discards) - allocate(symba_merger :: system%mergeadd_list) - allocate(symba_merger :: system%mergesub_list) + allocate(symba_merger :: system%pl_adds) + allocate(symba_merger :: system%pl_discards) allocate(symba_plplenc :: system%plplenc_list) allocate(symba_pltpenc :: system%pltpenc_list) end select diff --git a/src/symba/symba_fragmentation.f90 b/src/symba/symba_fragmentation.f90 index efdd8c0d7..a27d31e12 100644 --- a/src/symba/symba_fragmentation.f90 +++ b/src/symba/symba_fragmentation.f90 @@ -34,115 +34,117 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, select type(pl => system%pl) class is (symba_pl) - associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb) - ! Collisional fragments will be uniformly distributed around the pre-impact barycenter - nfrag = NFRAG_DISRUPT - allocate(m_frag(nfrag)) - allocate(rad_frag(nfrag)) - allocate(xb_frag(NDIM, nfrag)) - allocate(vb_frag(NDIM, nfrag)) - allocate(rot_frag(NDIM, nfrag)) - allocate(Ip_frag(NDIM, nfrag)) - - mtot = sum(mass(:)) - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot - - ! Get mass weighted mean of Ip and average density - Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot - vol(:) = 4._DP / 3._DP * PI * radius(:)**3 - avg_dens = mtot / sum(vol(:)) - - ! Distribute the mass among fragments, with a branch to check for the size of the second largest fragment - m_frag(1) = mass_res(1) - if (mass_res(2) > mass_res(1) / 3._DP) then - m_frag(2) = mass_res(2) - istart = 3 - else - istart = 2 - end if - ! Distribute remaining mass among the remaining bodies - do i = istart, nfrag - m_frag(i) = (mtot - sum(m_frag(1:istart - 1))) / (nfrag - istart + 1) - end do - - ! Distribute any residual mass if there is any and set the radius - m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) - rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) - - do i = 1, nfrag - Ip_frag(:, i) = Ip_new(:) - end do - - call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) - - if (lfailure) then - write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' - status = ACTIVE - nfrag = 0 - else - ! Populate the list of new bodies - write(*,'("Generating ",I2.0," fragments")') nfrag - status = DISRUPTION - - ! Add the family bodies to the subtraction list - nfamily = size(family(:)) - lmask(:) = .false. - lmask(family(:)) = .true. - pl%status(family(:)) = MERGED - nstart = mergesub_list%nbody + 1 - nend = mergesub_list%nbody + nfamily - call mergesub_list%append(pl, lmask) - ! Record how many bodies were subtracted in this event - mergesub_list%ncomp(nstart:nend) = nfamily - - allocate(plnew, mold=pl) - call plnew%setup(nfrag, param) - - plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] - system%maxid = system%maxid + nfrag - plnew%status(:) = ACTIVE - plnew%lcollision(:) = .false. - plnew%ldiscard(:) = .false. - plnew%xb(:,:) = xb_frag(:, :) - plnew%vb(:,:) = vb_frag(:, :) - do i = 1, nfrag - plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) - plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) + select type(pl_discards => system%pl_discards) + class is (symba_merger) + associate(pl_adds => system%pl_adds, cb => system%cb) + ! Collisional fragments will be uniformly distributed around the pre-impact barycenter + nfrag = NFRAG_DISRUPT + allocate(m_frag(nfrag)) + allocate(rad_frag(nfrag)) + allocate(xb_frag(NDIM, nfrag)) + allocate(vb_frag(NDIM, nfrag)) + allocate(rot_frag(NDIM, nfrag)) + allocate(Ip_frag(NDIM, nfrag)) + + mtot = sum(mass(:)) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot + + ! Get mass weighted mean of Ip and average density + Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot + vol(:) = 4._DP / 3._DP * PI * radius(:)**3 + avg_dens = mtot / sum(vol(:)) + + ! Distribute the mass among fragments, with a branch to check for the size of the second largest fragment + m_frag(1) = mass_res(1) + if (mass_res(2) > mass_res(1) / 3._DP) then + m_frag(2) = mass_res(2) + istart = 3 + else + istart = 2 + end if + ! Distribute remaining mass among the remaining bodies + do i = istart, nfrag + m_frag(i) = (mtot - sum(m_frag(1:istart - 1))) / (nfrag - istart + 1) end do - plnew%mass(:) = m_frag(:) - plnew%Gmass(:) = param%GU * m_frag(:) - plnew%density(:) = avg_dens - plnew%radius(:) = rad_frag(:) - plnew%info(:)%origin_type = "Disruption" - plnew%info(:)%origin_time = param%t + + ! Distribute any residual mass if there is any and set the radius + m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) + rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) + do i = 1, nfrag - plnew%info(i)%origin_xh(:) = plnew%xh(:,i) - plnew%info(i)%origin_vh(:) = plnew%vh(:,i) + Ip_frag(:, i) = Ip_new(:) end do - if (param%lrotation) then - plnew%Ip(:,:) = Ip_frag(:,:) - plnew%rot(:,:) = rot_frag(:,:) - end if - if (param%ltides) then - ibiggest = maxloc(pl%Gmass(family(:)), dim=1) - plnew%Q = pl%Q(ibiggest) - plnew%k2 = pl%k2(ibiggest) - plnew%tlag = pl%tlag(ibiggest) - end if - - ! Append the new merged body to the list and record how many we made - nstart = mergeadd_list%nbody + 1 - nend = mergeadd_list%nbody + plnew%nbody - call mergeadd_list%append(plnew) - mergeadd_list%ncomp(nstart:nend) = plnew%nbody - - call plnew%setup(0, param) - deallocate(plnew) - end if + + call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) + + if (lfailure) then + write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' + status = ACTIVE + nfrag = 0 + else + ! Populate the list of new bodies + write(*,'("Generating ",I2.0," fragments")') nfrag + status = DISRUPTION + + ! Add the family bodies to the subtraction list + nfamily = size(family(:)) + lmask(:) = .false. + lmask(family(:)) = .true. + pl%status(family(:)) = MERGED + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + nfamily + call pl_discards%append(pl, lmask) + ! Record how many bodies were subtracted in this event + pl_discards%ncomp(nstart:nend) = nfamily - end associate + allocate(plnew, mold=pl) + call plnew%setup(nfrag, param) + + plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] + system%maxid = system%maxid + nfrag + plnew%status(:) = ACTIVE + plnew%lcollision(:) = .false. + plnew%ldiscard(:) = .false. + plnew%xb(:,:) = xb_frag(:, :) + plnew%vb(:,:) = vb_frag(:, :) + do i = 1, nfrag + plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) + plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) + end do + plnew%mass(:) = m_frag(:) + plnew%Gmass(:) = param%GU * m_frag(:) + plnew%density(:) = avg_dens + plnew%radius(:) = rad_frag(:) + plnew%info(:)%origin_type = "Disruption" + plnew%info(:)%origin_time = param%t + do i = 1, nfrag + plnew%info(i)%origin_xh(:) = plnew%xh(:,i) + plnew%info(i)%origin_vh(:) = plnew%vh(:,i) + end do + if (param%lrotation) then + plnew%Ip(:,:) = Ip_frag(:,:) + plnew%rot(:,:) = rot_frag(:,:) + end if + if (param%ltides) then + ibiggest = maxloc(pl%Gmass(family(:)), dim=1) + plnew%Q = pl%Q(ibiggest) + plnew%k2 = pl%k2(ibiggest) + plnew%tlag = pl%tlag(ibiggest) + end if + + ! Append the new merged body to the list and record how many we made + nstart = pl_adds%nbody + 1 + nend = pl_adds%nbody + plnew%nbody + call pl_adds%append(plnew) + pl_adds%ncomp(nstart:nend) = plnew%nbody + + call plnew%setup(0, param) + deallocate(plnew) + end if + end associate + end select end select return @@ -179,126 +181,129 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m select type(pl => system%pl) class is (symba_pl) - associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb) - mtot = sum(mass(:)) - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot - lpure = .false. - - ! The largest body will stay untouched - if (mass(1) > mass(2)) then - jtarg = 1 - jproj = 2 - else - jtarg = 2 - jproj = 1 - end if - - if (mass_res(2) > 0.9_DP * mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched - write(*,*) 'Pure hit and run. No new fragments generated.' - nfrag = 0 - lpure = .true. - else ! Imperfect hit and run, so we'll keep the largest body and destroy the other - nfrag = NFRAG_DISRUPT - 1 + select type(pl_discards => system%pl_discards) + class is (symba_merger) + associate(pl_adds => system%pl_adds, cb => system%cb) + mtot = sum(mass(:)) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot lpure = .false. - allocate(m_frag(nfrag)) - allocate(id_frag(nfrag)) - allocate(rad_frag(nfrag)) - allocate(xb_frag(NDIM, nfrag)) - allocate(vb_frag(NDIM, nfrag)) - allocate(rot_frag(NDIM, nfrag)) - allocate(Ip_frag(NDIM, nfrag)) - m_frag(1) = mass(jtarg) - ibiggest = maxloc(pl%Gmass(family(:)), dim=1) - id_frag(1) = pl%id(ibiggest) - rad_frag(1) = radius(jtarg) - xb_frag(:, 1) = x(:, jtarg) - vb_frag(:, 1) = v(:, jtarg) - Ip_frag(:,1) = Ip(:, jtarg) - - ! Get mass weighted mean of Ip and average density - vol(:) = 4._DP / 3._DP * pi * radius(:)**3 - avg_dens = mass(jproj) / vol(jproj) - m_frag(2:nfrag) = (mtot - m_frag(1)) / (nfrag - 1) - rad_frag(2:nfrag) = (3 * m_frag(2:nfrag) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) - m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) - - do i = 1, nfrag - Ip_frag(:, i) = Ip(:, jproj) - end do - - ! Put the fragments on the circle surrounding the center of mass of the system - call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lpure) - if (lpure) then - write(*,*) 'Should have been a pure hit and run instead' - nfrag = 0 + + ! The largest body will stay untouched + if (mass(1) > mass(2)) then + jtarg = 1 + jproj = 2 else - write(*,'("Generating ",I2.0," fragments")') nfrag + jtarg = 2 + jproj = 1 end if - end if - if (lpure) then - status = ACTIVE - else - status = HIT_AND_RUN - - ! Add the family bodies to the subtraction list - nfamily = size(family(:)) - lmask(:) = .false. - lmask(family(:)) = .true. - pl%status(family(:)) = MERGED - nstart = mergesub_list%nbody + 1 - nend = mergesub_list%nbody + nfamily - call mergesub_list%append(pl, lmask) - ! Record how many bodies were subtracted in this event - mergesub_list%ncomp(nstart:nend) = nfamily - - allocate(plnew, mold=pl) - call plnew%setup(nfrag, param) - - plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] - system%maxid = system%maxid + nfrag - plnew%status(:) = ACTIVE - plnew%lcollision(:) = .false. - plnew%ldiscard(:) = .false. - plnew%xb(:,:) = xb_frag(:, :) - plnew%vb(:,:) = vb_frag(:, :) - do i = 1, nfrag - plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) - plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) - end do - plnew%mass(:) = m_frag(:) - plnew%Gmass(:) = param%GU * m_frag(:) - plnew%density(:) = avg_dens - plnew%radius(:) = rad_frag(:) - plnew%info(:)%origin_type = "Hit and run fragment" - plnew%info(:)%origin_time = param%t - do i = 1, nfrag - plnew%info(i)%origin_xh(:) = plnew%xh(:,i) - plnew%info(i)%origin_vh(:) = plnew%vh(:,i) - end do - if (param%lrotation) then - plnew%Ip(:,:) = Ip_frag(:,:) - plnew%rot(:,:) = rot_frag(:,:) - end if - if (param%ltides) then + + if (mass_res(2) > 0.9_DP * mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched + write(*,*) 'Pure hit and run. No new fragments generated.' + nfrag = 0 + lpure = .true. + else ! Imperfect hit and run, so we'll keep the largest body and destroy the other + nfrag = NFRAG_DISRUPT - 1 + lpure = .false. + allocate(m_frag(nfrag)) + allocate(id_frag(nfrag)) + allocate(rad_frag(nfrag)) + allocate(xb_frag(NDIM, nfrag)) + allocate(vb_frag(NDIM, nfrag)) + allocate(rot_frag(NDIM, nfrag)) + allocate(Ip_frag(NDIM, nfrag)) + m_frag(1) = mass(jtarg) ibiggest = maxloc(pl%Gmass(family(:)), dim=1) - plnew%Q = pl%Q(ibiggest) - plnew%k2 = pl%k2(ibiggest) - plnew%tlag = pl%tlag(ibiggest) + id_frag(1) = pl%id(ibiggest) + rad_frag(1) = radius(jtarg) + xb_frag(:, 1) = x(:, jtarg) + vb_frag(:, 1) = v(:, jtarg) + Ip_frag(:,1) = Ip(:, jtarg) + + ! Get mass weighted mean of Ip and average density + vol(:) = 4._DP / 3._DP * pi * radius(:)**3 + avg_dens = mass(jproj) / vol(jproj) + m_frag(2:nfrag) = (mtot - m_frag(1)) / (nfrag - 1) + rad_frag(2:nfrag) = (3 * m_frag(2:nfrag) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) + m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) + + do i = 1, nfrag + Ip_frag(:, i) = Ip(:, jproj) + end do + + ! Put the fragments on the circle surrounding the center of mass of the system + call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lpure) + if (lpure) then + write(*,*) 'Should have been a pure hit and run instead' + nfrag = 0 + else + write(*,'("Generating ",I2.0," fragments")') nfrag + end if end if - - ! Append the new merged body to the list and record how many we made - nstart = mergeadd_list%nbody + 1 - nend = mergeadd_list%nbody + plnew%nbody - call mergeadd_list%append(plnew) - mergeadd_list%ncomp(nstart:nend) = plnew%nbody - - call plnew%setup(0, param) - deallocate(plnew) + if (lpure) then + status = ACTIVE + else + status = HIT_AND_RUN + + ! Add the family bodies to the subtraction list + nfamily = size(family(:)) + lmask(:) = .false. + lmask(family(:)) = .true. + pl%status(family(:)) = MERGED + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + nfamily + call pl_discards%append(pl, lmask) + ! Record how many bodies were subtracted in this event + pl_discards%ncomp(nstart:nend) = nfamily + + allocate(plnew, mold=pl) + call plnew%setup(nfrag, param) + + plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] + system%maxid = system%maxid + nfrag + plnew%status(:) = ACTIVE + plnew%lcollision(:) = .false. + plnew%ldiscard(:) = .false. + plnew%xb(:,:) = xb_frag(:, :) + plnew%vb(:,:) = vb_frag(:, :) + do i = 1, nfrag + plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) + plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) + end do + plnew%mass(:) = m_frag(:) + plnew%Gmass(:) = param%GU * m_frag(:) + plnew%density(:) = avg_dens + plnew%radius(:) = rad_frag(:) + plnew%info(:)%origin_type = "Hit and run fragment" + plnew%info(:)%origin_time = param%t + do i = 1, nfrag + plnew%info(i)%origin_xh(:) = plnew%xh(:,i) + plnew%info(i)%origin_vh(:) = plnew%vh(:,i) + end do + if (param%lrotation) then + plnew%Ip(:,:) = Ip_frag(:,:) + plnew%rot(:,:) = rot_frag(:,:) + end if + if (param%ltides) then + ibiggest = maxloc(pl%Gmass(family(:)), dim=1) + plnew%Q = pl%Q(ibiggest) + plnew%k2 = pl%k2(ibiggest) + plnew%tlag = pl%tlag(ibiggest) + end if + + ! Append the new merged body to the list and record how many we made + nstart = pl_adds%nbody + 1 + nend = pl_adds%nbody + plnew%nbody + call pl_adds%append(plnew) + pl_adds%ncomp(nstart:nend) = plnew%nbody + + call plnew%setup(0, param) + deallocate(plnew) - end if - end associate + end if + end associate + end select end select return @@ -334,103 +339,105 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass, select type(pl => system%pl) class is (symba_pl) - associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb) - status = MERGED - write(*, '("Merging bodies ",99(I8,",",:))') pl%id(family(:)) - mass_new = sum(mass(:)) - - ! Merged body is created at the barycenter of the original bodies - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mass_new - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mass_new - - ! Get mass weighted mean of Ip and - vol(:) = 4._DP / 3._DP * PI * radius(:)**3 - volume_new = sum(vol(:)) - radius_new = (3 * volume_new / (4 * PI))**(1._DP / 3._DP) - - L_orb_old(:) = 0.0_DP - - ! Compute orbital angular momentum of pre-impact system - do i = 1, 2 - xc(:) = x(:, i) - xcom(:) - vc(:) = v(:, i) - vcom(:) - xcrossv(:) = xc(:) .cross. vc(:) - L_orb_old(:) = L_orb_old(:) + mass(i) * xcrossv(:) - end do + select type(pl_discards => system%pl_discards) + class is (symba_merger) + associate(pl_adds => system%pl_adds, cb => system%cb) + status = MERGED + write(*, '("Merging bodies ",99(I8,",",:))') pl%id(family(:)) + mass_new = sum(mass(:)) - if (param%lrotation) then - Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mass_new - L_spin_old(:) = L_spin(:,1) + L_spin(:,2) + ! Merged body is created at the barycenter of the original bodies + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mass_new + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mass_new + + ! Get mass weighted mean of Ip and + vol(:) = 4._DP / 3._DP * PI * radius(:)**3 + volume_new = sum(vol(:)) + radius_new = (3 * volume_new / (4 * PI))**(1._DP / 3._DP) - ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body - L_spin_new(:) = L_orb_old(:) + L_spin_old(:) - - ! Assume prinicpal axis rotation on 3rd Ip axis - rot_new(:) = L_spin_new(:) / (Ip_new(3) * mass_new * radius_new**2) - else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable - system%Lescape(:) = system%Lescape(:) + L_orb_old(:) - end if - - ! Keep track of the component of potential energy due to the pre-impact family for book-keeping - nfamily = size(family(:)) - pe = 0.0_DP - do j = 1, nfamily - do i = j + 1, nfamily - pe = pe - pl%mass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) + L_orb_old(:) = 0.0_DP + + ! Compute orbital angular momentum of pre-impact system + do i = 1, 2 + xc(:) = x(:, i) - xcom(:) + vc(:) = v(:, i) - vcom(:) + xcrossv(:) = xc(:) .cross. vc(:) + L_orb_old(:) = L_orb_old(:) + mass(i) * xcrossv(:) end do - end do - system%Ecollisions = system%Ecollisions + pe - system%Euntracked = system%Euntracked - pe - - ! Add the family bodies to the subtraction list - lmask(:) = .false. - lmask(family(:)) = .true. - pl%status(family(:)) = MERGED - nstart = mergesub_list%nbody + 1 - nend = mergesub_list%nbody + nfamily - call mergesub_list%append(pl, lmask) - ! Record how many bodies were subtracted in this event - mergesub_list%ncomp(nstart:nend) = nfamily + + if (param%lrotation) then + Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mass_new + L_spin_old(:) = L_spin(:,1) + L_spin(:,2) - ! Create the new merged body - allocate(plnew, mold=pl) - call plnew%setup(1, param) + ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body + L_spin_new(:) = L_orb_old(:) + L_spin_old(:) + + ! Assume prinicpal axis rotation on 3rd Ip axis + rot_new(:) = L_spin_new(:) / (Ip_new(3) * mass_new * radius_new**2) + else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable + system%Lescape(:) = system%Lescape(:) + L_orb_old(:) + end if + + ! Keep track of the component of potential energy due to the pre-impact family for book-keeping + nfamily = size(family(:)) + pe = 0.0_DP + do j = 1, nfamily + do i = j + 1, nfamily + pe = pe - pl%mass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) + end do + end do + system%Ecollisions = system%Ecollisions + pe + system%Euntracked = system%Euntracked - pe + + ! Add the family bodies to the subtraction list + lmask(:) = .false. + lmask(family(:)) = .true. + pl%status(family(:)) = MERGED + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + nfamily + call pl_discards%append(pl, lmask) + ! Record how many bodies were subtracted in this event + pl_discards%ncomp(nstart:nend) = nfamily - ! The merged body's name will be that of the largest of the two parents - ibiggest = maxloc(pl%Gmass(family(:)), dim=1) - plnew%id(1) = pl%id(family(ibiggest)) - plnew%status(1) = ACTIVE - plnew%lcollision = .false. - plnew%ldiscard = .false. - plnew%xb(:,1) = xcom(:) - plnew%vb(:,1) = vcom(:) - plnew%xh(:,1) = xcom(:) - cb%xb(:) - plnew%vh(:,1) = vcom(:) - cb%vb(:) - plnew%mass(1) = mass_new - plnew%Gmass(1) = param%GU * mass_new - plnew%density(1) = mass_new / volume_new - plnew%radius(1) = radius_new - plnew%info(1) = pl%info(family(ibiggest)) - if (param%lrotation) then - plnew%Ip(:,1) = Ip_new(:) - plnew%rot(:,1) = rot_new(:) - end if - if (param%ltides) then - plnew%Q = pl%Q(ibiggest) - plnew%k2 = pl%k2(ibiggest) - plnew%tlag = pl%tlag(ibiggest) - end if + ! Create the new merged body + allocate(plnew, mold=pl) + call plnew%setup(1, param) - ! Append the new merged body to the list and record how many we made - nstart = mergeadd_list%nbody + 1 - nend = mergeadd_list%nbody + plnew%nbody - call mergeadd_list%append(plnew) - mergeadd_list%ncomp(nstart:nend) = plnew%nbody + ! The merged body's name will be that of the largest of the two parents + ibiggest = maxloc(pl%Gmass(family(:)), dim=1) + plnew%id(1) = pl%id(family(ibiggest)) + plnew%status(1) = ACTIVE + plnew%lcollision = .false. + plnew%ldiscard = .false. + plnew%xb(:,1) = xcom(:) + plnew%vb(:,1) = vcom(:) + plnew%xh(:,1) = xcom(:) - cb%xb(:) + plnew%vh(:,1) = vcom(:) - cb%vb(:) + plnew%mass(1) = mass_new + plnew%Gmass(1) = param%GU * mass_new + plnew%density(1) = mass_new / volume_new + plnew%radius(1) = radius_new + plnew%info(1) = pl%info(family(ibiggest)) + if (param%lrotation) then + plnew%Ip(:,1) = Ip_new(:) + plnew%rot(:,1) = rot_new(:) + end if + if (param%ltides) then + plnew%Q = pl%Q(ibiggest) + plnew%k2 = pl%k2(ibiggest) + plnew%tlag = pl%tlag(ibiggest) + end if - call plnew%setup(0, param) - deallocate(plnew) + ! Append the new merged body to the list and record how many we made + nstart = pl_adds%nbody + 1 + nend = pl_adds%nbody + plnew%nbody + call pl_adds%append(plnew) + pl_adds%ncomp(nstart:nend) = plnew%nbody - end associate + call plnew%setup(0, param) + deallocate(plnew) + end associate + end select end select return @@ -468,111 +475,113 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family, select type(pl => system%pl) class is (symba_pl) - associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb) - ! Collisional fragments will be uniformly distributed around the pre-impact barycenter - nfrag = NFRAG_SUPERCAT - allocate(m_frag(nfrag)) - allocate(rad_frag(nfrag)) - allocate(xb_frag(NDIM, nfrag)) - allocate(vb_frag(NDIM, nfrag)) - allocate(rot_frag(NDIM, nfrag)) - allocate(Ip_frag(NDIM, nfrag)) - - mtot = sum(mass(:)) - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot - - ! Get mass weighted mean of Ip and average density - Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot - vol(:) = 4._DP / 3._DP * pi * radius(:)**3 - avg_dens = mtot / sum(vol(:)) - - ! If we are adding the first and largest fragment (lr), check to see if its mass is SMALLER than an equal distribution of - ! mass between all fragments. If so, we will just distribute the mass equally between the fragments - min_frag_mass = mtot / nfrag - if (mass_res(1) < min_frag_mass) then - m_frag(:) = min_frag_mass - else - m_frag(1) = mass_res(1) - m_frag(2:nfrag) = (mtot - mass_res(1)) / (nfrag - 1) - end if - ! Distribute any residual mass if there is any and set the radius - m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) - rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) - - do i = 1, nfrag - Ip_frag(:, i) = Ip_new(:) - end do + select type(pl_discards => system%pl_discards) + class is (symba_merger) + associate(pl_adds => system%pl_adds, cb => system%cb) + ! Collisional fragments will be uniformly distributed around the pre-impact barycenter + nfrag = NFRAG_SUPERCAT + allocate(m_frag(nfrag)) + allocate(rad_frag(nfrag)) + allocate(xb_frag(NDIM, nfrag)) + allocate(vb_frag(NDIM, nfrag)) + allocate(rot_frag(NDIM, nfrag)) + allocate(Ip_frag(NDIM, nfrag)) + + mtot = sum(mass(:)) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot + + ! Get mass weighted mean of Ip and average density + Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot + vol(:) = 4._DP / 3._DP * pi * radius(:)**3 + avg_dens = mtot / sum(vol(:)) + + ! If we are adding the first and largest fragment (lr), check to see if its mass is SMALLER than an equal distribution of + ! mass between all fragments. If so, we will just distribute the mass equally between the fragments + min_frag_mass = mtot / nfrag + if (mass_res(1) < min_frag_mass) then + m_frag(:) = min_frag_mass + else + m_frag(1) = mass_res(1) + m_frag(2:nfrag) = (mtot - mass_res(1)) / (nfrag - 1) + end if + ! Distribute any residual mass if there is any and set the radius + m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) + rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) + + do i = 1, nfrag + Ip_frag(:, i) = Ip_new(:) + end do - call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) - - if (lfailure) then - write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' - status = ACTIVE - nfrag = 0 - else - ! Populate the list of new bodies - write(*,'("Generating ",I2.0," fragments")') nfrag - status = SUPERCATASTROPHIC + call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) + + if (lfailure) then + write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' + status = ACTIVE + nfrag = 0 + else + ! Populate the list of new bodies + write(*,'("Generating ",I2.0," fragments")') nfrag + status = SUPERCATASTROPHIC - ! Add the family bodies to the subtraction list - nfamily = size(family(:)) - lmask(:) = .false. - lmask(family(:)) = .true. - pl%status(family(:)) = MERGED - nstart = mergesub_list%nbody + 1 - nend = mergesub_list%nbody + nfamily - call mergesub_list%append(pl, lmask) - ! Record how many bodies were subtracted in this event - mergesub_list%ncomp(nstart:nend) = nfamily + ! Add the family bodies to the subtraction list + nfamily = size(family(:)) + lmask(:) = .false. + lmask(family(:)) = .true. + pl%status(family(:)) = MERGED + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + nfamily + call pl_discards%append(pl, lmask) + ! Record how many bodies were subtracted in this event + pl_discards%ncomp(nstart:nend) = nfamily - allocate(plnew, mold=pl) - call plnew%setup(nfrag, param) - - plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] - system%maxid = system%maxid + nfrag - plnew%status(:) = ACTIVE - plnew%lcollision(:) = .false. - plnew%ldiscard(:) = .false. - plnew%xb(:,:) = xb_frag(:, :) - plnew%vb(:,:) = vb_frag(:, :) - do i = 1, nfrag - plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) - plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) - end do - plnew%mass(:) = m_frag(:) - plnew%Gmass(:) = param%GU * m_frag(:) - plnew%density(:) = avg_dens - plnew%radius(:) = rad_frag(:) - plnew%info(:)%origin_type = "Supercatastrophic" - plnew%info(:)%origin_time = param%t - do i = 1, nfrag - plnew%info(i)%origin_xh(:) = plnew%xh(:,i) - plnew%info(i)%origin_vh(:) = plnew%vh(:,i) - end do - if (param%lrotation) then - plnew%Ip(:,:) = Ip_frag(:,:) - plnew%rot(:,:) = rot_frag(:,:) - end if - if (param%ltides) then - ibiggest = maxloc(pl%Gmass(family(:)), dim=1) - plnew%Q = pl%Q(ibiggest) - plnew%k2 = pl%k2(ibiggest) - plnew%tlag = pl%tlag(ibiggest) + allocate(plnew, mold=pl) + call plnew%setup(nfrag, param) + + plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] + system%maxid = system%maxid + nfrag + plnew%status(:) = ACTIVE + plnew%lcollision(:) = .false. + plnew%ldiscard(:) = .false. + plnew%xb(:,:) = xb_frag(:, :) + plnew%vb(:,:) = vb_frag(:, :) + do i = 1, nfrag + plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) + plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) + end do + plnew%mass(:) = m_frag(:) + plnew%Gmass(:) = param%GU * m_frag(:) + plnew%density(:) = avg_dens + plnew%radius(:) = rad_frag(:) + plnew%info(:)%origin_type = "Supercatastrophic" + plnew%info(:)%origin_time = param%t + do i = 1, nfrag + plnew%info(i)%origin_xh(:) = plnew%xh(:,i) + plnew%info(i)%origin_vh(:) = plnew%vh(:,i) + end do + if (param%lrotation) then + plnew%Ip(:,:) = Ip_frag(:,:) + plnew%rot(:,:) = rot_frag(:,:) + end if + if (param%ltides) then + ibiggest = maxloc(pl%Gmass(family(:)), dim=1) + plnew%Q = pl%Q(ibiggest) + plnew%k2 = pl%k2(ibiggest) + plnew%tlag = pl%tlag(ibiggest) + end if + + ! Append the new merged body to the list and record how many we made + nstart = pl_adds%nbody + 1 + nend = pl_adds%nbody + plnew%nbody + call pl_adds%append(plnew) + pl_adds%ncomp(nstart:nend) = plnew%nbody + + call plnew%setup(0, param) + deallocate(plnew) end if - - ! Append the new merged body to the list and record how many we made - nstart = mergeadd_list%nbody + 1 - nend = mergeadd_list%nbody + plnew%nbody - call mergeadd_list%append(plnew) - mergeadd_list%ncomp(nstart:nend) = plnew%nbody - - call plnew%setup(0, param) - deallocate(plnew) - end if - - end associate + end associate + end select end select return diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 2e568dd7e..4401fb5db 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -225,54 +225,56 @@ module subroutine symba_io_write_discard(self, param) character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))' class(swiftest_body), allocatable :: pltemp - associate(pl => self%pl, npl => self%pl%nbody, mergesub_list => self%mergesub_list, mergeadd_list => self%mergeadd_list) + associate(pl => self%pl, npl => self%pl%nbody, pl_adds => self%pl_adds) if (self%tp_discards%nbody > 0) call io_write_discard(self, param) + select type(pl_discards => self%pl_discards) + class is (symba_merger) + if (pl_discards%nbody == 0) return + select case(param%out_stat) + case('APPEND') + open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', iostat = ierr) + case('NEW', 'REPLACE', 'UNKNOWN') + open(unit = LUN, file = param%discard_out, status = param%out_stat, form = 'FORMATTED', iostat = ierr) + case default + write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) + call util_exit(FAILURE) + end select + lfirst = .false. + if (param%lgr) then + call pl_discards%pv2v(param) + call pl_adds%pv2v(param) + end if - if (mergesub_list%nbody == 0) return - select case(param%out_stat) - case('APPEND') - open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', iostat = ierr) - case('NEW', 'REPLACE', 'UNKNOWN') - open(unit = LUN, file = param%discard_out, status = param%out_stat, form = 'FORMATTED', iostat = ierr) - case default - write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) - call util_exit(FAILURE) - end select - lfirst = .false. - if (param%lgr) then - call mergesub_list%pv2v(param) - call mergeadd_list%pv2v(param) - end if - - write(LUN, HDRFMT) param%t, mergesub_list%nbody, param%lbig_discard - iadd = 1 - isub = 1 - do while (iadd <= mergeadd_list%nbody) - nadd = mergeadd_list%ncomp(iadd) - nsub = mergesub_list%ncomp(isub) - do j = 1, nadd - if (iadd <= mergeadd_list%nbody) then - write(LUN, NAMEFMT) ADD, mergesub_list%id(iadd), mergesub_list%status(iadd) - write(LUN, VECFMT) mergeadd_list%xh(1, iadd), mergeadd_list%xh(2, iadd), mergeadd_list%xh(3, iadd) - write(LUN, VECFMT) mergeadd_list%vh(1, iadd), mergeadd_list%vh(2, iadd), mergeadd_list%vh(3, iadd) - else - exit - end if - iadd = iadd + 1 - end do - do j = 1, nsub - if (isub <= mergesub_list%nbody) then - write(LUN, NAMEFMT) SUB, mergesub_list%id(isub), mergesub_list%status(isub) - write(LUN, VECFMT) mergesub_list%xh(1, isub), mergesub_list%xh(2, isub), mergesub_list%xh(3, isub) - write(LUN, VECFMT) mergesub_list%vh(1, isub), mergesub_list%vh(2, isub), mergesub_list%vh(3, isub) - else - exit - end if - isub = isub + 1 + write(LUN, HDRFMT) param%t, pl_discards%nbody, param%lbig_discard + iadd = 1 + isub = 1 + do while (iadd <= pl_adds%nbody) + nadd = pl_adds%ncomp(iadd) + nsub = pl_discards%ncomp(isub) + do j = 1, nadd + if (iadd <= pl_adds%nbody) then + write(LUN, NAMEFMT) ADD, pl_discards%id(iadd), pl_discards%status(iadd) + write(LUN, VECFMT) pl_adds%xh(1, iadd), pl_adds%xh(2, iadd), pl_adds%xh(3, iadd) + write(LUN, VECFMT) pl_adds%vh(1, iadd), pl_adds%vh(2, iadd), pl_adds%vh(3, iadd) + else + exit + end if + iadd = iadd + 1 + end do + do j = 1, nsub + if (isub <= pl_discards%nbody) then + write(LUN, NAMEFMT) SUB, pl_discards%id(isub), pl_discards%status(isub) + write(LUN, VECFMT) pl_discards%xh(1, isub), pl_discards%xh(2, isub), pl_discards%xh(3, isub) + write(LUN, VECFMT) pl_discards%vh(1, isub), pl_discards%vh(2, isub), pl_discards%vh(3, isub) + else + exit + end if + isub = isub + 1 + end do end do - end do - close(LUN) + close(LUN) + end select end associate return diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 41e7a3a74..7065625b4 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -237,7 +237,7 @@ module subroutine symba_step_reset_system(self) ! Internals integer(I4B) :: i - associate(system => self, pltpenc_list => self%pltpenc_list, plplenc_list => self%plplenc_list, mergeadd_list => self%mergeadd_list, mergesub_list => self%mergesub_list) + associate(system => self, pltpenc_list => self%pltpenc_list, plplenc_list => self%plplenc_list, pl_adds => self%pl_adds, pl_discards => self%pl_discards) select type(pl => system%pl) class is (symba_pl) select type(tp => system%tp) @@ -265,8 +265,8 @@ module subroutine symba_step_reset_system(self) pltpenc_list%nenc = 0 end if - call mergeadd_list%resize(0) - call mergesub_list%resize(0) + call pl_adds%resize(0) + call pl_discards%resize(0) end select end select end associate diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 98c8889d8..2ebc11ebb 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -379,15 +379,17 @@ module subroutine symba_util_rearray_pl(self, system, param) class(symba_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(symba_parameters), intent(in) :: param !! Current run configuration parameters ! Internals - class(symba_pl), allocatable :: pl_discards !! The discarded body list. + class(symba_pl), allocatable :: tmp !! The discarded body list. - associate(pl => self, mergeadd_list => system%mergeadd_list) - allocate(pl_discards, mold=pl) - ! Remove the discards - call pl%spill(pl_discards, lspill_list=(pl%ldiscard(:) .or. pl%status(:) == INACTIVE), ldestructive=.true.) + associate(pl => self, pl_adds => system%pl_adds) + allocate(tmp, mold=pl) + ! Remove the discards and destroy the list, as the system already tracks pl_discards elsewhere + call pl%spill(tmp, lspill_list=(pl%ldiscard(:) .or. pl%status(:) == INACTIVE), ldestructive=.true.) + call tmp%setup(0,param) + deallocate(tmp) ! Add in any new bodies - call pl%append(mergeadd_list) + call pl%append(pl_adds) ! If there are still bodies in the system, sort by mass in descending order and re-index if (pl%nbody > 0) then @@ -397,9 +399,6 @@ module subroutine symba_util_rearray_pl(self, system, param) call pl%eucl_index() end if - ! Destroy the discarded body list, since we already have what we need in the mergesub_list - call pl_discards%setup(0,param) - deallocate(pl_discards) end associate return