From ac3491807d395d618774df4ba68f0012e33f8896 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 6 Aug 2021 15:27:47 -0400 Subject: [PATCH 1/7] 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 From 90faadf0a71236af13feaa135c1a99e5026ce19d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 6 Aug 2021 15:58:07 -0400 Subject: [PATCH 2/7] Refactored MTINY into GMTINY to be consistent in distinguishing G*mass and mass terms --- docs/src/fragmentation.f90 | 10 +++++----- docs/src/io.f90 | 2 +- docs/src/swiftest_classes.f90 | 4 ++-- docs/src/symba_classes.f90 | 8 ++++---- docs/src/symba_collision.f90 | 2 +- docs/src/symba_io.f90 | 12 ++++++------ docs/src/symba_setup.f90 | 2 +- docs/src/symba_util.f90 | 2 +- .../param.disruption_headon.in | 4 +--- .../param.disruption_off_axis.in | 2 +- examples/symba_energy_momentum/param.escape.in | 4 +--- examples/symba_energy_momentum/param.sun.in | 4 +--- .../param.supercatastrophic_headon.in | 4 +--- .../param.supercatastrophic_off_axis.in | 2 +- .../1pl_1pl_encounter/init_cond.py | 2 +- .../1pl_1pl_encounter/param.swiftest.in | 2 +- .../1pl_1tp_encounter/init_cond.py | 2 +- .../1pl_1tp_encounter/param.swiftest.in | 2 +- .../8pl_16tp_encounters/init_cond.py | 2 +- .../8pl_16tp_encounters/param.swiftest.in | 2 +- python/swiftest/swiftest/io.py | 18 +++++++++--------- src/fragmentation/fragmentation.f90 | 14 +++++++------- src/io/io.f90 | 2 +- src/modules/swiftest_classes.f90 | 4 ++-- src/modules/symba_classes.f90 | 8 ++++---- src/symba/symba_collision.f90 | 14 +++++++------- src/symba/symba_io.f90 | 12 ++++++------ src/symba/symba_setup.f90 | 2 +- src/symba/symba_util.f90 | 2 +- src/util/util_get_energy_momentum.f90 | 2 +- 30 files changed, 72 insertions(+), 80 deletions(-) diff --git a/docs/src/fragmentation.f90 b/docs/src/fragmentation.f90 index 460060183..c90f64cb4 100644 --- a/docs/src/fragmentation.f90 +++ b/docs/src/fragmentation.f90 @@ -369,7 +369,7 @@ subroutine calculate_system_energy(linclude_fragments) class is (symba_pl) select type(param) class is (symba_parameters) - plwksp%nplm = count(plwksp%Gmass > param%mtiny / mscale) + plwksp%nplm = count(plwksp%Gmass > param%Gmtiny / mscale) end select end select call tmpsys%pl%eucl_index() @@ -381,7 +381,7 @@ subroutine calculate_system_energy(linclude_fragments) class is (symba_pl) select type(param) class is (symba_parameters) - nplm = count(pl%mass > param%mtiny) + nplm = count(pl%mass > param%Gmtiny) end select end select if (lk_plpl) call pl%eucl_index() @@ -836,7 +836,7 @@ end subroutine fragmentation_initialize - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, Gmtiny, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine the collisional regime of two colliding bodies. @@ -857,7 +857,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v ! Arguments integer(I4B), intent(out) :: regime real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, Gmtiny real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 real(DP), intent(out) :: Qloss !! The residual energy after the collision ! Constants @@ -931,7 +931,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v Qloss = 0.0_DP U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27 - if ((m1 < mtiny).or.(m2 < mtiny)) then + if ((m1 < Gmtiny).or.(m2 < Gmtiny)) then regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime Mlr = Mtot Mslr = 0.0_DP diff --git a/docs/src/io.f90 b/docs/src/io.f90 index de1226951..1e0e8d626 100644 --- a/docs/src/io.f90 +++ b/docs/src/io.f90 @@ -493,7 +493,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) read(param_value, *) param%Ecollisions case("EUNTRACKED") read(param_value, *) param%Euntracked - case ("NPLMAX", "NTPMAX", "MTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters + case ("NPLMAX", "NTPMAX", "GMTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters case default write(iomsg,*) "Unknown parameter -> ",param_name iostat = -1 diff --git a/docs/src/swiftest_classes.f90 b/docs/src/swiftest_classes.f90 index 2455e77f2..051add1aa 100644 --- a/docs/src/swiftest_classes.f90 +++ b/docs/src/swiftest_classes.f90 @@ -469,11 +469,11 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, real(DP), intent(inout) :: Qloss end subroutine fragmentation_initialize - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, Gmtiny, Qloss) implicit none integer(I4B), intent(out) :: regime real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, Gmtiny real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 real(DP), intent(out) :: Qloss !! The residual energy after the collision end subroutine fragmentation_regime diff --git a/docs/src/symba_classes.f90 b/docs/src/symba_classes.f90 index 4cf69a3e1..6ecb27751 100644 --- a/docs/src/symba_classes.f90 +++ b/docs/src/symba_classes.f90 @@ -19,7 +19,7 @@ module symba_classes type, extends(swiftest_parameters) :: symba_parameters character(STRMAX) :: particle_file = PARTICLE_OUTFILE !! Name of output particle information file - real(DP) :: MTINY = -1.0_DP !! Smallest mass that is fully gravitating + real(DP) :: GMTINY = -1.0_DP !! Smallest mass that is fully gravitating integer(I4B), dimension(:), allocatable :: seed !! Random seeds logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. contains @@ -73,9 +73,9 @@ module symba_classes type, extends(helio_pl) :: symba_pl logical, dimension(:), allocatable :: lcollision !! flag indicating whether body has merged with another this time step logical, dimension(:), allocatable :: lencounter !! flag indicating whether body is part of an encounter this time step - logical, dimension(:), allocatable :: lmtiny !! flag indicating whether this body is below the MTINY cutoff value - integer(I4B) :: nplm !! number of bodies above the MTINY limit - integer(I8B) :: nplplm !! Number of body (all massive)-body (only those above MTINY) comparisons in the flattened upper triangular matrix + logical, dimension(:), allocatable :: lmtiny !! flag indicating whether this body is below the GMTINY cutoff value + integer(I4B) :: nplm !! number of bodies above the GMTINY limit + integer(I8B) :: nplplm !! Number of body (all massive)-body (only those above GMTINY) comparisons in the flattened upper triangular matrix integer(I4B), dimension(:), allocatable :: nplenc !! number of encounters with other planets this time step integer(I4B), dimension(:), allocatable :: ntpenc !! number of encounters with test particles this time step integer(I4B), dimension(:), allocatable :: levelg !! level at which this body should be moved diff --git a/docs/src/symba_collision.f90 b/docs/src/symba_collision.f90 index 952d59709..caffdf296 100644 --- a/docs/src/symba_collision.f90 +++ b/docs/src/symba_collision.f90 @@ -449,7 +449,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) v2_si(:) = plpl_collisions%v2(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The collective density of the parent and its children Mcb_si = cb%mass * param%MU2KG - mtiny_si = (param%MTINY / param%GU) * param%MU2KG + mtiny_si = (param%GMTINY / param%GU) * param%MU2KG mass_res(:) = 0.0_DP diff --git a/docs/src/symba_io.f90 b/docs/src/symba_io.f90 index 8e6420f61..713429365 100644 --- a/docs/src/symba_io.f90 +++ b/docs/src/symba_io.f90 @@ -71,8 +71,8 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms case ("FRAGMENTATION") call io_toupper(param_value) if (param_value == "YES" .or. param_value == "T") self%lfragmentation = .true. - case ("MTINY") - read(param_value, *) param%mtiny + case ("GMTINY") + read(param_value, *) param%Gmtiny case("SEED") read(param_value, *) nseeds_from_file ! Because the number of seeds can vary between compilers/systems, we need to make sure we can handle cases in which the input file has a different @@ -111,12 +111,12 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms write(*,*) "SEED: N,VAL = ",size(param%seed), param%seed(:) end if - if (self%mtiny < 0.0_DP) then - write(iomsg,*) "MTINY invalid or not set: ", self%mtiny + if (self%Gmtiny < 0.0_DP) then + write(iomsg,*) "GMTINY invalid or not set: ", self%Gmtiny iostat = -1 return else - write(*,*) "MTINY = ", self%mtiny + write(*,*) "GMTINY = ", self%Gmtiny end if if (.not.self%lclose) then @@ -167,7 +167,7 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms ! Special handling is required for writing the random number seed array as its size is not known until runtime ! For the "SEED" parameter line, the first value will be the size of the seed array and the rest will be the seed array elements write(param_name, Afmt) "PARTICLE_FILE"; write(param_value, Afmt) trim(adjustl(param%particle_file)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) - write(param_name, Afmt) "MTINY"; write(param_value, Rfmt) param%mtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "GMTINY"; write(param_value, Rfmt) param%Gmtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "FRAGMENTATION"; write(param_value, Lfmt) param%lfragmentation; write(unit, Afmt) adjustl(param_name), adjustl(param_value) if (param%lfragmentation) then write(param_name, Afmt) "SEED" diff --git a/docs/src/symba_setup.f90 b/docs/src/symba_setup.f90 index 021873a70..ab8b5543e 100644 --- a/docs/src/symba_setup.f90 +++ b/docs/src/symba_setup.f90 @@ -25,7 +25,7 @@ module subroutine symba_setup_initialize_system(self, param) call pl%sort("mass", ascending=.false.) select type(param) class is (symba_parameters) - pl%lmtiny(:) = pl%Gmass(:) > param%MTINY + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY pl%nplm = count(pl%lmtiny(:)) end select end select diff --git a/docs/src/symba_util.f90 b/docs/src/symba_util.f90 index 0517cff9a..911f85304 100644 --- a/docs/src/symba_util.f90 +++ b/docs/src/symba_util.f90 @@ -392,7 +392,7 @@ module subroutine symba_util_rearray_pl(self, system, param) ! If there are still bodies in the system, sort by mass in descending order and re-index if (pl%nbody > 0) then call pl%sort("mass", ascending=.false.) - pl%lmtiny(:) = pl%Gmass(:) > param%MTINY + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY pl%nplm = count(pl%lmtiny(:)) call pl%eucl_index() end if diff --git a/examples/symba_energy_momentum/param.disruption_headon.in b/examples/symba_energy_momentum/param.disruption_headon.in index de3c83bea..6dbe1f788 100644 --- a/examples/symba_energy_momentum/param.disruption_headon.in +++ b/examples/symba_energy_momentum/param.disruption_headon.in @@ -17,13 +17,11 @@ 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_headon.dat 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 -MTINY 1.0e-16 +GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.disruption_off_axis.in b/examples/symba_energy_momentum/param.disruption_off_axis.in index bb0915050..39303284e 100644 --- a/examples/symba_energy_momentum/param.disruption_off_axis.in +++ b/examples/symba_energy_momentum/param.disruption_off_axis.in @@ -22,7 +22,7 @@ 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 -MTINY 1.0e-16 +GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.escape.in b/examples/symba_energy_momentum/param.escape.in index 2b84eb719..99d572b75 100644 --- a/examples/symba_energy_momentum/param.escape.in +++ b/examples/symba_energy_momentum/param.escape.in @@ -19,13 +19,11 @@ CHK_RMIN 0.00465047 ! check for close solar encounters in AU CHK_RMAX 10000.0 ! discard outside of 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.escape.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT no ! Hill's sphere radii in input file -MTINY 1.0e-16 +GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.sun.in b/examples/symba_energy_momentum/param.sun.in index 65365b120..5e26e4cd3 100644 --- a/examples/symba_energy_momentum/param.sun.in +++ b/examples/symba_energy_momentum/param.sun.in @@ -19,13 +19,11 @@ CHK_RMIN 0.005 CHK_RMAX 1e2 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.escape.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT no ! Hill's sphere radii in input file -MTINY 1.0e-16 +GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.supercatastrophic_headon.in b/examples/symba_energy_momentum/param.supercatastrophic_headon.in index 19b15de7b..3ba223ad9 100644 --- a/examples/symba_energy_momentum/param.supercatastrophic_headon.in +++ b/examples/symba_energy_momentum/param.supercatastrophic_headon.in @@ -17,13 +17,11 @@ 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.supercatastrophic_headon.dat 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 -MTINY 1.0e-16 +GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in b/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in index 9cd214534..d033c76f1 100644 --- a/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in +++ b/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in @@ -23,7 +23,7 @@ ENC_OUT enc.supercatastrophic_off_axis.dat 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 -MTINY 1.0e-16 +GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py b/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py index 245f5fae0..6547c2802 100755 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py @@ -180,7 +180,7 @@ print(f'DU2M {DU2M}') print(f'TU2S {TU2S}') print(f'RHILL_PRESENT yes') -print(f'MTINY 1e-12') +print(f'GMTINY 1e-12') diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in b/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in index 3050dea4a..c69ee07f9 100644 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in @@ -30,4 +30,4 @@ MU2KG 1.988409870698051e+30 DU2M 149597870700.0 TU2S 31557600.0 RHILL_PRESENT yes -MTINY 1e-12 +GMTINY 1e-12 diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py b/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py index 6b9664542..ac55abb2b 100755 --- a/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py @@ -178,7 +178,7 @@ print(f'DU2M {DU2M}') print(f'TU2S {TU2S}') print(f'RHILL_PRESENT yes') -print(f'MTINY 1e-12') +print(f'GMTINY 1e-12') diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in index 32cfb89bd..9fb0bf743 100644 --- a/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in @@ -29,4 +29,4 @@ MU2KG 1.988409870698051e+30 DU2M 149597870700.0 TU2S 31557600.0 RHILL_PRESENT yes -MTINY 1e-12 +GMTINY 1e-12 diff --git a/examples/symba_swifter_comparison/8pl_16tp_encounters/init_cond.py b/examples/symba_swifter_comparison/8pl_16tp_encounters/init_cond.py index 707c80b51..546dcbb88 100755 --- a/examples/symba_swifter_comparison/8pl_16tp_encounters/init_cond.py +++ b/examples/symba_swifter_comparison/8pl_16tp_encounters/init_cond.py @@ -39,7 +39,7 @@ sim.param['GR'] = 'NO' sim.param['CHK_CLOSE'] = 'YES' sim.param['RHILL_PRESENT'] = 'YES' -sim.param['MTINY'] = 1.0e-12 +sim.param['GMTINY'] = 1.0e-12 sim.param['MU2KG'] = swiftest.MSun sim.param['TU2S'] = swiftest.JD2S diff --git a/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in b/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in index c72e9f4b4..f26f33592 100644 --- a/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in +++ b/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in @@ -32,4 +32,4 @@ ROTATION NO TIDES NO ENERGY NO GR NO -MTINY 1e-12 +GMTINY 1e-12 diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index c0cd7d61f..81840ba05 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -76,8 +76,8 @@ def read_swiftest_param(param_file_name, param): param['TIDES'] = param['TIDES'].upper() param['ENERGY'] = param['ENERGY'].upper() param['GR'] = param['GR'].upper() - if 'MTINY' in param: - param['MTINY'] = real2float(param['MTINY']) + if 'GMTINY' in param: + param['GMTINY'] = real2float(param['GMTINY']) except IOError: print(f"{param_file_name} not found.") return param @@ -900,7 +900,7 @@ def swift2swifter(swift_param, plname="", tpname="", conversion_questions={}): intxt = conversion_questions.get('BIG_DISCARD', None) if not intxt: - intxt = input("BIG_DISCARD: include data for all bodies > MTINY for each discard record? (y/N)> ") + intxt = input("BIG_DISCARD: include data for all bodies > GMTINY for each discard record? (y/N)> ") if intxt.upper() == 'Y': swifter_param['BIG_DISCARD'] = 'YES' else: @@ -1215,11 +1215,11 @@ def swifter2swiftest(swifter_param, plname="", tpname="", cbname="", conversion_ print(f"Cannot write to file {swiftest_param['CB_IN']}") return swifter_param - MTINY = conversion_questions.get('MTINY', None) - if not MTINY: - MTINY = input(f"Value of MTINY if this is a SyMBA simulation (enter nothing if this is not a SyMBA parameter file)> ") - if MTINY != '' and real2float(MTINY.strip()) > 0: - swiftest_param['MTINY'] = real2float(MTINY.strip()) + GMTINY = conversion_questions.get('GMTINY', None) + if not GMTINY: + GMTINY = input(f"Value of GMTINY if this is a SyMBA simulation (enter nothing if this is not a SyMBA parameter file)> ") + if GMTINY != '' and real2float(GMTINY.strip()) > 0: + swiftest_param['GMTINY'] = real2float(GMTINY.strip()) # Remove the unneeded parameters if 'C' in swiftest_param: @@ -1265,7 +1265,7 @@ def swift2swiftest(swift_param, plname="", tpname="", cbname="", conversion_ques def swiftest2swifter_param(swiftest_param, J2=0.0, J4=0.0): swifter_param = swiftest_param CBIN = swifter_param.pop("CB_IN", None) - MTINY = swifter_param.pop("MTINY", None) + GMTINY = swifter_param.pop("GMTINY", None) DISCARD_OUT = swifter_param.pop("DISCARD_OUT", None) MU2KG = swifter_param.pop("MU2KG", 1.0) DU2M = swifter_param.pop("DU2M", 1.0) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 460060183..c0cfd3c97 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -165,10 +165,9 @@ subroutine set_scale_factors() vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot ! Set scale factors - !! Because of the implied G, mass is actually G*mass with units of distance**3 / time**2 Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2))) rscale = sum(radius(:)) - mscale = sqrt(Escale * rscale) + mscale = sqrt(Escale * rscale / param%GU) vscale = sqrt(Escale / mscale) tscale = rscale / vscale Lscale = mscale * rscale * vscale @@ -346,6 +345,7 @@ subroutine calculate_system_energy(linclude_fragments) if (linclude_fragments) then ! Append the fragments if they are included ! Energy calculation requires the fragments to be in the system barcyentric frame, s tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag) + tmpsys%pl%Gmass(npl+1:npl_new) = param%GU * m_frag(1:nfrag) tmpsys%pl%radius(npl+1:npl_new) = rad_frag(1:nfrag) tmpsys%pl%xb(:,npl+1:npl_new) = xb_frag(:,1:nfrag) tmpsys%pl%vb(:,npl+1:npl_new) = vb_frag(:,1:nfrag) @@ -369,7 +369,7 @@ subroutine calculate_system_energy(linclude_fragments) class is (symba_pl) select type(param) class is (symba_parameters) - plwksp%nplm = count(plwksp%Gmass > param%mtiny / mscale) + plwksp%nplm = count(plwksp%Gmass > param%Gmtiny / mscale) end select end select call tmpsys%pl%eucl_index() @@ -381,7 +381,7 @@ subroutine calculate_system_energy(linclude_fragments) class is (symba_pl) select type(param) class is (symba_parameters) - nplm = count(pl%mass > param%mtiny) + nplm = count(pl%Gmass > param%Gmtiny) end select end select if (lk_plpl) call pl%eucl_index() @@ -836,7 +836,7 @@ end subroutine fragmentation_initialize - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, Gmtiny, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine the collisional regime of two colliding bodies. @@ -857,7 +857,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v ! Arguments integer(I4B), intent(out) :: regime real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, Gmtiny real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 real(DP), intent(out) :: Qloss !! The residual energy after the collision ! Constants @@ -931,7 +931,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v Qloss = 0.0_DP U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27 - if ((m1 < mtiny).or.(m2 < mtiny)) then + if ((m1 < Gmtiny).or.(m2 < Gmtiny)) then regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime Mlr = Mtot Mslr = 0.0_DP diff --git a/src/io/io.f90 b/src/io/io.f90 index 05313275c..ebfe2b1ef 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -493,7 +493,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) read(param_value, *) param%Ecollisions case("EUNTRACKED") read(param_value, *) param%Euntracked - case ("NPLMAX", "NTPMAX", "MTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters + case ("NPLMAX", "NTPMAX", "GMTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters case default write(iomsg,*) "Unknown parameter -> ",param_name iostat = -1 diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 25c18295a..7dc70a37f 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -470,11 +470,11 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, real(DP), intent(inout) :: Qloss end subroutine fragmentation_initialize - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, Gmtiny, Qloss) implicit none integer(I4B), intent(out) :: regime real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, Gmtiny real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 real(DP), intent(out) :: Qloss !! The residual energy after the collision end subroutine fragmentation_regime diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index ed495ecfa..dfe1c0326 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -19,7 +19,7 @@ module symba_classes type, extends(swiftest_parameters) :: symba_parameters character(STRMAX) :: particle_file = PARTICLE_OUTFILE !! Name of output particle information file - real(DP) :: MTINY = -1.0_DP !! Smallest mass that is fully gravitating + real(DP) :: GMTINY = -1.0_DP !! Smallest mass that is fully gravitating integer(I4B), dimension(:), allocatable :: seed !! Random seeds logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. contains @@ -73,9 +73,9 @@ module symba_classes type, extends(helio_pl) :: symba_pl logical, dimension(:), allocatable :: lcollision !! flag indicating whether body has merged with another this time step logical, dimension(:), allocatable :: lencounter !! flag indicating whether body is part of an encounter this time step - logical, dimension(:), allocatable :: lmtiny !! flag indicating whether this body is below the MTINY cutoff value - integer(I4B) :: nplm !! number of bodies above the MTINY limit - integer(I8B) :: nplplm !! Number of body (all massive)-body (only those above MTINY) comparisons in the flattened upper triangular matrix + logical, dimension(:), allocatable :: lmtiny !! flag indicating whether this body is below the GMTINY cutoff value + integer(I4B) :: nplm !! number of bodies above the GMTINY limit + integer(I8B) :: nplplm !! Number of body (all massive)-body (only those above GMTINY) comparisons in the flattened upper triangular matrix integer(I4B), dimension(:), allocatable :: nplenc !! number of encounters with other planets this time step integer(I4B), dimension(:), allocatable :: ntpenc !! number of encounters with test particles this time step integer(I4B), dimension(:), allocatable :: levelg !! level at which this body should be moved diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 952d59709..1910411b9 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -22,7 +22,7 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec logical, dimension(:), allocatable :: lcollision, lmask real(DP), dimension(NDIM) :: xr, vr integer(I4B) :: k - real(DP) :: rlim, mtot + real(DP) :: rlim, Gmtot logical :: isplpl if (self%nenc == 0) return @@ -55,8 +55,8 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec xr(:) = pl%xh(:, ind1(k)) - pl%xh(:, ind2(k)) vr(:) = pl%vb(:, ind1(k)) - pl%vb(:, ind2(k)) rlim = pl%radius(ind1(k)) + pl%radius(ind2(k)) - mtot = pl%Gmass(ind1(k)) + pl%Gmass(ind2(k)) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), mtot, rlim, dt, self%lvdotr(k)) + Gmtot = pl%Gmass(ind1(k)) + pl%Gmass(ind2(k)) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), Gmtot, rlim, dt, self%lvdotr(k)) end do else do concurrent(k = 1:nenc, lmask(k)) @@ -359,7 +359,7 @@ module subroutine symba_collision_make_family_pl(self, idx) p2 = pl%kin(idx(2))%parent if (p1 == p2) return ! This is a collision between to children of a shared parent. We will ignore it. - if (pl%Gmass(p1) > pl%Gmass(p2)) then + if (pl%mass(p1) > pl%mass(p2)) then index_parent = p1 index_child = p2 else @@ -449,7 +449,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) v2_si(:) = plpl_collisions%v2(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The collective density of the parent and its children Mcb_si = cb%mass * param%MU2KG - mtiny_si = (param%MTINY / param%GU) * param%MU2KG + mtiny_si = (param%GMTINY / param%GU) * param%MU2KG mass_res(:) = 0.0_DP @@ -463,8 +463,8 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) mass_res(1) = min(max(mlr, 0.0_DP), mtot) mass_res(2) = min(max(mslr, 0.0_DP), mtot) mass_res(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) - mass_res(:) = (mass_res(:) / param%MU2KG) * param%GU - Qloss = Qloss * (param%GU / param%MU2KG) * (param%TU2S / param%DU2M)**2 + mass_res(:) = (mass_res(:) / param%MU2KG) + Qloss = Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG select case (regime) case (COLLRESOLVE_REGIME_DISRUPTION) diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 4401fb5db..ba972db8b 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -71,8 +71,8 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms case ("FRAGMENTATION") call io_toupper(param_value) if (param_value == "YES" .or. param_value == "T") self%lfragmentation = .true. - case ("MTINY") - read(param_value, *) param%mtiny + case ("GMTINY") + read(param_value, *) param%Gmtiny case("SEED") read(param_value, *) nseeds_from_file ! Because the number of seeds can vary between compilers/systems, we need to make sure we can handle cases in which the input file has a different @@ -111,12 +111,12 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms write(*,*) "SEED: N,VAL = ",size(param%seed), param%seed(:) end if - if (self%mtiny < 0.0_DP) then - write(iomsg,*) "MTINY invalid or not set: ", self%mtiny + if (self%Gmtiny < 0.0_DP) then + write(iomsg,*) "GMTINY invalid or not set: ", self%Gmtiny iostat = -1 return else - write(*,*) "MTINY = ", self%mtiny + write(*,*) "GMTINY = ", self%Gmtiny end if if (.not.self%lclose) then @@ -167,7 +167,7 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms ! Special handling is required for writing the random number seed array as its size is not known until runtime ! For the "SEED" parameter line, the first value will be the size of the seed array and the rest will be the seed array elements write(param_name, Afmt) "PARTICLE_FILE"; write(param_value, Afmt) trim(adjustl(param%particle_file)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) - write(param_name, Afmt) "MTINY"; write(param_value, Rfmt) param%mtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "GMTINY"; write(param_value, Rfmt) param%Gmtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "FRAGMENTATION"; write(param_value, Lfmt) param%lfragmentation; write(unit, Afmt) adjustl(param_name), adjustl(param_value) if (param%lfragmentation) then write(param_name, Afmt) "SEED" diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index 021873a70..ab8b5543e 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -25,7 +25,7 @@ module subroutine symba_setup_initialize_system(self, param) call pl%sort("mass", ascending=.false.) select type(param) class is (symba_parameters) - pl%lmtiny(:) = pl%Gmass(:) > param%MTINY + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY pl%nplm = count(pl%lmtiny(:)) end select end select diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 2ebc11ebb..2ab088d02 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -394,7 +394,7 @@ module subroutine symba_util_rearray_pl(self, system, param) ! If there are still bodies in the system, sort by mass in descending order and re-index if (pl%nbody > 0) then call pl%sort("mass", ascending=.false.) - pl%lmtiny(:) = pl%Gmass(:) > param%MTINY + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY pl%nplm = count(pl%lmtiny(:)) call pl%eucl_index() end if diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index b004f84b7..fd331bcc3 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -83,7 +83,7 @@ module subroutine util_get_energy_momentum_system(self, param) ! Do the potential energy between pairs of massive bodies do k = 1, pl%nplpl associate(ik => pl%k_plpl(1, k), jk => pl%k_plpl(2, k)) - pepl(k) = -pl%mass(ik) * pl%mass(jk) / norm2(pl%xh(:, jk) - pl%xh(:, ik)) + pepl(k) = -param%GU * pl%mass(ik) * pl%mass(jk) / norm2(pl%xh(:, jk) - pl%xh(:, ik)) lstatpl(k) = (lstatus(ik) .and. lstatus(jk)) end associate end do From 9214b2ab7ccd9ea89124bba5f7bec066cb1f4dda Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 6 Aug 2021 15:58:39 -0400 Subject: [PATCH 3/7] Refactored MTINY into GMTINY to be consistent in distinguishing G*mass and mass terms --- .../symba_energy_momentum/param.supercatastrophic_off_axis.in | 2 -- 1 file changed, 2 deletions(-) diff --git a/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in b/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in index d033c76f1..49b8b0dd7 100644 --- a/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in +++ b/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in @@ -17,8 +17,6 @@ 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.supercatastrophic_off_axis.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded From 8d5455c231b0387381d26d49de8a4d0324c4564f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 6 Aug 2021 16:03:58 -0400 Subject: [PATCH 4/7] Undid some damage due to previous refactoring as this subroutine actually wants an mtiny value not Gmtiny --- src/fragmentation/fragmentation.f90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index c0cfd3c97..353867c0b 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -835,8 +835,7 @@ end subroutine restructure_failed_fragments end subroutine fragmentation_initialize - - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, Gmtiny, Qloss) + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine the collisional regime of two colliding bodies. @@ -931,7 +930,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v Qloss = 0.0_DP U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27 - if ((m1 < Gmtiny).or.(m2 < Gmtiny)) then + if ((m1 < mtiny).or.(m2 < mtiny)) then regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime Mlr = Mtot Mslr = 0.0_DP From 71d17198492e4b0381aef5e7f94b3b96ce1478b9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 6 Aug 2021 16:11:55 -0400 Subject: [PATCH 5/7] Undid some damage due to previous refactoring as this subroutine actually wants an mtiny value not Gmtiny. Fixed formatting of fragmentation interfaces --- src/fragmentation/fragmentation.f90 | 24 ++++++++++---------- src/modules/swiftest_classes.f90 | 34 ++++++++++++++--------------- 2 files changed, 29 insertions(+), 29 deletions(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 353867c0b..60572e99d 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -10,17 +10,17 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, use, intrinsic :: ieee_exceptions implicit none ! Arguments - class(swiftest_nbody_system), intent(inout) :: system - class(swiftest_parameters), intent(in) :: param - integer(I4B), dimension(:), intent(in) :: family - real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip - real(DP), dimension(:), intent(inout) :: mass, radius - integer(I4B), intent(inout) :: nfrag - real(DP), dimension(:), allocatable, intent(inout) :: m_frag, rad_frag - real(DP), dimension(:,:), allocatable, intent(inout) :: Ip_frag - real(DP), dimension(:,:), allocatable, intent(inout) :: xb_frag, vb_frag, rot_frag - logical, intent(out) :: lfailure ! Answers the question: Should this have been a merger instead? - real(DP), intent(inout) :: Qloss + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + integer(I4B), dimension(:), intent(in) :: family !! Index of bodies involved in the collision + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Two-body equivalent position, vector, spin momentum, and rotational inertia values for the collision + real(DP), dimension(:), intent(inout) :: mass, radius !! Two-body equivalent mass and radii for the bodies in the collision + integer(I4B), intent(inout) :: nfrag !! Number of fragments to generate + real(DP), dimension(:), allocatable, intent(inout) :: m_frag, rad_frag !! Distribution of fragment mass and radii + real(DP), dimension(:,:), allocatable, intent(inout) :: Ip_frag !! Fragment rotational inertia vectors + real(DP), dimension(:,:), allocatable, intent(inout) :: xb_frag, vb_frag, rot_frag !! Fragment barycentric position, barycentric velocity, and rotation vectors + real(DP), intent(inout) :: Qloss !! Energy lost during the collision + logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? ! Internals real(DP) :: mscale, rscale, vscale, tscale, Lscale, Escale ! Scale factors that reduce quantities to O(~1) in the collisional system real(DP) :: mtot @@ -856,7 +856,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v ! Arguments integer(I4B), intent(out) :: regime real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, Gmtiny + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 real(DP), intent(out) :: Qloss !! The residual energy after the collision ! Constants diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 7dc70a37f..e8d273f1c 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -457,26 +457,26 @@ module subroutine eucl_dist_index_plpl(self) module subroutine 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) implicit none - class(swiftest_nbody_system), intent(inout) :: system - class(swiftest_parameters), intent(in) :: param - integer(I4B), dimension(:), intent(in) :: family - real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip - real(DP), dimension(:), intent(inout) :: mass, radius - integer(I4B), intent(inout) :: nfrag - real(DP), dimension(:), allocatable, intent(inout) :: m_frag, rad_frag - real(DP), dimension(:,:), allocatable, intent(inout) :: Ip_frag - real(DP), dimension(:,:), allocatable, intent(inout) :: xb_frag, vb_frag, rot_frag - logical, intent(out) :: lfailure ! Answers the question: Should this have been a merger instead? - real(DP), intent(inout) :: Qloss + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + integer(I4B), dimension(:), intent(in) :: family !! Index of bodies involved in the collision + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Two-body equivalent position, vector, spin momentum, and rotational inertia values for the collision + real(DP), dimension(:), intent(inout) :: mass, radius !! Two-body equivalent mass and radii for the bodies in the collision + integer(I4B), intent(inout) :: nfrag !! Number of fragments to generate + real(DP), dimension(:), allocatable, intent(inout) :: m_frag, rad_frag !! Distribution of fragment mass and radii + real(DP), dimension(:,:), allocatable, intent(inout) :: Ip_frag !! Fragment rotational inertia vectors + real(DP), dimension(:,:), allocatable, intent(inout) :: xb_frag, vb_frag, rot_frag !! Fragment barycentric position, barycentric velocity, and rotation vectors + real(DP), intent(inout) :: Qloss !! Energy lost during the collision + logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? end subroutine fragmentation_initialize - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, Gmtiny, Qloss) + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) implicit none - integer(I4B), intent(out) :: regime - real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, Gmtiny - real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 - real(DP), intent(out) :: Qloss !! The residual energy after the collision + integer(I4B), intent(out) :: regime + real(DP), intent(out) :: Mlr, Mslr + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny + real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 + real(DP), intent(out) :: Qloss !! Energy lost during the collision end subroutine fragmentation_regime module pure subroutine gr_kick_getaccb_ns_body(self, system, param) From 77ee9ccdb892192c12aebf99ee7c6efdcce44323 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 6 Aug 2021 16:33:13 -0400 Subject: [PATCH 6/7] Changed id validation code to by type-bound to the system object, and include cb ids in the validation --- Makefile.Defines | 4 ++-- src/fragmentation/fragmentation.f90 | 8 ++++++-- src/modules/swiftest_classes.f90 | 9 +++++---- src/setup/setup.f90 | 4 +--- src/util/util_valid.f90 | 20 +++++++++++--------- 5 files changed, 25 insertions(+), 20 deletions(-) diff --git a/Makefile.Defines b/Makefile.Defines index 291f2c604..9fe6c3cbb 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -70,8 +70,8 @@ FFLAGS = $(IDEBUG) $(HEAPARR) FORTRAN = ifort #AR = xiar -#FORTRAN = gfortran -#FFLAGS = -ffree-line-length-none $(GDEBUG) $(GMEM) +FORTRAN = gfortran +FFLAGS = -ffree-line-length-none $(GDEBUG) $(GMEM) AR = ar # DO NOT include in CFLAGS the "-c" option to compile object only diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 60572e99d..99daff3b7 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -336,8 +336,12 @@ subroutine calculate_system_energy(linclude_fragments) call setup_construct_system(tmpsys, param) deallocate(tmpsys%cb) allocate(tmpsys%cb, source=cb) - allocate(ltmp(npl)) - ltmp(:) = .true. + allocate(ltmp, mold=pl%ldiscard) + ltmp(:) = .false. + ltmp(1:npl) = .true. + write(*,*) 'npl : ',npl + write(*,*) 'npl_new: ',npl_new + write(*,*) 'ltmp: ',ltmp call tmpsys%pl%setup(npl_new, param) call tmpsys%pl%fill(pl, ltmp) deallocate(ltmp) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index e8d273f1c..7f4fd140e 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -306,6 +306,7 @@ module swiftest_classes procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. procedure :: set_msys => util_set_msys !! Sets the value of msys from the masses of system bodies. procedure :: get_energy_and_momentum => util_get_energy_momentum_system !! Calculates the total system energy and momentum + procedure :: validate_ids => util_valid_id_system !! Validate the numerical ids passed to the system and save the maximum value end type swiftest_nbody_system type :: swiftest_encounter @@ -1292,11 +1293,11 @@ module subroutine util_spill_tp(self, discards, lspill_list, ldestructive) logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not end subroutine util_spill_tp - module subroutine util_valid(pl, tp) + module subroutine util_valid_id_system(self, param) implicit none - class(swiftest_pl), intent(in) :: pl - class(swiftest_tp), intent(in) :: tp - end subroutine util_valid + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine util_valid_id_system module subroutine util_version() implicit none diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 9346e8c12..8f96c48a1 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -127,14 +127,12 @@ module subroutine setup_initialize_system(self, param) call self%cb%initialize(param) call self%pl%initialize(param) call self%tp%initialize(param) - call util_valid(self%pl, self%tp) - self%maxid = maxval([self%pl%id(:), self%tp%id(:)]) + call self%validate_ids(param) call self%set_msys() call self%pl%set_mu(self%cb) call self%tp%set_mu(self%cb) call self%pl%eucl_index() if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb) - !if (param%lfirstenergy) then return end subroutine setup_initialize_system diff --git a/src/util/util_valid.f90 b/src/util/util_valid.f90 index c5923b38e..f05c81f35 100644 --- a/src/util/util_valid.f90 +++ b/src/util/util_valid.f90 @@ -2,7 +2,7 @@ use swiftest contains - module subroutine util_valid(pl, tp) + module subroutine util_valid_id_system(self, param) !! author: David A. Minton !! !! Validate massive body and test particle ids @@ -11,31 +11,33 @@ module subroutine util_valid(pl, tp) !! Adapted from David E. Kaufmann's Swifter routine: util_valid.f90 implicit none ! Arguments - class(swiftest_pl), intent(in) :: pl - class(swiftest_tp), intent(in) :: tp + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i integer(I4B), dimension(:), allocatable :: idarr - associate(npl => pl%nbody, ntp => tp%nbody) - allocate(idarr(npl+ntp)) + associate(cb => self%cb, pl => self%pl, npl => self%pl%nbody, tp => self%tp, ntp => self%tp%nbody) + allocate(idarr(1+npl+ntp)) + idarr(1) = cb%id do i = 1, npl - idarr(i) = pl%id(i) + idarr(1+i) = pl%id(i) end do do i = 1, ntp - idarr(npl+i) = tp%id(i) + idarr(1+npl+i) = tp%id(i) end do call util_sort(idarr) - do i = 1, npl + ntp - 1 + do i = 1, npl + ntp if (idarr(i) == idarr(i+1)) then write(*, *) "Swiftest error:" write(*, *) " more than one body/particle has id = ", idarr(i) call util_exit(FAILURE) end if end do + self%maxid = maxval(idarr) end associate return - end subroutine util_valid + end subroutine util_valid_id_system end submodule s_util_valid From adf3f1a94e62164dbdae7f28121b3e7ef69d3f87 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 6 Aug 2021 16:39:47 -0400 Subject: [PATCH 7/7] Fixed bad mask size on fill operation --- Makefile.Defines | 4 ++-- src/fragmentation/fragmentation.f90 | 7 +------ 2 files changed, 3 insertions(+), 8 deletions(-) diff --git a/Makefile.Defines b/Makefile.Defines index 9fe6c3cbb..291f2c604 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -70,8 +70,8 @@ FFLAGS = $(IDEBUG) $(HEAPARR) FORTRAN = ifort #AR = xiar -FORTRAN = gfortran -FFLAGS = -ffree-line-length-none $(GDEBUG) $(GMEM) +#FORTRAN = gfortran +#FFLAGS = -ffree-line-length-none $(GDEBUG) $(GMEM) AR = ar # DO NOT include in CFLAGS the "-c" option to compile object only diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 99daff3b7..ad28edf6c 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -336,15 +336,11 @@ subroutine calculate_system_energy(linclude_fragments) call setup_construct_system(tmpsys, param) deallocate(tmpsys%cb) allocate(tmpsys%cb, source=cb) - allocate(ltmp, mold=pl%ldiscard) + allocate(ltmp(npl_new)) ltmp(:) = .false. ltmp(1:npl) = .true. - write(*,*) 'npl : ',npl - write(*,*) 'npl_new: ',npl_new - write(*,*) 'ltmp: ',ltmp call tmpsys%pl%setup(npl_new, param) call tmpsys%pl%fill(pl, ltmp) - deallocate(ltmp) if (linclude_fragments) then ! Append the fragments if they are included ! Energy calculation requires the fragments to be in the system barcyentric frame, s @@ -359,7 +355,6 @@ subroutine calculate_system_energy(linclude_fragments) tmpsys%pl%rot(:,npl+1:npl_new) = rot_frag(:,1:nfrag) end if call tmpsys%pl%b2h(tmpsys%cb) - allocate(ltmp(npl_new)) ltmp(1:npl) = lexclude(1:npl) ltmp(npl+1:npl_new) = .false. call move_alloc(ltmp, lexclude)