diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index c9939d7db..8ba00c7cf 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -98,6 +98,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, L_frag_budget(:) = L_frag_orb(:) + L_frag_spin(:) f_spin = norm2(L_frag_spin(:)) / norm2(L_frag_budget(:)) + call reset_fragments() call define_coordinate_system() call construct_temporary_system() @@ -110,7 +111,6 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, ke_tot_deficit = 0.0_DP do while (try < MAXTRY) lfailure = .false. - call reset_fragments() call set_fragment_position_vectors() @@ -145,6 +145,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, if (.not.lfailure) exit call restructure_failed_fragments() + call reset_fragments() try = try + 1 end do call restore_scale_factors() @@ -326,8 +327,8 @@ subroutine define_coordinate_system() ! The cross product of the y- by z-axis will give us the x-axis x_col_unit(:) = y_col_unit(:) .cross. z_col_unit(:) + if (.not.any(x_frag(:,:) > 0.0_DP)) return rmag(:) = .mag. x_frag(:,:) - if (.not.any(rmag(:) > 0.0_DP)) return call random_number(L_sigma(:,:)) ! Randomize the tangential velocity direction. This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, ! otherwise we can get an ill-conditioned system diff --git a/src/helio/helio_coord.f90 b/src/helio/helio_coord.f90 deleted file mode 100644 index f40781810..000000000 --- a/src/helio/helio_coord.f90 +++ /dev/null @@ -1,113 +0,0 @@ -submodule (helio_classes) s_helio_coord - use swiftest -contains - - module subroutine helio_coord_vb2vh_pl(self, cb) - !! author: David A. Minton - !! - !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only) - !! - !! Adapted from David E. Kaufmann's Swifter routine coord_vb2vh.f90 - !! Adapted from Hal Levison's Swift routine coord_vb2vh.f - implicit none - ! Arguments - class(helio_pl), intent(inout) :: self !! Helio massive body object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - ! Internals - integer(I4B) :: i - - if (self%nbody == 0) return - - associate(pl => self, npl => self%nbody) - do i = 1, NDIM - cb%vb(i) = -sum(pl%Gmass(1:npl) * pl%vb(i, 1:npl)) / cb%Gmass - pl%vh(i, 1:npl) = pl%vb(i, 1:npl) - cb%vb(i) - end do - end associate - - return - end subroutine helio_coord_vb2vh_pl - - - module subroutine helio_coord_vb2vh_tp(self, vbcb) - !! author: David A. Minton - !! - !! Convert test particles from barycentric to heliocentric coordinates (velocity only) - !! - !! Adapted from David E. Kaufmann's Swifter routine coord_vb2vh_tp.f90 - !! Adapted from Hal Levison's Swift routine coord_vb2h_tp.f - implicit none - ! Arguments - class(helio_tp), intent(inout) :: self !! Helio massive body object - real(DP), dimension(:), intent(in) :: vbcb !! Barycentric velocity of the central body - - if (self%nbody == 0) return - - associate(tp => self, ntp => self%nbody) - where (tp%lmask(1:ntp)) - tp%vh(1, 1:ntp) = tp%vb(1, 1:ntp) - vbcb(1) - tp%vh(2, 1:ntp) = tp%vb(2, 1:ntp) - vbcb(2) - tp%vh(3, 1:ntp) = tp%vb(3, 1:ntp) - vbcb(3) - end where - end associate - - return - end subroutine helio_coord_vb2vh_tp - - - module subroutine helio_coord_vh2vb_pl(self, cb) - !! author: David A. Minton - !! - !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only) - !! - !! Adapted from David E. Kaufmann's Swifter routine coord_vh2vb.f90 - !! Adapted from Hal Levison's Swift routine coord_vh2b.f - implicit none - ! Arguments - class(helio_pl), intent(inout) :: self !! Helio massive body object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - ! Internals - integer(I4B) :: i - real(DP) :: Gmtot - - if (self%nbody == 0) return - - associate(pl => self, npl => self%nbody) - Gmtot = cb%Gmass + sum(pl%Gmass(1:npl)) - do i = 1, NDIM - cb%vb(i) = -sum(pl%Gmass(1:npl) * pl%vh(i, 1:npl)) / Gmtot - pl%vb(i, 1:npl) = pl%vh(i, 1:npl) + cb%vb(i) - end do - end associate - - return - end subroutine helio_coord_vh2vb_pl - - - module subroutine helio_coord_vh2vb_tp(self, vbcb) - !! author: David A. Minton - !! - !! Convert test particles from heliocentric to barycentric coordinates (velocity only) - !! - !! Adapted from David E. Kaufmann's Swifter routine coord_vh2vb_tp.f90 - !! Adapted from Hal Levison's Swift routine coord_vh2b_tp.f - implicit none - ! Arguments - class(helio_tp), intent(inout) :: self !! Helio massive body object - real(DP), dimension(:), intent(in) :: vbcb !! Barycentric velocity of the central body - - if (self%nbody == 0) return - - associate(tp => self, ntp => self%nbody) - where (tp%lmask(1:ntp)) - tp%vb(1, 1:ntp) = tp%vh(1, 1:ntp) + vbcb(1) - tp%vb(2, 1:ntp) = tp%vh(2, 1:ntp) + vbcb(2) - tp%vb(3, 1:ntp) = tp%vh(3, 1:ntp) + vbcb(3) - end where - end associate - - return - end subroutine helio_coord_vh2vb_tp - -end submodule s_helio_coord - diff --git a/src/helio/helio_setup.f90 b/src/helio/helio_setup.f90 new file mode 100644 index 000000000..2a8c24019 --- /dev/null +++ b/src/helio/helio_setup.f90 @@ -0,0 +1,22 @@ +submodule(helio_classes) s_helio_setup + use swiftest +contains + + module subroutine helio_setup_initialize_system(self, param) + !! author: David A. Minton + !! + !! Initialize a Helio nbody system from files, converting all heliocentric quantities to barycentric. + !! + implicit none + ! Arguments + class(helio_nbody_system), intent(inout) :: self !! Helio nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + + call whm_setup_initialize_system(self, param) + call self%pl%h2b(self%cb) + call self%tp%h2b(self%cb) + + return + end subroutine helio_setup_initialize_system + +end submodule s_helio_setup \ No newline at end of file diff --git a/src/io/io.f90 b/src/io/io.f90 index 3d1521680..c31718edb 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -40,7 +40,8 @@ module subroutine io_conservation_report(self, param, lterminal) open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "old", action = "write", position = "append", err = 667, iomsg = errmsg) end if end if - call pl%h2b(cb) + call pl%vb2vh(cb) + call pl%xh2xb(cb) call system%get_energy_and_momentum(param) ke_orbit_now = system%ke_orbit ke_spin_now = system%ke_spin @@ -1234,13 +1235,18 @@ module subroutine io_write_discard(self, param) character(*), parameter :: NPLFMT = '(I8)' character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))' class(swiftest_body), allocatable :: pltemp - character(len=STRMAX) :: errmsg + character(len=STRMAX) :: errmsg, out_stat if (param%discard_out == "") return associate(tp_discards => self%tp_discards, nsp => self%tp_discards%nbody, pl => self%pl, npl => self%pl%nbody) if (nsp == 0) return - select case(param%out_stat) + if (lfirst) then + out_stat = param%out_stat + else + out_stat = 'APPEND' + end if + select case(out_stat) case('APPEND') open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) case('NEW', 'REPLACE', 'UNKNOWN') diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index ba3226eb7..c41bf7649 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -14,7 +14,7 @@ module pure subroutine kick_getacch_int_pl(self) class(swiftest_pl), intent(inout) :: self ! Internals integer(I4B) :: k - real(DP) :: rji2, irij3, faci, facj + real(DP) :: rji2, irij3, faci, facj, rlim2 real(DP) :: dx, dy, dz associate(pl => self, npl => self%nbody, nplpl => self%nplpl) @@ -25,15 +25,18 @@ module pure subroutine kick_getacch_int_pl(self) dy = pl%xh(2, j) - pl%xh(2, i) dz = pl%xh(3, j) - pl%xh(3, i) rji2 = dx**2 + dy**2 + dz**2 - irij3 = 1.0_DP / (rji2 * sqrt(rji2)) - faci = pl%Gmass(i) * irij3 - facj = pl%Gmass(j) * irij3 - pl%ah(1, i) = pl%ah(1, i) + facj * dx - pl%ah(2, i) = pl%ah(2, i) + facj * dy - pl%ah(3, i) = pl%ah(3, i) + facj * dz - pl%ah(1, j) = pl%ah(1, j) - faci * dx - pl%ah(2, j) = pl%ah(2, j) - faci * dy - pl%ah(3, j) = pl%ah(3, j) - faci * dz + rlim2 = (pl%radius(i) + pl%radius(j))**2 + if (rji2 > rlim2) then + irij3 = 1.0_DP / (rji2 * sqrt(rji2)) + faci = pl%Gmass(i) * irij3 + facj = pl%Gmass(j) * irij3 + pl%ah(1, i) = pl%ah(1, i) + facj * dx + pl%ah(2, i) = pl%ah(2, i) + facj * dy + pl%ah(3, i) = pl%ah(3, i) + facj * dz + pl%ah(1, j) = pl%ah(1, j) - faci * dx + pl%ah(2, j) = pl%ah(2, j) - faci * dy + pl%ah(3, j) = pl%ah(3, j) - faci * dz + end if end if end associate end do diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index 89f4aa055..efa326a0f 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -15,7 +15,8 @@ module helio_classes !******************************************************************************************************************************** type, extends(whm_nbody_system) :: helio_nbody_system contains - procedure :: step => helio_step_system !! Advance the Helio nbody system forward in time by one step + procedure :: step => helio_step_system !! Advance the Helio nbody system forward in time by one step + procedure :: initialize => helio_setup_initialize_system !! Performs Helio-specific initilization steps, including converting to DH coordinates end type helio_nbody_system !******************************************************************************************************************************** @@ -35,8 +36,6 @@ module helio_classes !! Helio massive body particle class type, extends(swiftest_pl) :: helio_pl contains - procedure :: vh2vb => helio_coord_vh2vb_pl !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only) - procedure :: vb2vh => helio_coord_vb2vh_pl !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only) procedure :: drift => helio_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates procedure :: lindrift => helio_drift_linear_pl !! Method for linear drift of massive bodies due to barycentric momentum of Sun procedure :: accel_gr => helio_gr_kick_getacch_pl !! Acceleration term arising from the post-Newtonian correction @@ -53,8 +52,6 @@ module helio_classes !! Helio test particle class type, extends(swiftest_tp) :: helio_tp contains - procedure :: vh2vb => helio_coord_vh2vb_tp !! Convert test particles from heliocentric to barycentric coordinates (velocity only) - procedure :: vb2vh => helio_coord_vb2vh_tp !! Convert test particles from barycentric to heliocentric coordinates (velocity only) procedure :: lindrift => helio_drift_linear_tp !! Method for linear drift of massive bodies due to barycentric momentum of Sun procedure :: drift => helio_drift_tp !! Method for Danby drift in Democratic Heliocentric coordinates procedure :: accel_gr => helio_gr_kick_getacch_tp !! Acceleration term arising from the post-Newtonian correction @@ -65,32 +62,6 @@ module helio_classes end type helio_tp interface - module subroutine helio_coord_vb2vh_pl(self, cb) - use swiftest_classes, only : swiftest_cb - implicit none - class(helio_pl), intent(inout) :: self !! Helio massive body object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - end subroutine helio_coord_vb2vh_pl - - module subroutine helio_coord_vb2vh_tp(self, vbcb) - implicit none - class(helio_tp), intent(inout) :: self !! Helio massive body object - real(DP), dimension(:), intent(in) :: vbcb !! Barycentric velocity of the central body - end subroutine helio_coord_vb2vh_tp - - module subroutine helio_coord_vh2vb_pl(self, cb) - use swiftest_classes, only : swiftest_cb - implicit none - class(helio_pl), intent(inout) :: self !! Helio massive body object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - end subroutine helio_coord_vh2vb_pl - - module subroutine helio_coord_vh2vb_tp(self, vbcb) - implicit none - class(helio_tp), intent(inout) :: self !! Helio massive body object - real(DP), dimension(:), intent(in) :: vbcb !! Barycentric velocity of the central body - end subroutine helio_coord_vh2vb_tp - module subroutine helio_drift_body(self, system, param, dt) use swiftest_classes, only : swiftest_body, swiftest_nbody_system, swiftest_parameters implicit none @@ -206,6 +177,13 @@ module subroutine helio_kick_vb_tp(self, system, param, t, dt, lbeg) logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. end subroutine helio_kick_vb_tp + module subroutine helio_setup_initialize_system(self, param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(helio_nbody_system), intent(inout) :: self !! Helio nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine helio_setup_initialize_system + module subroutine helio_step_pl(self, system, param, t, dt) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 4cf180e84..aed6f23b0 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -147,8 +147,8 @@ module swiftest_classes logical, dimension(:), allocatable :: ldiscard !! Body should be discarded logical, dimension(:), allocatable :: lmask !! Logical mask used to select a subset of bodies when performing certain operations (drift, kick, accel, etc.) real(DP), dimension(:), allocatable :: mu !! G * (Mcb + [m]) - real(DP), dimension(:,:), allocatable :: xh !! Heliocentric position - real(DP), dimension(:,:), allocatable :: vh !! Heliocentric velocity + real(DP), dimension(:,:), allocatable :: xh !! Swiftestcentric position + real(DP), dimension(:,:), allocatable :: vh !! Swiftestcentric velocity real(DP), dimension(:,:), allocatable :: xb !! Barycentric position real(DP), dimension(:,:), allocatable :: vb !! Barycentric velocity real(DP), dimension(:,:), allocatable :: ah !! Total heliocentric acceleration @@ -226,7 +226,9 @@ module swiftest_classes procedure :: append => util_append_pl !! Appends elements from one structure to another procedure :: h2b => util_coord_h2b_pl !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity) procedure :: b2h => util_coord_b2h_pl !! Convert massive bodies from barycentric to heliocentric coordinates (position and velocity) - procedure :: xh2xb => util_coord_xh2xb_pl !! Convert position vectors of massive bodies from heliocentric to barycentric coordinates (position and velocity) + procedure :: vh2vb => util_coord_vh2vb_pl !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only) + procedure :: vb2vh => util_coord_vb2vh_pl !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only) + procedure :: xh2xb => util_coord_xh2xb_pl !! Convert massive bodies from heliocentric to barycentric coordinates (position only) procedure :: fill => util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) procedure :: resize => util_resize_pl !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small. procedure :: set_beg_end => util_set_beg_end_pl !! Sets the beginning and ending positions and velocities of planets. @@ -258,6 +260,9 @@ module swiftest_classes procedure :: append => util_append_tp !! Appends elements from one structure to another procedure :: h2b => util_coord_h2b_tp !! Convert test particles from heliocentric to barycentric coordinates (position and velocity) procedure :: b2h => util_coord_b2h_tp !! Convert test particles from barycentric to heliocentric coordinates (position and velocity) + procedure :: vb2vh => util_coord_vb2vh_tp !! Convert test particles from barycentric to heliocentric coordinates (velocity only) + procedure :: vh2vb => util_coord_vh2vb_tp !! Convert test particles from heliocentric to barycentric coordinates (velocity only) + procedure :: xh2xb => util_coord_xh2xb_tp !! Convert test particles from heliocentric to barycentric coordinates (position only) procedure :: fill => util_fill_tp !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) procedure :: get_peri => util_peri_tp !! Determine system pericenter passages for test particles procedure :: resize => util_resize_tp !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small. @@ -516,9 +521,9 @@ module pure subroutine gr_pseudovel2vel(param, mu, xh, pv, vh) implicit none class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body - real(DP), dimension(:), intent(in) :: xh !! Heliocentric position vector + real(DP), dimension(:), intent(in) :: xh !! Swiftestcentric position vector real(DP), dimension(:), intent(in) :: pv !! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32) - real(DP), dimension(:), intent(out) :: vh !! Heliocentric velocity vector + real(DP), dimension(:), intent(out) :: vh !! Swiftestcentric velocity vector end subroutine gr_pseudovel2vel module pure subroutine gr_pv2vh_body(self, param) @@ -531,8 +536,8 @@ module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv) implicit none class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body - real(DP), dimension(:), intent(in) :: xh !! Heliocentric position vector - real(DP), dimension(:), intent(in) :: vh !! Heliocentric velocity vector + real(DP), dimension(:), intent(in) :: xh !! Swiftestcentric position vector + real(DP), dimension(:), intent(in) :: vh !! Swiftestcentric velocity vector real(DP), dimension(:), intent(out) :: pv !! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32) end subroutine gr_vel2pseudovel @@ -687,8 +692,8 @@ module subroutine io_write_frame_encounter(iu, t, id1, id2, Gmass1, Gmass2, radi integer(I4B), intent(in) :: id1, id2 !! ids of the two encountering bodies real(DP), intent(in) :: Gmass1, Gmass2 !! G*mass of the two encountering bodies real(DP), intent(in) :: radius1, radius2 !! Radii of the two encountering bodies - real(DP), dimension(:), intent(in) :: xh1, xh2 !! Heliocentric position vectors of the two encountering bodies - real(DP), dimension(:), intent(in) :: vh1, vh2 !! Heliocentric velocity vectors of the two encountering bodies + real(DP), dimension(:), intent(in) :: xh1, xh2 !! Swiftestcentric position vectors of the two encountering bodies + real(DP), dimension(:), intent(in) :: vh1, vh2 !! Swiftestcentric velocity vectors of the two encountering bodies end subroutine io_write_frame_encounter module subroutine io_write_frame_system(self, iu, param) @@ -937,12 +942,42 @@ module subroutine util_coord_h2b_tp(self, cb) class(swiftest_cb), intent(in) :: cb !! Swiftest central body object end subroutine util_coord_h2b_tp + module subroutine util_coord_vb2vh_pl(self, cb) + implicit none + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + end subroutine util_coord_vb2vh_pl + + module subroutine util_coord_vb2vh_tp(self, vbcb) + implicit none + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + real(DP), dimension(:), intent(in) :: vbcb !! Barycentric velocity of the central body + end subroutine util_coord_vb2vh_tp + + module subroutine util_coord_vh2vb_pl(self, cb) + implicit none + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + end subroutine util_coord_vh2vb_pl + + module subroutine util_coord_vh2vb_tp(self, vbcb) + implicit none + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + real(DP), dimension(:), intent(in) :: vbcb !! Barycentric velocity of the central body + end subroutine util_coord_vh2vb_tp + module subroutine util_coord_xh2xb_pl(self, cb) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object end subroutine util_coord_xh2xb_pl + module subroutine util_coord_xh2xb_tp(self, cb) + implicit none + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_cb), intent(in) :: cb !! Swiftest central body object + end subroutine util_coord_xh2xb_tp + module subroutine util_copy_encounter(self, source) implicit none class(swiftest_encounter), intent(inout) :: self !! Encounter list diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index cbe626e43..454f454f5 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -72,23 +72,24 @@ module swiftest_globals integer(I4B), parameter :: HYPERBOLA = 1 !! Symbolic names for orbit types - hyperbola !> Symbolic names for body/particle status codes: - integer(I4B), parameter :: ACTIVE = 0 - integer(I4B), parameter :: INACTIVE = 1 - integer(I4B), parameter :: DISCARDED_RMAX = -1 - integer(I4B), parameter :: DISCARDED_RMIN = -2 - integer(I4B), parameter :: DISCARDED_RMAXU = -3 - integer(I4B), parameter :: DISCARDED_PERI = -4 - integer(I4B), parameter :: DISCARDED_PLR = -5 - integer(I4B), parameter :: DISCARDED_PLQ = -6 - integer(I4B), parameter :: DISCARDED_DRIFTERR = -7 - integer(I4B), parameter :: MERGED = -8 - integer(I4B), parameter :: DISRUPTION = -9 - integer(I4B), parameter :: SUPERCATASTROPHIC = -10 - integer(I4B), parameter :: GRAZE_AND_MERGE = -11 - integer(I4B), parameter :: HIT_AND_RUN = -12 - integer(I4B), parameter :: COLLISION = -13 - integer(I4B), parameter :: NEW_PARTICLE = -14 - integer(I4B), parameter :: OLD_PARTICLE = -15 + integer(I4B), parameter :: ACTIVE = 0 + integer(I4B), parameter :: INACTIVE = 1 + integer(I4B), parameter :: DISCARDED_RMAX = -1 + integer(I4B), parameter :: DISCARDED_RMIN = -2 + integer(I4B), parameter :: DISCARDED_RMAXU = -3 + integer(I4B), parameter :: DISCARDED_PERI = -4 + integer(I4B), parameter :: DISCARDED_PLR = -5 + integer(I4B), parameter :: DISCARDED_PLQ = -6 + integer(I4B), parameter :: DISCARDED_DRIFTERR = -7 + integer(I4B), parameter :: MERGED = -8 + integer(I4B), parameter :: DISRUPTION = -9 + integer(I4B), parameter :: SUPERCATASTROPHIC = -10 + integer(I4B), parameter :: GRAZE_AND_MERGE = -11 + integer(I4B), parameter :: HIT_AND_RUN_DISRUPT = -12 + integer(I4B), parameter :: HIT_AND_RUN_PURE = -13 + integer(I4B), parameter :: COLLISION = -14 + integer(I4B), parameter :: NEW_PARTICLE = -15 + integer(I4B), parameter :: OLD_PARTICLE = -16 !>Symbolic names for collisional outcomes from collresolve_resolve: integer(I4B), parameter :: COLLRESOLVE_REGIME_MERGE = 1 diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 9ba20d243..05698d44c 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -225,20 +225,24 @@ module subroutine symba_collision_resolve_mergers(self, system, param) class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions end subroutine symba_collision_resolve_mergers - module subroutine symba_collision_resolve_plplenc(self, system, param, t) + module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, irec) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Current simulation step size + integer(I4B), intent(in) :: irec !! Current recursion level end subroutine symba_collision_resolve_plplenc - module subroutine symba_collision_resolve_pltpenc(self, system, param, t) + module subroutine symba_collision_resolve_pltpenc(self, system, param, t, dt, irec) implicit none class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Current simulation step size + integer(I4B), intent(in) :: irec !! Current recursion level end subroutine symba_collision_resolve_pltpenc module subroutine symba_discard_pl(self, system, param) @@ -497,9 +501,10 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci) integer(I4B), value :: ireci !! input recursion level end subroutine symba_step_recur_system - module subroutine symba_step_reset_system(self) + module subroutine symba_step_reset_system(self, param) implicit none - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions end subroutine symba_step_reset_system end interface diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index 6d5b26394..0b9caebcc 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -73,8 +73,8 @@ module whm_classes type, extends(swiftest_nbody_system) :: whm_nbody_system contains !> Replace the abstract procedures with concrete ones - procedure :: initialize => whm_setup_initialize_system !! Performs WHM-specific initilization steps, like calculating the Jacobi masses - procedure :: step => whm_step_system !! Advance the WHM nbody system forward in time by one step + procedure :: initialize => whm_setup_initialize_system ! Performs WHM-specific initilization steps, like calculating the Jacobi masses + procedure :: step => whm_step_system !! Advance the WHM nbody system forward in time by one step end type whm_nbody_system interface diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index f1f42d3cf..f813e068a 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -186,15 +186,15 @@ module function symba_collision_casehitandrun(system, param, family, x, v, mass, end if end if if (lpure) then ! Reset these bodies back to being active so that nothing further is done to them - status = ACTIVE + status = HIT_AND_RUN_PURE select type(pl => system%pl) class is (symba_pl) - pl%status(family(:)) = status + pl%status(family(:)) = ACTIVE pl%ldiscard(family(:)) = .false. pl%lcollision(family(:)) = .false. end select else - status = HIT_AND_RUN + status = HIT_AND_RUN_DISRUPT call symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) end if @@ -871,7 +871,7 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, plnew%info(i)%origin_xh(:) = plnew%xh(:,i) plnew%info(i)%origin_vh(:) = plnew%vh(:,i) end do - case(HIT_AND_RUN) + case(HIT_AND_RUN_DISRUPT) plnew%info(1) = pl%info(ibiggest) plnew%status(1) = OLD_PARTICLE plnew%status(2:nfrag) = NEW_PARTICLE @@ -938,7 +938,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) 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, jtarg, jproj, regime + integer(I4B) :: i, jtarg, jproj, regime real(DP), dimension(2) :: radius_si, mass_si, density_si real(DP) :: mtiny_si, Mcb_si real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si @@ -946,7 +946,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) integer(I4B), parameter :: NRES = 3 !! Number of collisional product results real(DP), dimension(NRES) :: mass_res - associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) + associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) select type(pl => system%pl) class is (symba_pl) select type (cb => system%cb) @@ -969,10 +969,10 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) end if mass_si(:) = (mass(:)) * param%MU2KG !! The collective mass of the parent and its children radius_si(:) = radius(:) * param%DU2M !! The collective radius of the parent and its children - x1_si(:) = plpl_collisions%x1(:,i) * param%DU2M !! The position of the parent from inside the step (at collision) - v1_si(:) = plpl_collisions%v1(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) - x2_si(:) = plpl_collisions%x2(:,i) * param%DU2M !! The position of the parent from inside the step (at collision) - v2_si(:) = plpl_collisions%v2(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) + x1_si(:) = plplcollision_list%x1(:,i) * param%DU2M !! The position of the parent from inside the step (at collision) + v1_si(:) = plplcollision_list%v1(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) + x2_si(:) = plplcollision_list%x2(:,i) * param%DU2M !! The position of the parent from inside the step (at collision) + v2_si(:) = plplcollision_list%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%GMTINY / param%GU) * param%MU2KG @@ -994,13 +994,13 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) select case (regime) case (COLLRESOLVE_REGIME_DISRUPTION) - status = symba_collision_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + plplcollision_list%status(i) = symba_collision_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) case (COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - status = symba_collision_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + plplcollision_list%status(i) = symba_collision_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) case (COLLRESOLVE_REGIME_HIT_AND_RUN) - status = symba_collision_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + plplcollision_list%status(i) = symba_collision_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - status = symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) + plplcollision_list%status(i) = symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) case default write(*,*) "Error in symba_collision, unrecognized collision regime" call util_exit(FAILURE) @@ -1030,9 +1030,9 @@ module subroutine symba_collision_resolve_mergers(self, system, param) 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 + integer(I4B) :: i - associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) + associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) select type(pl => system%pl) class is (symba_pl) select type(cb => system%cb) @@ -1044,7 +1044,7 @@ module subroutine symba_collision_resolve_mergers(self, system, param) if (.not. lgoodcollision) cycle if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved - status = symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) + plplcollision_list%status(i) = symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) end do end select end select @@ -1054,7 +1054,7 @@ module subroutine symba_collision_resolve_mergers(self, system, param) end subroutine symba_collision_resolve_mergers - module subroutine symba_collision_resolve_plplenc(self, system, param, t) + module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, irec) !! author: David A. Minton !! !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the collision @@ -1065,8 +1065,11 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t) class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Current simulation step size + integer(I4B), intent(in) :: irec !! Current recursion level ! Internals real(DP) :: Eorbit_before, Eorbit_after + logical :: lplpl_collision associate(plplenc_list => self, plplcollision_list => system%plplcollision_list) select type(pl => system%pl) @@ -1074,26 +1077,44 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t) select type(param) class is (symba_parameters) if (plplcollision_list%nenc == 0) return ! No collisions to resolve - - write(*, *) "Collision between massive bodies detected at time t = ", t - call pl%vb2vh(system%cb) + ! Make sure that the heliocentric and barycentric coordinates are consistent with each other + call pl%vb2vh(system%cb) call pl%xh2xb(system%cb) - if (param%lfragmentation) then - call plplcollision_list%resolve_fragmentations(system, param) - else - call plplcollision_list%resolve_mergers(system, param) - end if - - ! Destroy the collision list now that the collisions are resolved - call plplcollision_list%setup(0) - + ! Get the energy before the collision is resolved if (param%lenergy) then call system%get_energy_and_momentum(param) Eorbit_before = system%te end if - call pl%rearray(system, param) + do + write(*, *) "Collision between massive bodies detected at time t = ", t + if (param%lfragmentation) then + call plplcollision_list%resolve_fragmentations(system, param) + else + call plplcollision_list%resolve_mergers(system, param) + end if + + ! Destroy the collision list now that the collisions are resolved + call plplcollision_list%setup(0) + + if ((system%pl_adds%nbody == 0) .and. (system%pl_discards%nbody == 0)) exit + + ! Save the add/discard information to file + call system%write_discard(param) + + ! Rearrange the arrays: Remove discarded bodies, add any new bodies, resort, and recompute all indices and encounter lists + call pl%rearray(system, param) + + ! Destroy the add/discard list so that we don't append the same body multiple times if another collision is detected + call system%pl_discards%setup(0, param) + call system%pl_adds%setup(0, param) + + ! Check whether or not any of the particles that were just added are themselves in a collision state. This will generate a new plplcollision_list + lplpl_collision = plplenc_list%collision_check(system, param, t, dt, irec) + + if (.not.lplpl_collision) exit + end do if (param%lenergy) then call system%get_energy_and_momentum(param) @@ -1109,7 +1130,7 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t) end subroutine symba_collision_resolve_plplenc - module subroutine symba_collision_resolve_pltpenc(self, system, param, t) + module subroutine symba_collision_resolve_pltpenc(self, system, param, t, dt, irec) !! author: David A. Minton !! !! Process the pl-tp collision list, then modifiy the massive bodies based on the outcome of the collision @@ -1120,12 +1141,13 @@ module subroutine symba_collision_resolve_pltpenc(self, system, param, t) class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions real(DP), intent(in) :: t !! Current simulation tim - + real(DP), intent(in) :: dt !! Current simulation step size + integer(I4B), intent(in) :: irec !! Current recursion level + + call system%tp%xh2xb(system%cb) call system%tp%discard(system, param) return end subroutine symba_collision_resolve_pltpenc - - end submodule s_symba_collision \ No newline at end of file diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index f6f7c0e6b..bf70f15dd 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -176,10 +176,6 @@ subroutine symba_discard_nonplpl(pl, system, param) ! First check for collisions with the central body associate(npl => pl%nbody, cb => system%cb) if (npl == 0) return - if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. & - (param%rmaxu >= 0.0_DP) .or. ((param%qmin >= 0.0_DP) .and. (param%qmin_coord == "BARY"))) then - call pl%h2b(cb) - end if if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) then call symba_discard_cb_pl(pl, system, param) end if @@ -290,6 +286,8 @@ module subroutine symba_discard_pl(self, system, param) select type(param) class is (symba_parameters) associate(pl => self, plplenc_list => system%plplenc_list, plplcollision_list => system%plplcollision_list) + call pl%vb2vh(system%cb) + call pl%xh2xb(system%cb) call plplenc_list%write(pl, pl, param) call symba_discard_nonplpl(self, system, param) diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 4bdc5c195..775efa962 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -329,7 +329,7 @@ module subroutine symba_io_write_discard(self, param) character(*), parameter :: NPLFMT = '(I8)' character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))' class(swiftest_body), allocatable :: pltemp - character(STRMAX) :: errmsg + character(STRMAX) :: errmsg, out_stat if (param%discard_out == "") return @@ -338,7 +338,12 @@ module subroutine symba_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) + if (lfirst) then + out_stat = param%out_stat + else + out_stat = 'APPEND' + end if + select case(out_stat) case('APPEND') open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) case('NEW', 'REPLACE', 'UNKNOWN') diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index 307c37675..67f8063b5 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -73,7 +73,7 @@ module subroutine symba_setup_initialize_system(self, param) ! Call parent method associate(system => self) - call whm_setup_initialize_system(system, param) + call helio_setup_initialize_system(system, param) call system%pltpenc_list%setup(0) call system%plplenc_list%setup(0) call system%plplcollision_list%setup(0) diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index ca5aa43cf..249c73a82 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -19,20 +19,23 @@ module subroutine symba_step_system(self, param, t, dt) ! Internals logical :: lencounter - call self%reset() select type(pl => self%pl) class is (symba_pl) select type(tp => self%tp) class is (symba_tp) - lencounter = pl%encounter_check(self, dt, 0) .or. tp%encounter_check(self, dt, 0) - if (lencounter) then - tp%lfirst = pl%lfirst - call self%interp(param, t, dt) - pl%lfirst = .true. - tp%lfirst = .true. - else - call helio_step_system(self, param, t, dt) - end if + select type(param) + class is (symba_parameters) + call self%reset(param) + lencounter = pl%encounter_check(self, dt, 0) .or. tp%encounter_check(self, dt, 0) + if (lencounter) then + tp%lfirst = pl%lfirst + call self%interp(param, t, dt) + pl%lfirst = .true. + tp%lfirst = .true. + else + call helio_step_system(self, param, t, dt) + end if + end select end select end select @@ -197,8 +200,8 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci) lplpl_collision = plplenc_list%collision_check(system, param, t+dtl, dtl, ireci) lpltp_collision = pltpenc_list%collision_check(system, param, t+dtl, dtl, ireci) - if (lplpl_collision) call plplenc_list%resolve_collision(system, param, t+dtl) - if (lpltp_collision) call pltpenc_list%resolve_collision(system, param, t+dtl) + if (lplpl_collision) call plplenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) + if (lpltp_collision) call pltpenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) end if call self%set_recur_levels(ireci) @@ -212,7 +215,7 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci) end subroutine symba_step_recur_system - module subroutine symba_step_reset_system(self) + module subroutine symba_step_reset_system(self, param) !! author: David A. Minton !! !! Resets pl, tp,and encounter structures at the start of a new step @@ -221,7 +224,8 @@ module subroutine symba_step_reset_system(self) !! Adapted from Hal Levison's Swift routine symba5_step.f implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions ! Internals integer(I4B) :: i @@ -254,8 +258,8 @@ module subroutine symba_step_reset_system(self) system%pltpenc_list%nenc = 0 end if - call system%pl_adds%resize(0) - call system%pl_discards%resize(0) + call system%pl_adds%setup(0, param) + call system%pl_discards%setup(0, param) end select end select end associate diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 12919bed0..a80c8e646 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -424,25 +424,25 @@ module subroutine symba_util_rearray_pl(self, system, param) associate(pl => self, pl_adds => system%pl_adds) - allocate(tmp, mold=self) + ! Deallocate any temporary variables + if (allocated(pl%xbeg)) deallocate(pl%xbeg) + if (allocated(pl%xend)) deallocate(pl%xend) + ! Remove the discards and destroy the list, as the system already tracks pl_discards elsewhere allocate(lmask, source=pl%ldiscard(:)) lmask(:) = lmask(:) .or. pl%status(:) == INACTIVE + allocate(tmp, mold=self) call pl%spill(tmp, lspill_list=lmask, ldestructive=.true.) - deallocate(lmask) call tmp%setup(0,param) deallocate(tmp) + deallocate(lmask) - ! Deallocate any temporary variables - if (allocated(pl%xbeg)) deallocate(pl%xbeg) - if (allocated(pl%xend)) deallocate(pl%xend) + ! Store the original plplenc list so we don't remove any of the original encounters + allocate(plplenc_old, source=system%plplenc_list) + call plplenc_old%copy(system%plplenc_list) ! Add in any new bodies if (pl_adds%nbody > 0) then - ! First store the original plplenc list so we don't remove any of the original encounters - allocate(plplenc_old, source=system%plplenc_list) - call plplenc_old%copy(system%plplenc_list) - ! Append the adds to the main pl object call pl%append(pl_adds, lsource_mask=[(.true., i=1, pl_adds%nbody)]) @@ -450,47 +450,65 @@ module subroutine symba_util_rearray_pl(self, system, param) lmask(:) = pl%status(1:pl%nbody) == NEW_PARTICLE call symba_io_dump_particle_info(system, param, plidx=pack([(i, i=1, pl%nbody)], lmask)) - where(pl%status(:) /= INACTIVE) pl%status(:) = ACTIVE - call pl%sort("mass", ascending=.false.) - pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY - pl%nplm = count(pl%lmtiny(:)) - - ! Reindex the bodies and calculate the level 0 encounter list - call pl%eucl_index() - lencounter = pl%encounter_check(system, param%dt, 0) - select type(tp => system%tp) - class is (symba_tp) - lencounter = tp%encounter_check(system, param%dt, 0) - end select - - associate(idnew1 => system%plplenc_list%id1, idnew2 => system%plplenc_list%id2, idold1 => plplenc_old%id1, idold2 => plplenc_old%id2) - do k = 1, system%plplenc_list%nenc - if ((idnew1(k) == idold1(k)) .and. (idnew2(k) == idold2(k))) then - ! This is an encounter we already know about, so save the old information - system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) - system%plplenc_list%status(k) = plplenc_old%status(k) - system%plplenc_list%x1(:,k) = plplenc_old%x1(:,k) - system%plplenc_list%x2(:,k) = plplenc_old%x2(:,k) - system%plplenc_list%v1(:,k) = plplenc_old%v1(:,k) - system%plplenc_list%v2(:,k) = plplenc_old%v2(:,k) - system%plplenc_list%t(k) = plplenc_old%t(k) - system%plplenc_list%level(k) = plplenc_old%level(k) - else if (((idnew1(k) == idold2(k)) .and. (idnew2(k) == idold1(k)))) then - ! This is an encounter we already know about, but with the order reversed, so save the old information - system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) - system%plplenc_list%status(k) = plplenc_old%status(k) - system%plplenc_list%x1(:,k) = plplenc_old%x2(:,k) - system%plplenc_list%x2(:,k) = plplenc_old%x1(:,k) - system%plplenc_list%v1(:,k) = plplenc_old%v2(:,k) - system%plplenc_list%v2(:,k) = plplenc_old%v1(:,k) - system%plplenc_list%t(k) = plplenc_old%t(k) - system%plplenc_list%level(k) = plplenc_old%level(k) - end if - end do - end associate + deallocate(lmask) end if + ! Reset all of the status flags for this body + where(pl%status(:) /= INACTIVE) + pl%status(:) = ACTIVE + pl%ldiscard(:) = .false. + pl%lcollision(:) = .false. + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY + pl%lmask(:) = .true. + elsewhere + pl%ldiscard(:) = .true. + pl%lmask(:) = .false. + end where + + pl%nplm = count(pl%lmtiny(1:pl%nbody) .and. pl%lmask(1:pl%nbody)) + + ! Reindex the new list of bodies + call pl%sort("mass", ascending=.false.) + call pl%eucl_index() + + ! Reset the kinship trackers + pl%kin(1:pl%nbody)%nchild = 0 + pl%kin(1:pl%nbody)%parent = [(i, i=1, pl%nbody)] + + ! Re-build the encounter list + lencounter = pl%encounter_check(system, param%dt, 0) + select type(tp => system%tp) + class is (symba_tp) + lencounter = tp%encounter_check(system, param%dt, 0) + end select + + associate(idnew1 => system%plplenc_list%id1, idnew2 => system%plplenc_list%id2, idold1 => plplenc_old%id1, idold2 => plplenc_old%id2) + do k = 1, system%plplenc_list%nenc + if ((idnew1(k) == idold1(k)) .and. (idnew2(k) == idold2(k))) then + ! This is an encounter we already know about, so save the old information + system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) + system%plplenc_list%status(k) = plplenc_old%status(k) + system%plplenc_list%x1(:,k) = plplenc_old%x1(:,k) + system%plplenc_list%x2(:,k) = plplenc_old%x2(:,k) + system%plplenc_list%v1(:,k) = plplenc_old%v1(:,k) + system%plplenc_list%v2(:,k) = plplenc_old%v2(:,k) + system%plplenc_list%t(k) = plplenc_old%t(k) + system%plplenc_list%level(k) = plplenc_old%level(k) + else if (((idnew1(k) == idold2(k)) .and. (idnew2(k) == idold1(k)))) then + ! This is an encounter we already know about, but with the order reversed, so save the old information + system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) + system%plplenc_list%status(k) = plplenc_old%status(k) + system%plplenc_list%x1(:,k) = plplenc_old%x2(:,k) + system%plplenc_list%x2(:,k) = plplenc_old%x1(:,k) + system%plplenc_list%v1(:,k) = plplenc_old%v2(:,k) + system%plplenc_list%v2(:,k) = plplenc_old%v1(:,k) + system%plplenc_list%t(k) = plplenc_old%t(k) + system%plplenc_list%level(k) = plplenc_old%level(k) + end if + end do + end associate + end associate return diff --git a/src/util/util_append.f90 b/src/util/util_append.f90 index 88755870f..971e3c950 100644 --- a/src/util/util_append.f90 +++ b/src/util/util_append.f90 @@ -245,8 +245,6 @@ module subroutine util_append_pl(self, source, lsource_mask) call util_append_body(self, source, lsource_mask) end associate - - call self%eucl_index() class default write(*,*) "Invalid object passed to the append method. Source must be of class swiftest_pl or its descendents" call util_exit(FAILURE) diff --git a/src/util/util_coord.f90 b/src/util/util_coord.f90 index 4daad0e53..b4a77231c 100644 --- a/src/util/util_coord.f90 +++ b/src/util/util_coord.f90 @@ -42,46 +42,10 @@ module subroutine util_coord_h2b_pl(self, cb) end subroutine util_coord_h2b_pl - module subroutine util_coord_xh2xb_pl(self, cb) - !! author: David A. Minton - !! - !! Convert position vectors of massive bodies from heliocentric to barycentric coordinates (position and velocity) - !! - !! Adapted from David E. Kaufmann's Swifter routine coord_h2b.f90 - !! Adapted from Hal Levison's Swift routine coord_h2b.f - implicit none - ! Arguments - class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - ! Internals - integer(I4B) :: i - real(DP) :: Gmtot - real(DP), dimension(NDIM) :: xtmp - - if (self%nbody == 0) return - associate(pl => self, npl => self%nbody) - Gmtot = cb%Gmass - xtmp(:) = 0.0_DP - do i = 1, npl - if (pl%status(i) == INACTIVE) cycle - Gmtot = Gmtot + pl%Gmass(i) - xtmp(:) = xtmp(:) + pl%Gmass(i) * pl%xh(:,i) - end do - cb%xb(:) = -xtmp(:) / Gmtot - do i = 1, npl - if (pl%status(i) == INACTIVE) cycle - pl%xb(:,i) = pl%xh(:,i) + cb%xb(:) - end do - end associate - - return - end subroutine util_coord_xh2xb_pl - - module subroutine util_coord_h2b_tp(self, cb) !! author: David A. Minton !! - !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity) + !! Convert test particles from heliocentric to barycentric coordinates (position and velocity) !! !! Adapted from David E. Kaufmann's Swifter routine coord_h2b_tp.f90 !! Adapted from Hal Levison's Swift routine coord_h2b_tp.f @@ -156,5 +120,174 @@ module subroutine util_coord_b2h_tp(self, cb) return end subroutine util_coord_b2h_tp + + + module subroutine util_coord_vb2vh_pl(self, cb) + !! author: David A. Minton + !! + !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only) + !! + !! Adapted from David E. Kaufmann's Swifter routine coord_vb2vh.f90 + !! Adapted from Hal Levison's Swift routine coord_vb2vh.f + implicit none + ! Arguments + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + ! Internals + integer(I4B) :: i + + if (self%nbody == 0) return + + associate(pl => self, npl => self%nbody) + do i = 1, NDIM + cb%vb(i) = -sum(pl%Gmass(1:npl) * pl%vb(i, 1:npl)) / cb%Gmass + pl%vh(i, 1:npl) = pl%vb(i, 1:npl) - cb%vb(i) + end do + end associate + + return + end subroutine util_coord_vb2vh_pl + + + module subroutine util_coord_vb2vh_tp(self, vbcb) + !! author: David A. Minton + !! + !! Convert test particles from barycentric to heliocentric coordinates (velocity only) + !! + !! Adapted from David E. Kaufmann's Swifter routine coord_vb2vh_tp.f90 + !! Adapted from Hal Levison's Swift routine coord_vb2h_tp.f + implicit none + ! Arguments + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + real(DP), dimension(:), intent(in) :: vbcb !! Barycentric velocity of the central body + + if (self%nbody == 0) return + + associate(tp => self, ntp => self%nbody) + where (tp%lmask(1:ntp)) + tp%vh(1, 1:ntp) = tp%vb(1, 1:ntp) - vbcb(1) + tp%vh(2, 1:ntp) = tp%vb(2, 1:ntp) - vbcb(2) + tp%vh(3, 1:ntp) = tp%vb(3, 1:ntp) - vbcb(3) + end where + end associate + + return + end subroutine util_coord_vb2vh_tp + + + module subroutine util_coord_vh2vb_pl(self, cb) + !! author: David A. Minton + !! + !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only) + !! + !! Adapted from David E. Kaufmann's Swifter routine coord_vh2vb.f90 + !! Adapted from Hal Levison's Swift routine coord_vh2b.f + implicit none + ! Arguments + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + ! Internals + integer(I4B) :: i + real(DP) :: Gmtot + + if (self%nbody == 0) return + + associate(pl => self, npl => self%nbody) + Gmtot = cb%Gmass + sum(pl%Gmass(1:npl)) + do i = 1, NDIM + cb%vb(i) = -sum(pl%Gmass(1:npl) * pl%vh(i, 1:npl)) / Gmtot + pl%vb(i, 1:npl) = pl%vh(i, 1:npl) + cb%vb(i) + end do + end associate + + return + end subroutine util_coord_vh2vb_pl + + + module subroutine util_coord_vh2vb_tp(self, vbcb) + !! author: David A. Minton + !! + !! Convert test particles from heliocentric to barycentric coordinates (velocity only) + !! + !! Adapted from David E. Kaufmann's Swifter routine coord_vh2vb_tp.f90 + !! Adapted from Hal Levison's Swift routine coord_vh2b_tp.f + implicit none + ! Arguments + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + real(DP), dimension(:), intent(in) :: vbcb !! Barycentric velocity of the central body + + if (self%nbody == 0) return + + associate(tp => self, ntp => self%nbody) + where (tp%lmask(1:ntp)) + tp%vb(1, 1:ntp) = tp%vh(1, 1:ntp) + vbcb(1) + tp%vb(2, 1:ntp) = tp%vh(2, 1:ntp) + vbcb(2) + tp%vb(3, 1:ntp) = tp%vh(3, 1:ntp) + vbcb(3) + end where + end associate + + return + end subroutine util_coord_vh2vb_tp + + module subroutine util_coord_xh2xb_pl(self, cb) + !! author: David A. Minton + !! + !! Convert position vectors of massive bodies from heliocentric to barycentric coordinates (position only) + !! + !! Adapted from David E. Kaufmann's Swifter routine coord_h2b.f90 + !! Adapted from Hal Levison's Swift routine coord_h2b.f + implicit none + ! Arguments + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + ! Internals + integer(I4B) :: i + real(DP) :: Gmtot + real(DP), dimension(NDIM) :: xtmp + + if (self%nbody == 0) return + associate(pl => self, npl => self%nbody) + Gmtot = cb%Gmass + xtmp(:) = 0.0_DP + do i = 1, npl + if (pl%status(i) == INACTIVE) cycle + Gmtot = Gmtot + pl%Gmass(i) + xtmp(:) = xtmp(:) + pl%Gmass(i) * pl%xh(:,i) + end do + cb%xb(:) = -xtmp(:) / Gmtot + do i = 1, npl + if (pl%status(i) == INACTIVE) cycle + pl%xb(:,i) = pl%xh(:,i) + cb%xb(:) + end do + end associate + + return + end subroutine util_coord_xh2xb_pl + + + module subroutine util_coord_xh2xb_tp(self, cb) + !! author: David A. Minton + !! + !! Convert test particles from heliocentric to barycentric coordinates (position only) + !! + !! Adapted from David E. Kaufmann's Swifter routine coord_h2b_tp.f90 + !! Adapted from Hal Levison's Swift routine coord_h2b_tp.f + implicit none + ! Arguments + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_cb), intent(in) :: cb !! Swiftest central body object + ! Internals + integer(I4B) :: i + + if (self%nbody == 0) return + associate(tp => self, ntp => self%nbody) + do concurrent (i = 1:ntp, tp%status(i) /= INACTIVE) + tp%xb(:, i) = tp%xh(:, i) + cb%xb(:) + end do + end associate + + return + end subroutine util_coord_xh2xb_tp + end submodule s_util_coord \ No newline at end of file diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 7535ea67b..9a4060dc1 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -93,7 +93,7 @@ module subroutine util_get_energy_momentum_system(self, param) end do do concurrent (k = 1:pl%nplpl, lstatpl(k)) - pepl(k) = -pl%Gmass(indi(k)) * pl%mass(indj(k)) / norm2(pl%xb(:, indi(k)) - pl%xb(:, indj(k))) + pepl(k) = -(pl%Gmass(indi(k)) * pl%mass(indj(k))) / norm2(pl%xb(:, indi(k)) - pl%xb(:, indj(k))) end do end associate diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 index eaed16c14..ad44fbefb 100644 --- a/src/whm/whm_setup.f90 +++ b/src/whm/whm_setup.f90 @@ -73,7 +73,7 @@ module subroutine whm_setup_initialize_system(self, param) !! implicit none ! Arguments - class(whm_nbody_system), intent(inout) :: self !! Swiftest system object + class(whm_nbody_system), intent(inout) :: self !! WHM nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters call setup_initialize_system(self, param)