diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 315b098a8..fbbfa15fe 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -70,7 +70,7 @@ module rmvs_classes procedure :: encounter_check => rmvs_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body procedure :: accel => rmvs_kick_getacch_tp !! Calculates either the standard or modified version of the acceleration depending if the !! if the test particle is undergoing a close encounter or not - procedure :: setup => rmvs_setup_tp !! Constructor method - Allocates space for number of particles + procedure :: setup => rmvs_setup_tp !! Constructor method - Allocates space for the input number of bodiess procedure :: append => rmvs_util_append_tp !! Appends elements from one structure to another procedure :: fill => rmvs_util_fill_tp !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) procedure :: resize => rmvs_util_resize_tp !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small. @@ -93,7 +93,7 @@ module rmvs_classes class(rmvs_nbody_system), dimension(:), allocatable :: planetocentric !! Planetocentric version of the massive body objects (one for each massive body) logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations contains - procedure :: setup => rmvs_setup_pl !! Constructor method - Allocates space for number of particles + procedure :: setup => rmvs_setup_pl !! Constructor method - Allocates space for the input number of bodiess procedure :: append => rmvs_util_append_pl !! Appends elements from one structure to another procedure :: fill => rmvs_util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) procedure :: resize => rmvs_util_resize_pl !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small. diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 45b5acfbc..b0daf9496 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -91,16 +91,24 @@ module symba_classes procedure :: drift => symba_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level procedure :: encounter_check => symba_encounter_check_pl !! Checks if massive bodies are going through close encounters with each other procedure :: accel => symba_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies - procedure :: setup => symba_setup_pl !! Constructor method - Allocates space for number of particle + procedure :: setup => symba_setup_pl !! Constructor method - Allocates space for the input number of bodies procedure :: append => symba_util_append_pl !! Appends elements from one structure to another procedure :: fill => symba_util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) procedure :: get_peri => symba_util_peri_pl !! Determine system pericenter passages for massive bodies - procedure :: resize => symba_util_resize_pl !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small. + procedure :: resize => symba_util_resize_pl !! Checks the current size of a SyMBA massive body against the requested size and resizes it if it is too small. procedure :: sort => symba_util_sort_pl !! Sorts body arrays by a sortable componen procedure :: rearrange => symba_util_sort_rearrange_pl !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods procedure :: spill => symba_util_spill_pl !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) end type symba_pl + type, extends(symba_pl) :: symba_merger + integer(I4B), dimension(:), allocatable :: ncomp + contains + procedure :: append => symba_util_append_merger !! Appends elements from one structure to another + procedure :: resize => symba_util_resize_merger !! Checks the current size of a SyMBA merger list against the requested size and resizes it if it is too small. + procedure :: setup => symba_setup_merger !! Constructor method - Allocates space for the input number of bodies + end type symba_merger + !******************************************************************************************************************************** ! symba_tp class definitions and method interfaces !******************************************************************************************************************************* @@ -113,7 +121,7 @@ module symba_classes procedure :: drift => symba_drift_tp !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level procedure :: encounter_check => symba_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body procedure :: accel => symba_kick_getacch_tp !! Compute heliocentric accelerations of test particles - procedure :: setup => symba_setup_tp !! Constructor method - Allocates space for number of particle + procedure :: setup => symba_setup_tp !! Constructor method - Allocates space for the input number of bodies procedure :: append => symba_util_append_tp !! Appends elements from one structure to another procedure :: fill => symba_util_fill_tp !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) procedure :: resize => symba_util_resize_tp !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small. @@ -151,8 +159,8 @@ module symba_classes ! symba_nbody_system class definitions and method interfaces !******************************************************************************************************************************** type, extends(helio_nbody_system) :: symba_nbody_system - class(symba_pl), allocatable :: mergeadd_list !! List of added bodies in mergers or collisions - class(symba_pl), allocatable :: mergesub_list !! List of subtracted bodies in mergers or collisions + 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_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 @@ -266,6 +274,16 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc logical :: lany_encounter !! Returns true if there is at least one close encounter end function symba_encounter_check_tp + module function symba_fragmentation_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) result(status) + implicit none + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision + real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(in) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collisio + integer(I4B) :: status !! Status flag assigned to this outcome + end function symba_fragmentation_casemerge + module subroutine symba_io_write_discard(self, param) use swiftest_classes, only : swiftest_parameters implicit none @@ -357,11 +375,26 @@ module subroutine symba_io_write_frame_info(self, iu, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine symba_io_write_frame_info + module subroutine symba_setup_initialize_system(self, param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_nbody_system), intent(inout) :: self !! SyMBA system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine symba_setup_initialize_system + + module subroutine symba_setup_merger(self, n, param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_merger), intent(inout) :: self !! SyMBA merger list object + integer(I4B), intent(in) :: n !! Number of particles to allocate space for + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine symba_setup_merger + module subroutine symba_setup_pl(self, n, param) use swiftest_classes, only : swiftest_parameters implicit none - class(symba_pl), intent(inout) :: self !! SyMBA massive body object - integer(I4B), intent(in) :: n !! Number of particles to allocate space for + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + integer(I4B), intent(in) :: n !! Number of particles to allocate space for class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine symba_setup_pl @@ -371,19 +404,6 @@ module subroutine symba_setup_pltpenc(self,n) integer(I4B), intent(in) :: n !! Number of encounters to allocate space for end subroutine symba_setup_pltpenc - module subroutine symba_setup_plplenc(self,n) - implicit none - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-tp encounter structure - integer(I4B), intent(in) :: n !! Number of encounters to allocate space for - end subroutine symba_setup_plplenc - - module subroutine symba_setup_initialize_system(self, param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_nbody_system), intent(inout) :: self !! SyMBA system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine symba_setup_initialize_system - module subroutine symba_setup_tp(self, n, param) use swiftest_classes, only : swiftest_parameters implicit none @@ -448,6 +468,14 @@ end subroutine symba_util_append_arr_kin end interface interface + module subroutine symba_util_append_merger(self, source, lsource_mask) + use swiftest_classes, only : swiftest_body + implicit none + class(symba_merger), intent(inout) :: self !! SyMBA massive body object + class(swiftest_body), intent(in) :: source !! Source object to append + logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + end subroutine symba_util_append_merger + module subroutine symba_util_append_pl(self, source, lsource_mask) use swiftest_classes, only : swiftest_body implicit none @@ -522,6 +550,12 @@ end subroutine symba_util_resize_arr_kin end interface interface + module subroutine symba_util_resize_merger(self, nnew) + implicit none + class(symba_merger), intent(inout) :: self !! SyMBA merger list object + integer(I4B), intent(in) :: nnew !! New size neded + end subroutine symba_util_resize_merger + module subroutine symba_util_resize_pl(self, nnew) implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index e581e52b1..a79f52bca 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -46,7 +46,7 @@ module whm_classes procedure :: sort => whm_util_sort_pl !! Sort a WHM massive body object in-place. procedure :: rearrange => whm_util_sort_rearrange_pl !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods procedure :: spill => whm_util_spill_pl !!"Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) - procedure :: setup => whm_setup_pl !! Constructor method - Allocates space for number of particles + procedure :: setup => whm_setup_pl !! Constructor method - Allocates space for the input number of bodiess procedure :: step => whm_step_pl !! Steps the body forward one stepsize end type whm_pl diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index ca5f38c6e..edb641907 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_pl :: system%mergeadd_list) - allocate(symba_pl :: system%mergesub_list) + allocate(symba_merger :: system%mergeadd_list) + allocate(symba_merger :: system%mergesub_list) allocate(symba_plplenc :: system%plplenc_list) allocate(symba_pltpenc :: system%pltpenc_list) end select diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index efae0ba39..0f8b4d519 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -208,8 +208,8 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v fam_size = 2 + sum(nchild(:)) allocate(family(fam_size)) family = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] - fam_size = count(pl%status(family(:)) == ACTIVE) - family = pack(family(:), pl%status(family(:)) == ACTIVE) + fam_size = count(pl%status(family(:)) == COLLISION) + family = pack(family(:), pl%status(family(:)) == COLLISION) L_spin(:,:) = 0.0_DP Ip(:,:) = 0.0_DP @@ -423,12 +423,12 @@ module subroutine symba_collision_resolve_mergers(self, system, param) class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions ! Internals - integer(I4B) :: i - logical :: lgoodcollision - integer(I4B), dimension(:), allocatable :: family !! List of indices of all bodies inovlved in the collision - integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - real(DP), dimension(NDIM,2) :: x, v, L_spin, Ip !! Output values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(2) :: mass, radius !! Output values that represent a 2-body equivalent of a possibly 2+ body collision + integer(I4B), dimension(:), allocatable :: family !! List of indices of all bodies inovlved in the collision + integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision + real(DP), dimension(NDIM,2) :: x, v, L_spin, Ip !! Output values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(2) :: mass, radius !! Output values that represent a 2-body equivalent of a possibly 2+ body collision + logical :: lgoodcollision + integer(I4B) :: i, status associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) select type(pl => system%pl) @@ -439,7 +439,7 @@ module subroutine symba_collision_resolve_mergers(self, system, param) lgoodcollision = symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) if (.not. lgoodcollision) cycle if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved - !call symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) + status = symba_fragmentation_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) end do end select end associate diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 2fa789d46..1c5b4a732 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -153,6 +153,7 @@ module subroutine symba_discard_pl(self, system, param) call symba_discard_nonplpl(self, system, param) call plplenc_list%scrub_non_collision(system, param) if (plplenc_list%nenc == 0) return ! No collisions to resolve + write(*, *) "Collision detected at time t = ",param%t call pl%h2b(system%cb) if (param%lfragmentation) then diff --git a/src/symba/symba_fragmentation.f90 b/src/symba/symba_fragmentation.f90 new file mode 100644 index 000000000..ccc79892f --- /dev/null +++ b/src/symba/symba_fragmentation.f90 @@ -0,0 +1,135 @@ +submodule (symba_classes) s_symba_fragmentation + use swiftest +contains + + module function symba_fragmentation_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) result(status) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Merge planets. + !! + !! Adapted from David E. Kaufmann's Swifter routines symba_merge_pl.f90 and symba_discard_merge_pl.f90 + !! + !! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision + real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(in) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + ! Result + integer(I4B) :: status !! Status flag assigned to this outcome + ! Internals + integer(I4B) :: i, j, ibiggest, nfamily, nstart, nend + real(DP) :: mass_new, radius_new, volume_new, pe + real(DP), dimension(NDIM) :: xcom, vcom, xc, vc, xcrossv + real(DP), dimension(2) :: vol + real(DP), dimension(NDIM) :: L_orb_old, L_spin_old + real(DP), dimension(NDIM) :: L_spin_new, rot_new, Ip_new + logical, dimension(system%pl%nbody) :: lmask + class(symba_pl), allocatable :: plnew + + 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 + + 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) + + ! 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%Gmass(i) * pl%Gmass(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 = 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 + + ! Create the new merged body + allocate(plnew, mold=pl) + call plnew%setup(1, param) + + ! 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%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 + pl%Ip(:,1) = Ip_new(:) + pl%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 + + ! 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 associate + end select + + return + + end function symba_fragmentation_casemerge + +end submodule s_symba_fragmentation diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index f35e45408..d3091155d 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -215,7 +215,7 @@ module subroutine symba_io_write_discard(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals integer(I4B), parameter :: LUN = 40 - integer(I4B) :: i, ierr + integer(I4B) :: iadd, isub, j, ierr, nsub, nadd logical, save :: lfirst = .true. real(DP), dimension(:,:), allocatable :: vh character(*), parameter :: HDRFMT = '(E23.16, 1X, I8, 1X, L1)' @@ -245,33 +245,35 @@ module subroutine symba_io_write_discard(self, param) end if write(LUN, HDRFMT) param%t, mergesub_list%nbody, param%lbig_discard - do i = 1, mergesub_list%nbody - write(LUN, NAMEFMT) SUB, mergesub_list%id(i), mergesub_list%status(i) - write(LUN, VECFMT) mergesub_list%xh(1, i), mergesub_list%xh(2, i), mergesub_list%xh(3, i) - write(LUN, VECFMT) mergesub_list%vh(1, i), mergesub_list%vh(2, i), mergesub_list%vh(3, i) + 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) SUB, 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 + end do end do - ! This is incomplete until the mergeadd_list methods are completed - ! if (param%lbig_discard) then - ! if (param%lgr) then - ! allocate(pltemp, source = pl) - ! call pltemp%pv2v(param) - ! allocate(vh, source = pltemp%vh) - ! deallocate(pltemp) - ! else - ! allocate(vh, source = pl%vh) - ! end if - - ! write(LUN, NPLFMT) npl - ! do i = 1, npl - ! write(LUN, PLNAMEFMT) pl%id(i), pl%Gmass(i), pl%radius(i) - ! write(LUN, VECFMT) pl%xh(1, i), pl%xh(2, i), pl%xh(3, i) - ! write(LUN, VECFMT) vh(1, i), vh(2, i), vh(3, i) - ! end do - ! deallocate(vh) - ! end if - close(LUN) - end associate + close(LUN) + end associate return end subroutine symba_io_write_discard diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index 524420609..021873a70 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -35,6 +35,32 @@ module subroutine symba_setup_initialize_system(self, param) end subroutine symba_setup_initialize_system + module subroutine symba_setup_merger(self, n, param) + !! author: David A. Minton + !! + !! Allocate SyMBA test particle structure + !! + !! Equivalent in functionality to David E. Kaufmann's Swifter routine symba_setup.f90 + implicit none + ! Arguments + class(symba_merger), intent(inout) :: self !! SyMBA merger list object + integer(I4B), intent(in) :: n !! Number of particles to allocate space for + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter + ! Internals + integer(I4B) :: i + + !> Call allocation method for parent class. In this case, helio_pl does not have its own setup method so we use the base method for swiftest_pl + call symba_setup_pl(self, n, param) + if (n <= 0) return + + if (allocated(self%ncomp)) deallocate(self%ncomp) + allocate(self%ncomp(n)) + self%ncomp(:) = 0 + + return + end subroutine symba_setup_merger + + module subroutine symba_setup_pl(self, n, param) !! author: David A. Minton !! diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 56eaacc2c..7063e0910 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -116,6 +116,37 @@ module subroutine symba_util_append_pl(self, source, lsource_mask) end subroutine symba_util_append_pl + module subroutine symba_util_append_merger(self, source, lsource_mask) + !! author: David A. Minton + !! + !! Append components from one massive body object to another. + !! This method will automatically resize the destination body if it is too small + implicit none + ! Arguments + class(symba_merger), intent(inout) :: self !! SyMBA massive body object + class(swiftest_body), intent(in) :: source !! Source object to append + logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + ! Internals + integer(I4B), dimension(:), allocatable :: ncomp_tmp !! Temporary placeholder for ncomp incase we are appending a symba_pl object to a symba_merger + + select type(source) + class is (symba_merger) + call symba_util_append_pl(self, source, lsource_mask) + call util_append(self%ncomp, source%ncomp, lsource_mask) + class is (symba_pl) + call symba_util_append_pl(self, source, lsource_mask) + allocate(ncomp_tmp, mold=source%id) + ncomp_tmp(:) = 0 + call util_append(self%ncomp, ncomp_tmp, lsource_mask) + class default + write(*,*) "Invalid object passed to the append method. Source must be of class symba_pl or its descendents!" + call util_exit(FAILURE) + end select + + return + end subroutine symba_util_append_merger + + module subroutine symba_util_append_tp(self, source, lsource_mask) !! author: David A. Minton !! @@ -406,10 +437,27 @@ module subroutine symba_util_resize_arr_kin(arr, nnew) end subroutine symba_util_resize_arr_kin + module subroutine symba_util_resize_merger(self, nnew) + !! author: David A. Minton + !! + !! Checks the current size of a SyMBA merger list against the requested size and resizes it if it is too small. + implicit none + ! Arguments + class(symba_merger), intent(inout) :: self !! SyMBA massive body object + integer(I4B), intent(in) :: nnew !! New size neded + + call symba_util_resize_pl(self, nnew) + + call util_resize(self%ncomp, nnew) + + return + end subroutine symba_util_resize_merger + + module subroutine symba_util_resize_pl(self, nnew) !! author: David A. Minton !! - !! Checks the current size of a massive body object against the requested size and resizes it if it is too small. + !! Checks the current size of a SyMBA massive body object against the requested size and resizes it if it is too small. implicit none ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA massive body object