diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index fafe319f8..c9603a2a8 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -73,16 +73,16 @@ rot_vectors = {"disruption_headon" : [np.array([0.0, 0.0, 0.0]), np.array([0.0, 0.0, 0.0])], - "disruption_off_axis": [np.array([0.0, 0.0, -6.0e4]), - np.array([0.0, 0.0, 1.0e5])], + "disruption_off_axis": [np.array([0.0, 0.0, -6.0e3]), + np.array([0.0, 0.0, 1.0e4])], "supercatastrophic_headon": [np.array([0.0, 0.0, 0.0]), np.array([0.0, 0.0, 0.0])], - "supercatastrophic_off_axis": [np.array([0.0, 0.0, -6.0e4]), - np.array([0.0, 0.0, 1.0e5])], - "hitandrun_disrupt" : [np.array([0.0, 0.0, 6.0e4]), - np.array([0.0, 0.0, 1.0e5])], - "hitandrun_pure" : [np.array([0.0, 0.0, 6.0e4]), - np.array([0.0, 0.0, 1.0e5])] + "supercatastrophic_off_axis": [np.array([0.0, 0.0, -6.0e3]), + np.array([0.0, 0.0, 1.0e4])], + "hitandrun_disrupt" : [np.array([0.0, 0.0, 6.0e3]), + np.array([0.0, 0.0, 1.0e4])], + "hitandrun_pure" : [np.array([0.0, 0.0, 6.0e3]), + np.array([0.0, 0.0, 1.0e4])] } body_Gmass = {"disruption_headon" : [1e-7, 1e-10], diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index b55fe3d98..e3f818d7e 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -118,7 +118,7 @@ subroutine abstract_io_dump_param(self, param_file_name) import base_parameters implicit none class(base_parameters),intent(in) :: self !! Output collection of parameters - character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) + character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) end subroutine abstract_io_dump_param subroutine abstract_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 9a23ee4f4..cc0363b84 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -45,7 +45,7 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) class(base_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! The time of the collision ! Internals - integer(I4B) :: i,j,nfrag + integer(I4B) :: i,j,nimp real(DP), dimension(NDIM) :: vcom, rnorm select type(nbody_system) @@ -62,8 +62,8 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) allocate(before%pl, source=pl) end select - nfrag = size(impactors%id(:)) - do i = 1, nfrag + nimp = size(impactors%id(:)) + do i = 1, nimp j = impactors%id(i) vcom(:) = pl%vb(:,j) - impactors%vbcom(:) rnorm(:) = .unit. (impactors%rb(:,2) - impactors%rb(:,1)) @@ -183,7 +183,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) associate(fragments => nbody_system%collider%fragments) ! Calculate the initial energy of the nbody_system without the collisional family - call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) + call self%get_energy_and_momentum(nbody_system, param, phase="before") ! The new body's metadata will be taken from the largest of the two impactor bodies, so we need ! its index in the main pl structure @@ -216,7 +216,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) fragments%vc(:,1) = 0.0_DP ! Get the energy of the system after the collision - call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) + call self%get_energy_and_momentum(nbody_system, param, phase="after") ! Update any encounter lists that have the removed bodies in them so that they instead point to the new body do k = 1, nbody_system%plpl_encounter%nenc diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 50e2bb75a..5b6e041f5 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -121,14 +121,10 @@ module collision real(DP) :: ke_spin_tot !! Spin kinetic energy of all fragments real(DP) :: pe !! Potential energy of all fragments real(DP) :: be !! Binding energy of all fragments - real(DP) :: E_budget !! Kinetic energy budget for computing fragment trajectories - real(DP), dimension(NDIM) :: L_budget !! Angular momentum budget for computing fragment trajectories real(DP), dimension(nbody) :: ke_orbit !! Orbital kinetic energy of each individual fragment real(DP), dimension(nbody) :: ke_spin !! Spin kinetic energy of each individual fragment contains procedure :: reset => collision_util_reset_fragments !! Deallocates all allocatable arrays and sets everything else to 0 - procedure :: get_angular_momentum => collision_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments - procedure :: get_energy => collision_util_get_energy !! Calcualtes the current kinetic energy of the fragments procedure :: set_coordinate_system => collision_util_set_coordinate_fragments !! Sets the coordinate system of the fragments final :: collision_final_fragments !! Finalizer deallocates all allocatables end type collision_fragments @@ -146,27 +142,37 @@ module collision integer(I4B) :: status !! Status flag to pass to the collision list once the collision has been resolved integer(I4B) :: collision_id !! ID number of this collision event - ! For the following variables, index 1 refers to the *entire* n-body nbody_system in its pre-collisional state and index 2 refers to the nbody_system in its post-collisional state - real(DP), dimension(NDIM,2) :: L_orbit !! Before/after orbital angular momentum - real(DP), dimension(NDIM,2) :: L_spin !! Before/after spin angular momentum - real(DP), dimension(NDIM,2) :: L_total !! Before/after total nbody_system angular momentum + ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1 + real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor + real(DP) :: mscale = 1.0_DP !! Mass scale factor + real(DP) :: tscale = 1.0_DP !! Time scale factor + real(DP) :: vscale = 1.0_DP !! Velocity scale factor (a convenience unit that is derived from dscale and tscale) + real(DP) :: Escale = 1.0_DP !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale) + real(DP) :: Lscale = 1.0_DP !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale) + + ! For the following variables, index 1 refers to the *entire* n-body system in its pre-collisional state and index 2 refers to the system in its post-collisional state + real(DP), dimension(NDIM,2) :: L_orbit !! Before/after orbital angular momentum + real(DP), dimension(NDIM,2) :: L_spin !! Before/after spin angular momentum + real(DP), dimension(NDIM,2) :: L_total !! Before/after total nbody_system angular momentum real(DP), dimension(2) :: ke_orbit !! Before/after orbital kinetic energy real(DP), dimension(2) :: ke_spin !! Before/after spin kinetic energy real(DP), dimension(2) :: pe !! Before/after potential energy real(DP), dimension(2) :: be !! Before/after binding energy real(DP), dimension(2) :: te !! Before/after total system energy + contains + procedure :: generate => collision_generate_basic !! Merges the impactors to make a single final body + procedure :: hitandrun => collision_generate_hitandrun !! Merges the impactors to make a single final body + procedure :: merge => collision_generate_merge !! Merges the impactors to make a single final body + procedure :: add_fragments => collision_util_add_fragments_to_collider !! Add fragments to nbody_system + procedure :: get_energy_and_momentum => collision_util_get_energy_and_momentum !! Calculates total nbody_system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) + procedure :: reset => collision_util_reset_system !! Deallocates all allocatables procedure :: setup => collision_util_setup_collider !! Initializer for the encounter collision system and the before/after snapshots procedure :: setup_impactors => collision_util_setup_impactors_collider !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones procedure :: setup_fragments => collision_util_setup_fragments_collider !! Initializer for the fragments of the collision system. - procedure :: add_fragments => collision_util_add_fragments_to_collider !! Add fragments to nbody_system - procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total nbody_system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) - procedure :: reset => collision_util_reset_system !! Deallocates all allocatables - procedure :: set_budgets => collision_util_set_budgets !! Sets the energy and momentum budgets of the fragments based on the collider value procedure :: set_coordinate_system => collision_util_set_coordinate_collider !! Sets the coordinate system of the collisional system - procedure :: generate => collision_generate_basic !! Merges the impactors to make a single final body - procedure :: hitandrun => collision_generate_hitandrun !! Merges the impactors to make a single final body - procedure :: merge => collision_generate_merge !! Merges the impactors to make a single final body + procedure :: set_natural_scale => collision_util_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system. + procedure :: set_original_scale => collision_util_set_original_scale_factors !! Restores dimenional quantities back to the original system units end type collision_basic @@ -380,26 +386,23 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters end subroutine collision_util_add_fragments_to_collider - module subroutine collision_util_get_angular_momentum(self) - implicit none - class(collision_fragments(*)), intent(inout) :: self !! Fragment system object - end subroutine collision_util_get_angular_momentum - - module subroutine collision_util_get_energy(self) + module subroutine collision_util_construct_after_system(collider, nbody_system, param, after_system) + !! Author: David A. Minton + !! + !! Constructs a temporary internal system consisting of active bodies and additional fragments. This internal temporary system is used to calculate system energy with and without fragments implicit none - class(collision_fragments(*)), intent(inout) :: self !! Fragment system object - end subroutine collision_util_get_energy + ! Arguments + class(collision_basic), intent(inout) :: collider !! Collision system object + class(base_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object + class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(base_nbody_system), allocatable, intent(out) :: after_system !! Output temporary swiftest nbody system object + end subroutine collision_util_construct_after_system module subroutine collision_util_reset_fragments(self) implicit none class(collision_fragments(*)), intent(inout) :: self end subroutine collision_util_reset_fragments - module subroutine collision_util_set_budgets(self) - implicit none - class(collision_basic), intent(inout) :: self !! Collision system object - end subroutine collision_util_set_budgets - module subroutine collision_util_set_coordinate_collider(self) implicit none class(collision_basic), intent(inout) :: self !! collisional system @@ -444,14 +447,14 @@ module subroutine collision_util_get_idvalues_snapshot(self, idvals) integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot end subroutine collision_util_get_idvalues_snapshot - module subroutine collision_util_get_energy_momentum(self, nbody_system, param, lbefore) + module subroutine collision_util_get_energy_and_momentum(self, nbody_system, param, phase) use base, only : base_nbody_system, base_parameters implicit none - class(collision_basic), intent(inout) :: self !! Encounter collision system object + class(collision_basic), intent(inout) :: self !! Encounter collision system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current swiftest run configuration parameters - logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the nbody_system, with impactors included and fragments excluded or vice versa - end subroutine collision_util_get_energy_momentum + character(len=*), intent(in) :: phase !! One of "before" or "after", indicating which phase of the calculation this needs to be done + end subroutine collision_util_get_energy_and_momentum module subroutine collision_util_index_map(self) implicit none @@ -476,6 +479,17 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) real(DP), intent(in), optional :: t !! Time of snapshot if different from nbody_system time character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision. end subroutine collision_util_snapshot + + module subroutine collision_util_set_natural_scale_factors(self) + implicit none + class(collision_basic), intent(inout) :: self !! collision system object + end subroutine collision_util_set_natural_scale_factors + + module subroutine collision_util_set_original_scale_factors(self) + implicit none + class(collision_basic), intent(inout) :: self !! collision system object + end subroutine collision_util_set_original_scale_factors + end interface contains diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index a6d5c5330..834f15de0 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -312,13 +312,13 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) implicit none ! Arguments class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters with Swiftest additions - real(DP), intent(in) :: t !! Time of collision - integer(I4B), intent(in) :: status !! Status flag to assign to adds + class(base_parameters), intent(inout) :: param !! Current run configuration parameters with Swiftest additions + real(DP), intent(in) :: t !! Time of collision + integer(I4B), intent(in) :: status !! Status flag to assign to adds ! Internals integer(I4B) :: i, ibiggest, ismallest, iother, nimpactors, nfrag - logical, dimension(:), allocatable :: lmask - class(swiftest_pl), allocatable :: plnew, plsub + logical, dimension(:), allocatable :: lmask + class(swiftest_pl), allocatable :: plnew, plsub character(*), parameter :: FRAGFMT = '("Newbody",I0.7)' character(len=NAMELEN) :: newname, origin_type real(DP) :: volume @@ -488,7 +488,8 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) real(DP), intent(in) :: dt !! Current simulation step size integer(I4B), intent(in) :: irec !! Current recursion level ! Internals - real(DP) :: E_orbit_before, E_orbit_after + real(DP) :: E_before, E_after, dLmag + real(DP), dimension(NDIM) :: L_orbit_before, L_orbit_after, L_spin_before, L_spin_after, L_before, L_after, dL_orbit, dL_spin, dL logical :: lplpl_collision character(len=STRMAX) :: timestr integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision @@ -515,7 +516,10 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) ! Get the energy before the collision is resolved if (param%lenergy) then call nbody_system%get_energy_and_momentum(param) - E_orbit_before = nbody_system%te + E_before = nbody_system%te + L_orbit_before(:) = nbody_system%L_orbit(:) + L_spin_before(:) = nbody_system%L_spin(:) + L_before(:) = nbody_system%L_total(:) end if do loop = 1, MAXCASCADE @@ -584,8 +588,17 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) if (param%lenergy) then call nbody_system%get_energy_and_momentum(param) - E_orbit_after = nbody_system%te - nbody_system%E_collisions = nbody_system%E_collisions + (E_orbit_after - E_orbit_before) + E_after = nbody_system%te + nbody_system%E_collisions = nbody_system%E_collisions + (E_after - E_before) + + L_orbit_after(:) = nbody_system%L_orbit(:) + L_spin_after(:) = nbody_system%L_spin(:) + L_after(:) = nbody_system%L_total(:) + + dL_orbit = L_orbit_after(:) - L_orbit_before(:) + dL_spin = L_spin_after(:) - L_spin_before(:) + dL = L_after(:) - L_before(:) + dLmag = .mag.dL end if end associate diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index fda192a5b..5848f35ff 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -17,7 +17,7 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p !! Adds fragments to the temporary system pl object implicit none ! Arguments - class(collision_basic), intent(in) :: self !! Collision system system object + class(collision_basic), intent(in) :: self !! Collision system system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters ! Internals @@ -61,13 +61,77 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p end subroutine collision_util_add_fragments_to_collider + module subroutine collision_util_construct_after_system(collider, nbody_system, param, after_system) + !! Author: David A. Minton + !! + !! Constructs a temporary internal system consisting of active bodies with impacting bodies replaced with additional fragments. This is used to compute the energy and momentum + !! of the system after the collision has been resolved + implicit none + ! Arguments + class(collision_basic), intent(inout) :: collider !! Collision system object + class(base_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object + class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(base_nbody_system), allocatable, intent(out) :: after_system !! Output temporary swiftest nbody system object + ! Internals + integer(I4B) :: i, status + class(swiftest_nbody_system), allocatable :: tmpsys + class(swiftest_parameters), allocatable :: tmpparam + + select type(nbody_system) + class is (swiftest_nbody_system) + select type(param) + class is (swiftest_parameters) + associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, cb => nbody_system%cb) + + ! Update barycentric vector values + do concurrent(i = 1:nfrag) + fragments%rb(:,i) = fragments%rc(:,i) + impactors%rbcom(:) + fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) + end do + + ! Set up a new temporary system based on the original + allocate(tmpparam, source=param) + call swiftest_util_setup_construct_system(tmpsys, tmpparam) + + if (allocated(tmpsys%cb)) deallocate(tmpsys%cb) + allocate(tmpsys%cb, source=cb) + + if (allocated(tmpsys%pl)) deallocate(tmpsys%pl) + allocate(tmpsys%pl, source=pl) + + if (allocated(tmpsys%collider)) deallocate(tmpsys%collider) + allocate(tmpsys%collider, source=collider) + call tmpsys%collider%set_original_scale() + + select case(impactors%regime) + case(COLLRESOLVE_REGIME_DISRUPTION) + status = DISRUPTED + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + status = SUPERCATASTROPHIC + case(COLLRESOLVE_REGIME_HIT_AND_RUN) + status = HIT_AND_RUN_DISRUPT + end select + + call collision_resolve_mergeaddsub(tmpsys, tmpparam, nbody_system%t, status) + call tmpsys%pl%rearray(tmpsys, tmpparam) + + call move_alloc(tmpsys, after_system) + + end associate + end select + end select + + return + end subroutine collision_util_construct_after_system + + module subroutine collision_util_get_idvalues_snapshot(self, idvals) !! author: David A. Minton !! !! Returns an array of all id values saved in this snapshot implicit none ! Arguments - class(collision_snapshot), intent(in) :: self !! Fraggle snapshot object + class(collision_snapshot), intent(in) :: self !! Collision snapshot object integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot ! Internals integer(I4B) :: npl_before, ntp_before, npl_after, ntp_after, ntot, nlo, nhi @@ -114,141 +178,73 @@ module subroutine collision_util_get_idvalues_snapshot(self, idvals) end subroutine collision_util_get_idvalues_snapshot - module subroutine collision_util_get_angular_momentum(self) - !! Author: David A. Minton - !! - !! Calculates the current angular momentum of the fragments - implicit none - ! Arguments - class(collision_fragments(*)), intent(inout) :: self !! Fraggle fragment system object - ! Internals - integer(I4B) :: i - - associate(fragments => self, nfrag => self%nbody) - - do i = 1, nfrag - fragments%L_orbit(:,i) = fragments%mass(i) * (fragments%rc(:,i) .cross. fragments%vc(:, i)) - fragments%L_spin(:,i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i) * fragments%rot(:,i) - end do - - fragments%L_orbit_tot(:) = sum(fragments%L_orbit, dim=2) - fragments%L_spin_tot(:) = sum(fragments%L_spin, dim=2) - end associate - - return - end subroutine collision_util_get_angular_momentum - - - module subroutine collision_util_get_energy(self) - !! Author: David A. Minton - !! - !! Calculates the current energy of the fragments - implicit none - ! Argument - class(collision_fragments(*)), intent(inout) :: self !! Fragment system object - ! Internals - integer(I4B) :: i,j - real(DP), dimension(self%nbody) :: pepl - - associate(fragments => self, nfrag => self%nbody) - - do concurrent(i = 1:nfrag) - fragments%ke_orbit(i) = fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) - fragments%ke_spin(i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) ) - end do - - fragments%ke_orbit(:) = fragments%ke_orbit(:) / 2 - fragments%ke_spin(:) = fragments%ke_spin(:) / 2 - fragments%ke_orbit_tot = sum(fragments%ke_orbit(:)) - fragments%ke_spin_tot = sum(fragments%ke_spin(:)) - - fragments%pe = 0.0_DP - do i = 1, nfrag - do concurrent(j = i+1:nfrag) - pepl(j) = - (fragments%Gmass(i) * fragments%mass(j)) / .mag.(fragments%rc(:,i) - fragments%rc(:,j)) - end do - fragments%pe = fragments%pe + sum(pepl(i+1:nfrag)) - end do - fragments%be = -sum(3*fragments%Gmass(:)*fragments%mass(:)/(5*fragments%radius(:))) - - end associate - - return - end subroutine collision_util_get_energy - - - module subroutine collision_util_get_energy_momentum(self, nbody_system, param, lbefore) + module subroutine collision_util_get_energy_and_momentum(self, nbody_system, param, phase) !! Author: David A. Minton !! - !! Calculates total system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) + !! Calculates total system energy in either the pre-collision outcome state (phase = "before") or the post-collision outcome state (lbefore = .false.) !! This subrourtine works by building a temporary internal massive body object out of the non-excluded bodies and optionally with fragments appended. !! This will get passed to the energy calculation subroutine so that energy is computed exactly the same way is it is in the main program. - !! This will temporarily expand the massive body object in a temporary system object called tmpsys to feed it into symba_energy + !! This will temporarily expand the massive body object in a temporary system object called after_system to feed it into symba_energy implicit none ! Arguments - class(collision_basic), intent(inout) :: self !! Encounter collision system object - class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current swiftest run configuration parameters - logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the nbody_system, with impactors included and fragments excluded or vice versa + class(collision_basic), intent(inout) :: self !! Encounter collision system object + class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + character(len=*), intent(in) :: phase !! One of "before" or "after", indicating which phase of the calculation this needs to be done ! Internals - integer(I4B) :: stage,i,j - real(DP), dimension(NDIM) :: L_orbit, L_spin - real(DP) :: ke_orbit, ke_spin, pe, be + class(base_nbody_system), allocatable :: after_system + integer(I4B) :: i select type(nbody_system) class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody, pl => nbody_system%pl, cb => nbody_system%cb) - - - if (lbefore) then - L_orbit(:) = sum(impactors%L_orbit(:,:), dim=2) - L_spin(:) = sum(impactors%L_spin(:,:), dim=2) - ke_orbit = 0.0_DP - ke_spin = 0.0_DP - do concurrent(i = 1:2) - ke_orbit = ke_orbit + impactors%mass(i) * dot_product(impactors%vc(:,i), impactors%vc(:,i)) - ke_spin = ke_spin + impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i) * dot_product(impactors%rot(:,i), impactors%rot(:,i)) - end do - ke_orbit = ke_orbit / 2 - ke_spin = ke_spin / 2 - - pe = -(impactors%Gmass(1) * impactors%mass(2)) / .mag.(impactors%rc(:,2) - impactors%rc(:,1)) - be = -sum(3*impactors%Gmass(:)*impactors%mass(:)/(5*impactors%radius(:))) - - else - call fragments%get_angular_momentum() - L_orbit(:) = fragments%L_orbit_tot(:) - L_spin(:) = fragments%L_spin_tot(:) - - call fragments%get_energy() - ke_orbit = fragments%ke_orbit_tot - ke_spin = fragments%ke_spin_tot - pe = fragments%pe - be = fragments%be - + if (phase == "before") then + call nbody_system%get_energy_and_momentum(param) + self%L_orbit(:,1) = nbody_system%L_orbit(:) / self%Lscale + self%L_spin(:,1) = nbody_system%L_spin(:) / self%Lscale + self%L_total(:,1) = nbody_system%L_total(:) / self%Lscale + self%ke_orbit(1) = nbody_system%ke_orbit / self%Escale + self%ke_spin(1) = nbody_system%ke_spin / self%Escale + self%pe(1) = nbody_system%pe / self%Escale + self%be(1) = nbody_system%be / self%Escale + self%te(1) = nbody_system%te / self%Escale + else if (phase == "after") then + call collision_util_construct_after_system(self, nbody_system, param, after_system) + select type(after_system) + class is (swiftest_nbody_system) + call after_system%get_energy_and_momentum(param) + self%L_orbit(:,2) = after_system%L_orbit(:) / self%Lscale + self%L_spin(:,2) = after_system%L_spin(:) / self%Lscale + self%L_total(:,2) = after_system%L_total(:) / self%Lscale + self%ke_orbit(2) = after_system%ke_orbit / self%Escale + self%ke_spin(2) = after_system%ke_spin / self%Escale + self%pe(2) = after_system%pe / self%Escale + self%be(2) = after_system%be / self%Escale + self%te(2) = after_system%te / self%Escale + + do concurrent(i = 1:nfrag) + fragments%ke_orbit(i) = 0.5_DP * fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) + fragments%ke_spin(i) = 0.5_DP * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i), fragments%rot(:,i)) + fragments%L_orbit(:,i) = fragments%mass(i) * fragments%rc(:,i) .cross. fragments%vc(:,i) + fragments%L_spin(:,i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * fragments%rot(:,i) + end do + call swiftest_util_get_potential_energy(nfrag, [(.true., i = 1, nfrag)], after_system%cb%Gmass, fragments%Gmass, fragments%mass, fragments%rb, fragments%pe) + fragments%be = sum(-3*fragments%Gmass(:)*fragments%mass(:)/(5*fragments%radius(:))) + fragments%L_orbit_tot(:) = sum(fragments%L_orbit(:,:),dim=2) + fragments%L_spin_tot(:) = sum(fragments%L_spin(:,:),dim=2) + fragments%ke_orbit_tot = sum(fragments%ke_orbit(:)) + fragments%ke_spin_tot = sum(fragments%ke_spin(:)) + end select end if - ! Calculate the current fragment energy and momentum balances - if (lbefore) then - stage = 1 - else - stage = 2 - end if - self%L_orbit(:,stage) = L_orbit(:) - self%L_spin(:,stage) = L_spin(:) - self%L_total(:,stage) = L_orbit(:) + L_spin(:) - self%ke_orbit(stage) = ke_orbit - self%ke_spin(stage) = ke_spin - self%pe(stage) = pe - self%be(stage) = be - self%te(stage) = ke_orbit + ke_spin + pe + be + end associate end select end select return - end subroutine collision_util_get_energy_momentum + end subroutine collision_util_get_energy_and_momentum module subroutine collision_util_index_map(self) @@ -373,25 +369,6 @@ module subroutine collision_util_reset_system(self) end subroutine collision_util_reset_system - module subroutine collision_util_set_budgets(self) - !! author: David A. Minton - !! - !! Sets the energy and momentum budgets of the fragments based on the collider values and the before/after values of energy and momentum - implicit none - ! Arguments - class(collision_basic), intent(inout) :: self !! Fraggle collision system object - - associate(impactors => self%impactors, fragments => self%fragments) - - fragments%L_budget(:) = self%L_total(:,1) - fragments%E_budget = self%te(1) - impactors%Qloss - - end associate - - return - end subroutine collision_util_set_budgets - - module subroutine collision_util_set_coordinate_collider(self) !! author: David A. Minton !! @@ -568,10 +545,8 @@ module subroutine collision_util_setup_fragments_collider(self, nfrag) self%fragments%vmag(:) = 0.0_DP self%fragments%L_orbit_tot(:) = 0.0_DP self%fragments%L_spin_tot(:) = 0.0_DP - self%fragments%L_budget(:) = 0.0_DP self%fragments%ke_orbit_tot = 0.0_DP self%fragments%ke_spin_tot = 0.0_DP - self%fragments%E_budget = 0.0_DP return end subroutine collision_util_setup_fragments_collider @@ -580,7 +555,7 @@ end subroutine collision_util_setup_fragments_collider module subroutine collision_util_shift_vector_to_origin(m_frag, vec_frag) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! - !! Adjusts the position or velocity of the fragments as needed to align them with the origin + !! Adjusts the position or velocity of the fragments as needed to align them with the center of mass origin implicit none ! Arguments real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses @@ -737,4 +712,130 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) end subroutine collision_util_snapshot + module subroutine collision_util_set_natural_scale_factors(self) + !! author: David A. Minton + !! + !! Scales dimenional quantities to ~O(1) with respect to the collisional system. + !! This scaling makes it it easier to converge on a solution without having floating point issues + implicit none + ! Arguments + class(collision_basic), intent(inout) :: self !! Collision system object + ! Internals + integer(I4B) :: i + real(DP) :: vesc + + associate(collider => self, fragments => self%fragments, impactors => self%impactors) + ! Set primary scale factors (mass, length, and time) based on the impactor properties at the time of collision + collider%mscale = minval(fragments%mass(:)) + collider%dscale = minval(fragments%radius(:)) + + vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) + collider%tscale = collider%dscale / vesc + + ! Set secondary scale factors for convenience + collider%vscale = collider%dscale / collider%tscale + collider%Escale = collider%mscale * collider%vscale**2 + collider%Lscale = collider%mscale * collider%dscale * collider%vscale + + ! Scale all dimensioned quantities of impactors and fragments + impactors%rbcom(:) = impactors%rbcom(:) / collider%dscale + impactors%vbcom(:) = impactors%vbcom(:) / collider%vscale + impactors%rbimp(:) = impactors%rbimp(:) / collider%dscale + impactors%rb(:,:) = impactors%rb(:,:) / collider%dscale + impactors%vb(:,:) = impactors%vb(:,:) / collider%vscale + impactors%rc(:,:) = impactors%rc(:,:) / collider%dscale + impactors%vc(:,:) = impactors%vc(:,:) / collider%vscale + impactors%mass(:) = impactors%mass(:) / collider%mscale + impactors%Gmass(:) = impactors%Gmass(:) / (collider%dscale**3/collider%tscale**2) + impactors%Mcb = impactors%Mcb / collider%mscale + impactors%radius(:) = impactors%radius(:) / collider%dscale + impactors%L_spin(:,:) = impactors%L_spin(:,:) / collider%Lscale + impactors%L_orbit(:,:) = impactors%L_orbit(:,:) / collider%Lscale + + do concurrent(i = 1:2) + impactors%rot(:,i) = impactors%L_spin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) + end do + + fragments%mtot = fragments%mtot / collider%mscale + fragments%mass(:) = fragments%mass(:) / collider%mscale + fragments%Gmass(:) = fragments%Gmass(:) / (collider%dscale**3/collider%tscale**2) + fragments%radius(:) = fragments%radius(:) / collider%dscale + impactors%Qloss = impactors%Qloss / collider%Escale + end associate + + return + end subroutine collision_util_set_natural_scale_factors + + + module subroutine collision_util_set_original_scale_factors(self) + !! author: David A. Minton + !! + !! Restores dimenional quantities back to the system units + use, intrinsic :: ieee_exceptions + implicit none + ! Arguments + class(collision_basic), intent(inout) :: self !! Fragment system object + ! Internals + integer(I4B) :: i + logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes + + call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily + call ieee_set_halting_mode(IEEE_ALL,.false.) + + associate(collider => self, fragments => self%fragments, impactors => self%impactors) + + ! Restore scale factors + impactors%rbcom(:) = impactors%rbcom(:) * collider%dscale + impactors%vbcom(:) = impactors%vbcom(:) * collider%vscale + impactors%rbimp(:) = impactors%rbimp(:) * collider%dscale + + impactors%mass = impactors%mass * collider%mscale + impactors%Gmass(:) = impactors%Gmass(:) * (collider%dscale**3/collider%tscale**2) + impactors%Mcb = impactors%Mcb * collider%mscale + impactors%mass_dist = impactors%mass_dist * collider%mscale + impactors%radius = impactors%radius * collider%dscale + impactors%rb = impactors%rb * collider%dscale + impactors%vb = impactors%vb * collider%vscale + impactors%rc = impactors%rc * collider%dscale + impactors%vc = impactors%vc * collider%vscale + impactors%L_spin = impactors%L_spin * collider%Lscale + impactors%L_orbit = impactors%L_orbit * collider%Lscale + do concurrent(i = 1:2) + impactors%rot(:,i) = impactors%L_spin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) + end do + + fragments%mtot = fragments%mtot * collider%mscale + fragments%mass(:) = fragments%mass(:) * collider%mscale + fragments%Gmass(:) = fragments%Gmass(:) * (collider%dscale**3/collider%tscale**2) + fragments%radius(:) = fragments%radius(:) * collider%dscale + fragments%rot(:,:) = fragments%rot(:,:) / collider%tscale + fragments%rc(:,:) = fragments%rc(:,:) * collider%dscale + fragments%vc(:,:) = fragments%vc(:,:) * collider%vscale + fragments%rb(:,:) = fragments%rb(:,:) * collider%dscale + fragments%vb(:,:) = fragments%vb(:,:) * collider%vscale + + impactors%Qloss = impactors%Qloss * collider%Escale + + collider%L_orbit(:,:) = collider%L_orbit(:,:) * collider%Lscale + collider%L_spin(:,:) = collider%L_spin(:,:) * collider%Lscale + collider%L_total(:,:) = collider%L_total(:,:) * collider%Lscale + collider%ke_orbit(:) = collider%ke_orbit(:) * collider%Escale + collider%ke_spin(:) = collider%ke_spin(:) * collider%Escale + collider%pe(:) = collider%pe(:) * collider%Escale + collider%be(:) = collider%be(:) * collider%Escale + collider%te(:) = collider%te(:) * collider%Escale + + collider%mscale = 1.0_DP + collider%dscale = 1.0_DP + collider%vscale = 1.0_DP + collider%tscale = 1.0_DP + collider%Lscale = 1.0_DP + collider%Escale = 1.0_DP + end associate + call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) + + return + end subroutine collision_util_set_original_scale_factors + + end submodule s_collision_util \ No newline at end of file diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index cadda571f..df6fb6d5f 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -21,9 +21,9 @@ module subroutine fraggle_generate(self, nbody_system, param, t) implicit none ! Arguments class(collision_fraggle), intent(inout) :: self - class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Time of collision + class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Time of collision ! Internals integer(I4B) :: i, ibiggest, nfrag character(len=STRMAX) :: message @@ -55,7 +55,6 @@ module subroutine fraggle_generate(self, nbody_system, param, t) ! Populate the list of new bodies nfrag = fragments%nbody write(message, *) nfrag - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") select case(impactors%regime) case(COLLRESOLVE_REGIME_DISRUPTION) status = DISRUPTED @@ -127,13 +126,12 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call self%set_natural_scale() call fragments%reset() lfailure_local = .false. - call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) - call self%set_budgets() + call self%get_energy_and_momentum(nbody_system, param, phase="before") self%fail_scale = fail_scale_initial call fraggle_generate_pos_vec(self) - call fraggle_generate_rot_vec(self) - call fraggle_generate_vel_vec(self,lfailure_local) - call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) + call fraggle_generate_rot_vec(self, nbody_system, param) + call fraggle_generate_vel_vec(self, nbody_system, param, lfailure_local) + call self%get_energy_and_momentum(nbody_system, param, phase="after") call self%set_original_scale() dE = self%te(2) - self%te(1) @@ -165,12 +163,12 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) !! implicit none ! Arguments - class(collision_fraggle), intent(inout) :: self !! Collision system object + class(collision_fraggle), intent(inout) :: self !! Collision system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(base_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Time of collision ! Result - integer(I4B) :: status !! Status flag assigned to this outcome + integer(I4B) :: status !! Status flag assigned to this outcome ! Internals integer(I4B) :: i, ibiggest, jtarg, jproj, nfrag logical :: lpure @@ -239,6 +237,7 @@ module subroutine fraggle_generate_pos_vec(collider) real(DP), dimension(NDIM,2) :: fragment_cloud_center real(DP), dimension(2) :: fragment_cloud_radius, rdistance logical, dimension(collider%fragments%nbody) :: loverlap + real(DP), dimension(collider%fragments%nbody) :: mass_rscale integer(I4B) :: i, j, loop logical :: lcat, lhitandrun integer(I4B), parameter :: MAXLOOP = 20000 @@ -254,6 +253,13 @@ module subroutine fraggle_generate_pos_vec(collider) ! Place the first two bodies at the centers of the two fragment clouds, but be sure they are sufficiently far apart to avoid overlap loverlap(:) = .true. rdistance(:) = rdistance_scale_factor * .mag.impactors%vc(:,:) + ! Give the fragment positions a random value that is scaled with fragment mass so that the more massive bodies tend to be closer to the impact point + ! Later, velocities will be scaled such that the farther away a fragment is placed from the impact point, the higher will its velocity be. + call random_number(mass_rscale) + mass_rscale(:) = (mass_rscale(:) + 1.0_DP) / 2 + mass_rscale(:) = mass_rscale(:) * (fragments%mtot / fragments%mass(:))**(0.125_DP) ! The power is arbitrary. It just gives the velocity a small mass dependence + mass_rscale(:) = mass_rscale(:) / maxval(mass_rscale(:)) + do loop = 1, MAXLOOP if (.not.any(loverlap(:))) exit fragment_cloud_center(:,1) = impactors%rc(:,1) - rdistance(1) * impactors%bounce_unit(:) @@ -269,8 +275,9 @@ module subroutine fraggle_generate_pos_vec(collider) fragment_cloud_radius(1) = .mag.(fragment_cloud_center(:,1) - impactors%rbimp(:)) fragment_cloud_radius(2) = .mag.(fragment_cloud_center(:,2) - impactors%rbimp(:)) end if + do concurrent(i = 1:nfrag, loverlap(i)) - if (i < 3) then + if (i <= 2) then fragments%rc(:,i) = fragment_cloud_center(:,i) else ! Make a random cloud @@ -278,12 +285,12 @@ module subroutine fraggle_generate_pos_vec(collider) ! Make the fragment cloud symmertic about 0 fragments%rc(:,i) = 2 *(fragments%rc(:,i) - 0.5_DP) - + j = fragments%origin_body(i) ! Scale the cloud size - fragments%rc(:,i) = fragment_cloud_radius(j) * fragments%rc(:,i) - + fragments%rc(:,i) = fragment_cloud_radius(j) * mass_rscale(:) * .unit.fragments%rc(:,i) + ! Shift to the cloud center coordinates fragments%rc(:,i) = fragments%rc(:,i) + fragment_cloud_center(:,j) end if @@ -307,27 +314,25 @@ module subroutine fraggle_generate_pos_vec(collider) fragments%rb(:,i) = fragments%rc(:,i) + impactors%rbcom(:) end do - impactors%rbcom(:) = 0.0_DP - do concurrent(i = 1:nfrag) - impactors%rbcom(:) = impactors%rbcom(:) + fragments%mass(i) * fragments%rb(:,i) - end do - impactors%rbcom(:) = impactors%rbcom(:) / fragments%mtot end associate return end subroutine fraggle_generate_pos_vec - module subroutine fraggle_generate_rot_vec(collider) + module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Calculates the spins of a collection of fragments such that they conserve angular momentum implicit none ! Arguments - class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object + class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object + class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals real(DP), dimension(NDIM) :: Lbefore, Lafter, L_spin, rotdir real(DP) :: v_init, v_final, mass_init, mass_final, rotmag + real(DP), parameter :: random_scale_factor = 0.1_DP !! The relative scale factor to apply to the random component of the rotation vector integer(I4B) :: i associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) @@ -350,7 +355,7 @@ module subroutine fraggle_generate_rot_vec(collider) rotdir = rotdir - 0.5_DP rotdir = .unit. rotdir call random_number(rotmag) - rotmag = 0.01_DP * rotmag * .mag. (Lbefore(:) - Lafter(:)) / ((nfrag - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) + rotmag = random_scale_factor * rotmag * .mag. (Lbefore(:) - Lafter(:)) / ((nfrag - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) fragments%rot(:,i) = rotmag * rotdir end do @@ -360,7 +365,7 @@ module subroutine fraggle_generate_rot_vec(collider) end subroutine fraggle_generate_rot_vec - module subroutine fraggle_generate_vel_vec(collider, lfailure) + module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailure) !! Author: David A. Minton !! !! Generates an initial velocity distribution. For disruptions, the velocity magnitude is set to be @@ -368,163 +373,151 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) !! 2x the escape velocity of the smallest of the two bodies. implicit none ! Arguments - class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object - logical, intent(out) :: lfailure !! Did the velocity computation fail? + class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object + class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + logical, intent(out) :: lfailure !! Did the velocity computation fail? ! Internals - integer(I4B) :: i, j, loop, try, istart, n, ndof, nfrag + integer(I4B) :: i, j, loop, try, istart, nfrag logical :: lhitandrun, lsupercat - real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual, vshear, vunit - real(DP) :: vmag, vesc, rotmag, E_residual, ke_per_dof, ke_tot, E_residual_min + real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual + real(DP) :: vmag, vesc, E_residual, ke_avail, ke_remove, E_residual_min, fscale, f_spin, f_orbit integer(I4B), dimension(collider%fragments%nbody) :: vsign - real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale, ke_avail + real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove + real(DP), parameter :: vmin_initial_factor = 1.5_DP ! For the initial "guess" of fragment velocities, this is the maximum velocity relative to escape velocity that the fragments will have + real(DP), parameter :: vmax_initial_factor = 3.0_DP ! For the initial "guess" of fragment velocities, this is the maximum velocity relative to escape velocity that the fragments will have integer(I4B), parameter :: MAXLOOP = 100 - integer(I4B), parameter :: MAXTRY = 100 - real(DP), parameter :: TOL = 1e-6 - class(collision_fragments(:)), allocatable :: fragments + integer(I4B), parameter :: MAXTRY = 20 + real(DP), parameter :: TOL = 1.0_DP + class(collision_fraggle), allocatable :: collider_local + character(len=STRMAX) :: message associate(impactors => collider%impactors) nfrag = collider%fragments%nbody lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) lsupercat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - allocate(fragments, source=collider%fragments) - - ! The fragments will be divided into two "clouds" based on identified origin body. - ! These clouds will collectively travel like two impactors bouncing off of each other. - where(fragments%origin_body(:) == 1) - vsign(:) = -1 - elsewhere - vsign(:) = 1 - end where - - ! The minimum fragment velocity will be set by the escape velocity - if (lhitandrun) then - vesc = sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) - else - vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) - end if - - ! Scale the magnitude of the velocity by the distance from the impact point - ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct - do concurrent(i = 1:nfrag) - rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) - vscale(i) = .mag. rimp(:) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1))) - end do - vscale(:) = vscale(:)/minval(vscale(:)) + allocate(collider_local, source=collider) + associate(fragments => collider_local%fragments) - ! Give the fragment velocities a random value that is scaled with fragment mass - call random_number(mass_vscale) - mass_vscale(:) = (mass_vscale(:) + 1.0_DP) / 2 - mass_vscale(:) = mass_vscale(:) * (fragments%mtot / fragments%mass(:))**(0.125_DP) ! The power is arbitrary. It just gives the velocity a small mass dependence - mass_vscale(:) = mass_vscale(:) / minval(mass_vscale(:)) + ! The fragments will be divided into two "clouds" based on identified origin body. + ! These clouds will collectively travel like two impactors bouncing off of each other. + where(fragments%origin_body(:) == 1) + vsign(:) = -1 + elsewhere + vsign(:) = 1 + end where - ! Set the velocities of all fragments using all of the scale factors determined above - do concurrent(i = 1:nfrag) - j = fragments%origin_body(i) - vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j)) + ! The minimum fragment velocity will be set by the escape velocity if (lhitandrun) then - if (i == 1) then - fragments%vc(:,1) = impactors%vc(:,1) - else - vmag = .mag.impactors%vc(:,2) / (maxval(mass_vscale(:) * maxval(vscale(:)))) - fragments%vc(:,i) = vmag * mass_vscale(i) * vscale(i) * impactors%bounce_unit(:) * vsign(i) + vrot(:) - end if + vesc = sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) else - ! Add more velocity dispersion to disruptions vs hit and runs. - vmag = vesc * vscale(i) * mass_vscale(i) - rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) - vimp_unit(:) = .unit. rimp(:) - fragments%vc(:,i) = vmag * (impactors%bounce_unit(:) + vimp_unit(:)) * vsign(i) + vrot(:) + vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) end if - end do - if (lhitandrun) then - istart = 2 - else - istart = 1 - end if - call fragments%set_coordinate_system() - E_residual_min = -huge(1.0_DP) - outer: do try = 1, MAXTRY - do loop = 1, MAXLOOP - call fragments%get_energy() - E_residual = fragments%E_budget - (fragments%ke_orbit_tot + fragments%ke_spin_tot + fragments%pe + fragments%be) - if ((abs(E_residual) < abs(E_residual_min)) .or. ((E_residual >= 0.0_DP) .and. (E_residual_min < 0.0_DP))) then ! This is our best case so far. Save it for posterity - if (allocated(collider%fragments)) deallocate(collider%fragments) - allocate(collider%fragments, source=fragments) - E_residual_min = E_residual - if ((E_residual > 0.0_DP) .and. (E_residual < TOL * fragments%E_budget)) exit outer + ! Scale the magnitude of the velocity by the distance from the impact point + ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct + do concurrent(i = 1:nfrag) + rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) + vscale(i) = .mag. rimp(:) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1))) + end do + + ! Set the velocity scale factor to span from vmin/vesc to vmax/vesc + vscale(:) = vscale(:)/minval(vscale(:)) + fscale = log(vmax_initial_factor - vmin_initial_factor + 1.0_DP) / log(maxval(vscale(:))) + vscale(:) = vscale(:)**fscale + vmin_initial_factor - 1.0_DP + + ! Set the velocities of all fragments using all of the scale factors determined above + do concurrent(i = 1:nfrag) + j = fragments%origin_body(i) + vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j)) + if (lhitandrun) then + if (i == 1) then + fragments%vc(:,1) = impactors%vc(:,1) + else + vmag = .mag.impactors%vc(:,2) / maxval(vscale(:)) + fragments%vc(:,i) = vmag * vscale(i) * impactors%bounce_unit(:) * vsign(i) + vrot(:) + end if + else + ! Add more velocity dispersion to disruptions vs hit and runs. + vmag = vesc * vscale(i) + rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) + vimp_unit(:) = .unit. rimp(:) + fragments%vc(:,i) = (impactors%bounce_unit(:) + vimp_unit(:)) * vsign(i) + fragments%vc(:,i) = vmag * .unit.fragments%vc(:,i) !+ vrot(:) end if - ! Make sure we don't take away too much orbital kinetic energy, otherwise the fragment can't escape - ke_avail(:) = fragments%ke_orbit(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%rmag(:) - ke_tot = 0.0_DP - ke_per_dof = -E_residual - do i = 1, 2*(nfrag - istart + 1) - n = count(ke_avail(istart:nfrag) > -E_residual/i) - if (E_residual < 0.0_DP) n = n + count(fragments%ke_spin(istart:nfrag) > -E_residual/i) - if (abs(n * ke_per_dof) > ke_tot) then - ke_per_dof = -E_residual/i - ke_tot = n * ke_per_dof - ndof = i - if (abs(ke_tot) > abs(E_residual)) then - ke_tot = -E_residual - ke_per_dof = ke_tot/n - exit - end if + end do + + ! Hit and run collisions should only affect the runner + if (lhitandrun) then + istart = 2 + else + istart = 1 + end if + + ! Every time the collision-frame velocities are altered, we need to be sure to shift everything back to the center-of-mass frame + call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) + call fragments%set_coordinate_system() + E_residual_min = huge(1.0_DP) + lfailure = .false. + outer: do try = 1, MAXTRY + do loop = 1, MAXLOOP + call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") + ! Check for any residual angular momentum, and if there is any, put it into spin of the largest body + E_residual = collider_local%te(2) - collider_local%te(1) + E_residual = E_residual + impactors%Qloss + L_residual(:) = collider_local%L_total(:,2) - collider_local%L_total(:,1) + fragments%L_spin(:,1) = fragments%L_spin(:,1) - L_residual(:) + fragments%rot(:,1) = fragments%L_spin(:,1) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) + + call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") + E_residual = collider_local%te(2) + impactors%Qloss - collider_local%te(1) + L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) + if ((abs(E_residual) < abs(E_residual_min)) .or. ((E_residual <= 0.0_DP) .and. (E_residual_min >= 0.0_DP))) then ! This is our best case so far. Save it for posterity + if (allocated(collider%fragments)) deallocate(collider%fragments) + allocate(collider%fragments, source=fragments) + E_residual_min = E_residual + if (abs(E_residual) <= tiny(0.0_DP)) exit outer + if ((E_residual < 0.0_DP) .and. (impactors%Qloss / abs(E_residual)) > TOL) exit outer end if - end do - do concurrent(i = istart:nfrag, ke_avail(i) > ke_per_dof) - vmag = max(fragments%vmag(i)**2 - 2*ke_per_dof/fragments%mass(i),vesc**2) - fragments%vmag(i) = sqrt(vmag) - fragments%vc(:,i) = fragments%vmag(i) * .unit.fragments%vc(:,i) - end do - do concurrent(i = istart:nfrag, fragments%ke_spin(i) > ke_per_dof) - rotmag = fragments%rotmag(i)**2 - 2*ke_per_dof/(fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) - rotmag = max(rotmag, 0.0_DP) - fragments%rotmag(i) = sqrt(rotmag) - fragments%rot(:,i) = fragments%rotmag(i) * .unit.fragments%rot(:,i) - end do - call fragments%set_coordinate_system() - if (lsupercat) then - ! Put some of the residual angular momentum into velocity shear. Not too much, or we get some weird trajectories - call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%L_orbit_tot(:) + fragments%L_spin_tot(:)) - do concurrent(i = istart:nfrag) - vunit(:) = .unit. (Lresidual(:) .cross. fragments%r_unit(:,i)) - vshear(:) = vunit(:) * (.mag.Lresidual(:) / ((nfrag-istart+1)*fragments%mass(i) * fragments%rmag(i))) - fragments%vc(:,i) = fragments%vc(:,i) + vshear(:) - end do - end if - ! Check for any residual angular momentum, and if there is any, put it into spin - call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%L_orbit_tot(:) + fragments%L_spin_tot(:)) - do concurrent(i = 1:nfrag) - fragments%L_spin(:,i) = fragments%L_spin(:,i) + Lresidual(:) / nfrag - fragments%rot(:,i) = fragments%L_spin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) - end do + ke_avail = fragments%ke_orbit_tot - 0.5_DP * fragments%mtot * vesc**2 - call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%L_orbit_tot(:) + fragments%L_spin_tot(:)) + ! Remove a constant amount of velocity from the bodies so we don't shift the center of mass and screw up the momentum + f_spin = (fragments%ke_spin_tot )/ (fragments%ke_spin_tot + fragments%ke_orbit_tot) + f_orbit = 1.0_DP - f_spin - end do - ! We didn't converge. Try another configuration and see if we get a better result - call fraggle_generate_pos_vec(collider) - call fraggle_generate_rot_vec(collider) - collider%fail_scale = collider%fail_scale*1.01_DP - end do outer - lfailure = E_residual < 0.0_DP + ke_remove = min(f_orbit * E_residual, ke_avail) + fscale = sqrt((fragments%ke_orbit_tot - f_orbit * ke_remove)/fragments%ke_orbit_tot) + fragments%vc(:,:) = fscale * fragments%vc(:,:) - do concurrent(i = 1:nfrag) - collider%fragments%vb(:,i) = collider%fragments%vc(:,i) + impactors%vbcom(:) - end do + ke_remove = min(E_residual - ke_remove, fragments%ke_spin_tot) + ke_rot_remove(:) = ke_remove * (fragments%ke_spin(:) / fragments%ke_spin_tot) + where(ke_rot_remove(:) > fragments%ke_spin(:)) ke_rot_remove(:) = fragments%ke_spin(:) + do concurrent(i = 1:nfrag) + fscale = sqrt((fragments%ke_spin(i) - ke_rot_remove(i))/fragments%ke_spin(i)) + fragments%rot(:,i) = fscale * fragments%rot(:,i) + end do - impactors%vbcom(:) = 0.0_DP - do concurrent(i = 1:nfrag) - impactors%vbcom(:) = impactors%vbcom(:) + collider%fragments%mass(i) * collider%fragments%vb(:,i) - end do - impactors%vbcom(:) = impactors%vbcom(:) / collider%fragments%mtot + ! Update the unit vectors and magnitudes for the fragments based on their new orbits and rotations + call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) + call fragments%set_coordinate_system() + + end do + ! We didn't converge. Try another configuration and see if we get a better result + call fraggle_generate_pos_vec(collider_local) + call fraggle_generate_rot_vec(collider_local, nbody_system, param) + collider_local%fail_scale = collider_local%fail_scale*1.01_DP + end do outer + lfailure = E_residual < 0.0_DP + if (lfailure) then + write(message,*) "Fraggle failed to converge after ",try*loop," steps. This collision may have added energy." + else + write(message,*) "Fraggle succeeded after ",try*loop," steps." + end if + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + end associate end associate return end subroutine fraggle_generate_vel_vec diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index b50aa3ac3..b0eb38a76 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -25,29 +25,18 @@ module fraggle type, extends(collision_basic) :: collision_fraggle - ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1 - real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor - real(DP) :: mscale = 1.0_DP !! Mass scale factor - real(DP) :: tscale = 1.0_DP !! Time scale factor - real(DP) :: vscale = 1.0_DP !! Velocity scale factor (a convenience unit that is derived from dscale and tscale) - real(DP) :: Escale = 1.0_DP !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale) - real(DP) :: Lscale = 1.0_DP !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale) real(DP) :: fail_scale !! Scale factor to apply to distance values in the position model when overlaps occur. contains - procedure :: disrupt => fraggle_generate_disrupt !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. - procedure :: generate => fraggle_generate !! A simple disruption models that does not constrain energy loss in collisions - procedure :: hitandrun => fraggle_generate_hitandrun - procedure :: set_mass_dist => fraggle_util_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type - procedure :: set_natural_scale => fraggle_util_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system. - procedure :: set_original_scale => fraggle_util_set_original_scale_factors !! Restores dimenional quantities back to the original system units - procedure :: setup_fragments => fraggle_util_setup_fragments_system !! Initializer for the fragments of the collision system. - procedure :: reset => fraggle_util_reset_system !! Deallocates all allocatables - final :: fraggle_final_system !! Finalizer will deallocate all allocatables + procedure :: disrupt => fraggle_generate_disrupt !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. + procedure :: generate => fraggle_generate !! A simple disruption models that does not constrain energy loss in collisions + procedure :: hitandrun => fraggle_generate_hitandrun !! Generates either a pure hit and run, or one in which the runner is disrupted + procedure :: set_mass_dist => fraggle_util_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type + procedure :: setup_fragments => fraggle_util_setup_fragments_system !! Initializer for the fragments of the collision system. + procedure :: reset => fraggle_util_reset_system !! Deallocates all allocatables + final :: fraggle_final_system !! Finalizer will deallocate all allocatables end type collision_fraggle interface - - module subroutine fraggle_generate(self, nbody_system, param, t) implicit none class(collision_fraggle), intent(inout) :: self !! Fraggle fragment system object @@ -78,15 +67,19 @@ module subroutine fraggle_generate_pos_vec(collider) class(collision_fraggle), intent(inout) :: collider !! Fraggle ollision system object end subroutine fraggle_generate_pos_vec - module subroutine fraggle_generate_rot_vec(collider) + module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) implicit none - class(collision_fraggle), intent(inout) :: collider !! Collision system object + class(collision_fraggle), intent(inout) :: collider !! Collision system object + class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine fraggle_generate_rot_vec - module subroutine fraggle_generate_vel_vec(collider, lfailure) + module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailure) implicit none - class(collision_fraggle), intent(inout) :: collider !! Collision system object - logical, intent(out) :: lfailure !! Did the velocity computation fail? + class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object + class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + logical, intent(out) :: lfailure !! Did the velocity computation fail? end subroutine fraggle_generate_vel_vec module subroutine fraggle_util_setup_fragments_system(self, nfrag) @@ -110,17 +103,6 @@ module subroutine fraggle_util_set_mass_dist(self, param) class(collision_fraggle), intent(inout) :: self !! Fraggle collision object class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters end subroutine fraggle_util_set_mass_dist - - module subroutine fraggle_util_set_natural_scale_factors(self) - implicit none - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - end subroutine fraggle_util_set_natural_scale_factors - - module subroutine fraggle_util_set_original_scale_factors(self) - implicit none - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - end subroutine fraggle_util_set_original_scale_factors - end interface contains diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index ffe57632f..667ab6a15 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -73,6 +73,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) real(DP), dimension(2) :: volume real(DP), dimension(NDIM) :: Ip_avg real(DP) :: mfrag, mremaining, min_mfrag, mtot, mcumul, G + real(DP), dimension(:), allocatable :: mass real(DP), parameter :: BETA = 2.85_DP integer(I4B), parameter :: NFRAGMAX = 100 !! Maximum number of fragments that can be generated integer(I4B), parameter :: NFRAGMIN = 7 !! Minimum number of fragments that can be generated (set by the fraggle_generate algorithm for constraining momentum and energy) @@ -80,6 +81,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) integer(I4B), parameter :: iMlr = 1 integer(I4B), parameter :: iMslr = 2 integer(I4B), parameter :: iMrem = 3 + integer(I4B), dimension(:), allocatable :: ind logical :: flipper associate(impactors => self%impactors) @@ -89,7 +91,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) G = impactors%Gmass(1) / impactors%mass(1) Ip_avg(:) = (impactors%mass(1) * impactors%Ip(:,1) + impactors%mass(2) * impactors%Ip(:,2)) / mtot - if (impactors%mass(1) > impactors%mass(2)) then + if (impactors%mass(1) >= impactors%mass(2)) then jtarg = 1 jproj = 2 else @@ -147,16 +149,17 @@ module subroutine fraggle_util_set_mass_dist(self, param) select type(fragments => self%fragments) class is (collision_fragments(*)) fragments%mtot = mtot + allocate(mass, mold=fragments%mass) ! Make the first two bins the same as the Mlr and Mslr values that came from collision_regime - fragments%mass(1) = impactors%mass_dist(iMlr) - fragments%mass(2) = impactors%mass_dist(iMslr) + mass(1) = impactors%mass_dist(iMlr) + mass(2) = impactors%mass_dist(iMslr) ! Distribute the remaining mass the 3:nfrag bodies following the model SFD given by slope BETA mremaining = impactors%mass_dist(iMrem) do i = iMrem, nfrag mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) - fragments%mass(i) = mfrag + mass(i) = mfrag mremaining = mremaining - mfrag end do @@ -166,12 +169,18 @@ module subroutine fraggle_util_set_mass_dist(self, param) else ! If the remainder is postiive, this means that the number of fragments required by the SFD is larger than our upper limit set by computational expediency. istart = iMslr ! We will increase the mass of the 2:nfrag bodies to compensate, which ensures that the second largest fragment remains the second largest end if - mfrag = 1._DP + mremaining / sum(fragments%mass(istart:nfrag)) - fragments%mass(istart:nfrag) = fragments%mass(istart:nfrag) * mfrag + mfrag = 1._DP + mremaining / sum(mass(istart:nfrag)) + mass(istart:nfrag) = mass(istart:nfrag) * mfrag ! There may still be some small residual due to round-off error. If so, simply add it to the last bin of the mass distribution. - mremaining = fragments%mtot - sum(fragments%mass(1:nfrag)) - fragments%mass(nfrag) = fragments%mass(nfrag) + mremaining + mremaining = fragments%mtot - sum(mass(1:nfrag)) + mass(nfrag) = mass(nfrag) + mremaining + + ! Sort the distribution in descending order by mass so that the largest fragment is always the first + call swiftest_util_sort(-mass, ind) + call swiftest_util_sort_rearrange(mass, ind, nfrag) + fragments%mass(:) = mass(:) + deallocate(mass) fragments%Gmass(:) = G * fragments%mass(:) @@ -223,132 +232,6 @@ module subroutine fraggle_util_set_mass_dist(self, param) end subroutine fraggle_util_set_mass_dist - module subroutine fraggle_util_set_natural_scale_factors(self) - !! author: David A. Minton - !! - !! Scales dimenional quantities to ~O(1) with respect to the collisional system. - !! This scaling makes it easier for the non-linear minimization to converge on a solution - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - ! Internals - integer(I4B) :: i - real(DP) :: vesc - - associate(collider => self, fragments => self%fragments, impactors => self%impactors) - ! Set primary scale factors (mass, length, and time) based on the impactor properties at the time of collision - collider%mscale = minval(fragments%mass(:)) - collider%dscale = minval(fragments%radius(:)) - - vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) - collider%tscale = collider%dscale / vesc - - ! Set secondary scale factors for convenience - collider%vscale = collider%dscale / collider%tscale - collider%Escale = collider%mscale * collider%vscale**2 - collider%Lscale = collider%mscale * collider%dscale * collider%vscale - - ! Scale all dimensioned quantities of impactors and fragments - impactors%rbcom(:) = impactors%rbcom(:) / collider%dscale - impactors%vbcom(:) = impactors%vbcom(:) / collider%vscale - impactors%rbimp(:) = impactors%rbimp(:) / collider%dscale - impactors%rb(:,:) = impactors%rb(:,:) / collider%dscale - impactors%vb(:,:) = impactors%vb(:,:) / collider%vscale - impactors%rc(:,:) = impactors%rc(:,:) / collider%dscale - impactors%vc(:,:) = impactors%vc(:,:) / collider%vscale - impactors%mass(:) = impactors%mass(:) / collider%mscale - impactors%Gmass(:) = impactors%Gmass(:) / (collider%dscale**3/collider%tscale**2) - impactors%Mcb = impactors%Mcb / collider%mscale - impactors%radius(:) = impactors%radius(:) / collider%dscale - impactors%L_spin(:,:) = impactors%L_spin(:,:) / collider%Lscale - impactors%L_orbit(:,:) = impactors%L_orbit(:,:) / collider%Lscale - - do concurrent(i = 1:2) - impactors%rot(:,i) = impactors%L_spin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) - end do - - fragments%mtot = fragments%mtot / collider%mscale - fragments%mass(:) = fragments%mass(:) / collider%mscale - fragments%Gmass(:) = fragments%Gmass(:) / (collider%dscale**3/collider%tscale**2) - fragments%radius(:) = fragments%radius(:) / collider%dscale - impactors%Qloss = impactors%Qloss / collider%Escale - end associate - - return - end subroutine fraggle_util_set_natural_scale_factors - - - module subroutine fraggle_util_set_original_scale_factors(self) - !! author: David A. Minton - !! - !! Restores dimenional quantities back to the system units - use, intrinsic :: ieee_exceptions - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: self !! Fraggle fragment system object - ! Internals - integer(I4B) :: i - logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes - - call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily - call ieee_set_halting_mode(IEEE_ALL,.false.) - - associate(collider => self, fragments => self%fragments, impactors => self%impactors) - - ! Restore scale factors - impactors%rbcom(:) = impactors%rbcom(:) * collider%dscale - impactors%vbcom(:) = impactors%vbcom(:) * collider%vscale - impactors%rbimp(:) = impactors%rbimp(:) * collider%dscale - - impactors%mass = impactors%mass * collider%mscale - impactors%Gmass(:) = impactors%Gmass(:) * (collider%dscale**3/collider%tscale**2) - impactors%Mcb = impactors%Mcb * collider%mscale - impactors%mass_dist = impactors%mass_dist * collider%mscale - impactors%radius = impactors%radius * collider%dscale - impactors%rb = impactors%rb * collider%dscale - impactors%vb = impactors%vb * collider%vscale - impactors%rc = impactors%rc * collider%dscale - impactors%vc = impactors%vc * collider%vscale - impactors%L_spin = impactors%L_spin * collider%Lscale - impactors%L_orbit = impactors%L_orbit * collider%Lscale - do concurrent(i = 1:2) - impactors%rot(:,i) = impactors%L_spin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) - end do - - fragments%mtot = fragments%mtot * collider%mscale - fragments%mass(:) = fragments%mass(:) * collider%mscale - fragments%Gmass(:) = fragments%Gmass(:) * (collider%dscale**3/collider%tscale**2) - fragments%radius(:) = fragments%radius(:) * collider%dscale - fragments%rot(:,:) = fragments%rot(:,:) / collider%tscale - fragments%rc(:,:) = fragments%rc(:,:) * collider%dscale - fragments%vc(:,:) = fragments%vc(:,:) * collider%vscale - fragments%rb(:,:) = fragments%rb(:,:) * collider%dscale - fragments%vb(:,:) = fragments%vb(:,:) * collider%vscale - - impactors%Qloss = impactors%Qloss * collider%Escale - - collider%L_orbit(:,:) = collider%L_orbit(:,:) * collider%Lscale - collider%L_spin(:,:) = collider%L_spin(:,:) * collider%Lscale - collider%L_total(:,:) = collider%L_total(:,:) * collider%Lscale - collider%ke_orbit(:) = collider%ke_orbit(:) * collider%Escale - collider%ke_spin(:) = collider%ke_spin(:) * collider%Escale - collider%pe(:) = collider%pe(:) * collider%Escale - collider%be(:) = collider%be(:) * collider%Escale - collider%te(:) = collider%te(:) * collider%Escale - - collider%mscale = 1.0_DP - collider%dscale = 1.0_DP - collider%vscale = 1.0_DP - collider%tscale = 1.0_DP - collider%Lscale = 1.0_DP - collider%Escale = 1.0_DP - end associate - call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) - - return - end subroutine fraggle_util_set_original_scale_factors - - module subroutine fraggle_util_setup_fragments_system(self, nfrag) !! author: David A. Minton !! diff --git a/src/swiftest/swiftest_driver.f90 b/src/swiftest/swiftest_driver.f90 index 0e9f369df..d75b8a561 100644 --- a/src/swiftest/swiftest_driver.f90 +++ b/src/swiftest/swiftest_driver.f90 @@ -71,7 +71,7 @@ program swiftest_driver associate (system_history => param%system_history) ! If this is a new run, compute energy initial conditions (if energy tracking is turned on) and write the initial conditions to file. - call nbody_system%display_run_information(param, integration_timer, phase="FIRST") + call nbody_system%display_run_information(param, integration_timer, phase="first") if (param%lenergy) then if (param%lrestart) then call nbody_system%get_t0_values(param) @@ -120,7 +120,7 @@ program swiftest_driver ! Dump any remaining history if it exists call nbody_system%dump(param) call system_history%dump(param) - call nbody_system%display_run_information(param, integration_timer, phase="LAST") + call nbody_system%display_run_information(param, integration_timer, phase="last") end associate end associate diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index db3ebe555..0c5dbb85d 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -203,7 +203,7 @@ module subroutine swiftest_io_display_run_information(self, param, integration_t class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters type(walltimer), intent(inout) :: integration_timer !! Object used for computing elapsed wall time - character(len=*), optional, intent(in) :: phase !! One of "FIRST" or "LAST" + character(len=*), optional, intent(in) :: phase !! One of "first" or "last" ! Internals integer(I4B) :: phase_val real(DP) :: tfrac !! Fraction of total simulation time completed @@ -218,9 +218,9 @@ module subroutine swiftest_io_display_run_information(self, param, integration_t phase_val = 1 if (present(phase)) then - if (phase == "FIRST") then + if (phase == "first") then phase_val = 0 - else if (phase == "LAST") then + else if (phase == "last") then phase_val = -1 end if end if diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index fd393dbc1..872981289 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -392,7 +392,7 @@ module swiftest procedure :: init_particle_info => swiftest_util_setup_initialize_particle_info_system !! Initialize the nbody_system from input files ! procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. procedure :: set_msys => swiftest_util_set_msys !! Sets the value of msys from the masses of nbody_system bodies. - procedure :: get_energy_and_momentum => swiftest_util_get_energy_momentum_system !! Calculates the total nbody_system energy and momentum + procedure :: get_energy_and_momentum => swiftest_util_get_energy_and_momentum_system !! Calculates the total nbody_system energy and momentum procedure :: get_idvals => swiftest_util_get_idvalues_system !! Returns an array of all id values in use in the nbody_system procedure :: rescale => swiftest_util_rescale_system !! Rescales the nbody_system into a new set of units procedure :: validate_ids => swiftest_util_valid_id_system !! Validate the numerical ids passed to the nbody_system and save the maximum value @@ -580,7 +580,7 @@ module subroutine swiftest_io_display_run_information(self, param, integration_t class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters type(walltimer), intent(inout) :: integration_timer !! Object used for computing elapsed wall time - character(len=*), optional, intent(in) :: phase !! One of "FIRST" or "LAST" + character(len=*), optional, intent(in) :: phase !! One of "first" or "last" end subroutine swiftest_io_display_run_information module subroutine swiftest_io_dump_param(self, param_file_name) @@ -1319,7 +1319,6 @@ end subroutine swiftest_util_fill_arr_logical end interface interface - pure module subroutine swiftest_util_flatten_eucl_ij_to_k(n, i, j, k) !$omp declare simd(swiftest_util_flatten_eucl_ij_to_k) implicit none @@ -1350,6 +1349,46 @@ module subroutine swiftest_util_flatten_eucl_pltp(self, pl, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine + module subroutine swiftest_util_get_energy_and_momentum_system(self, param) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine swiftest_util_get_energy_and_momentum_system + + module subroutine swiftest_util_get_idvalues_system(self, idvals) + implicit none + class(swiftest_nbody_system), intent(in) :: self !! Encounter snapshot object + integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot + end subroutine swiftest_util_get_idvalues_system + end interface + + interface swiftest_util_get_potential_energy + module subroutine swiftest_util_get_potential_energy_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mass, rb, pe) + implicit none + integer(I4B), intent(in) :: npl + integer(I8B), intent(in) :: nplpl + integer(I4B), dimension(:,:), intent(in) :: k_plpl + logical, dimension(:), intent(in) :: lmask + real(DP), intent(in) :: GMcb + real(DP), dimension(:), intent(in) :: Gmass + real(DP), dimension(:), intent(in) :: mass + real(DP), dimension(:,:), intent(in) :: rb + real(DP), intent(out) :: pe + end subroutine swiftest_util_get_potential_energy_flat + + module subroutine swiftest_util_get_potential_energy_triangular(npl, lmask, GMcb, Gmass, mass, rb, pe) + implicit none + integer(I4B), intent(in) :: npl + logical, dimension(:), intent(in) :: lmask + real(DP), intent(in) :: GMcb + real(DP), dimension(:), intent(in) :: Gmass + real(DP), dimension(:), intent(in) :: mass + real(DP), dimension(:,:), intent(in) :: rb + real(DP), intent(out) :: pe + end subroutine swiftest_util_get_potential_energy_triangular + end interface + + interface module subroutine swiftest_util_get_vals_storage(self, idvals, tvals) class(swiftest_storage(*)), intent(in) :: self !! Swiftest storage object integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values in all snapshots @@ -1471,18 +1510,6 @@ module subroutine swiftest_util_resize_tp(self, nnew) integer(I4B), intent(in) :: nnew !! New size neded end subroutine swiftest_util_resize_tp - module subroutine swiftest_util_get_energy_momentum_system(self, param) - implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine swiftest_util_get_energy_momentum_system - - module subroutine swiftest_util_get_idvalues_system(self, idvals) - implicit none - class(swiftest_nbody_system), intent(in) :: self !! Encounter snapshot object - integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot - end subroutine swiftest_util_get_idvalues_system - module subroutine swiftest_util_set_beg_end_pl(self, rbeg, rend, vbeg) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index aa3f4c4d0..30e3080ad 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1221,7 +1221,7 @@ module subroutine swiftest_util_flatten_eucl_pltp(self, pl, param) end subroutine swiftest_util_flatten_eucl_pltp - module subroutine swiftest_util_get_energy_momentum_system(self, param) + module subroutine swiftest_util_get_energy_and_momentum_system(self, param) !! author: David A. Minton !! !! Compute total nbody_system angular momentum vector and kinetic, potential and total nbody_system energy @@ -1298,9 +1298,9 @@ module subroutine swiftest_util_get_energy_momentum_system(self, param) end if if (param%lflatten_interactions) then - call swiftest_util_get_energy_potential_flat(npl, pl%nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%rb, nbody_system%pe) + call swiftest_util_get_potential_energy(npl, pl%nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%rb, nbody_system%pe) else - call swiftest_util_get_energy_potential_triangular(npl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%rb, nbody_system%pe) + call swiftest_util_get_potential_energy(npl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%rb, nbody_system%pe) end if ! Potential energy from the oblateness term @@ -1329,10 +1329,10 @@ module subroutine swiftest_util_get_energy_momentum_system(self, param) end associate return - end subroutine swiftest_util_get_energy_momentum_system + end subroutine swiftest_util_get_energy_and_momentum_system - subroutine swiftest_util_get_energy_potential_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mass, rb, pe) + module subroutine swiftest_util_get_potential_energy_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mass, rb, pe) !! author: David A. Minton !! !! Compute total nbody_system potential energy @@ -1381,10 +1381,10 @@ subroutine swiftest_util_get_energy_potential_flat(npl, nplpl, k_plpl, lmask, GM pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), lmask(1:npl)) return - end subroutine swiftest_util_get_energy_potential_flat + end subroutine swiftest_util_get_potential_energy_flat - subroutine swiftest_util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, rb, pe) + module subroutine swiftest_util_get_potential_energy_triangular(npl, lmask, GMcb, Gmass, mass, rb, pe) !! author: David A. Minton !! !! Compute total nbody_system potential energy @@ -1427,41 +1427,7 @@ subroutine swiftest_util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass pe = pe + sum(pecb(1:npl), lmask(1:npl)) return - end subroutine swiftest_util_get_energy_potential_triangular - - - module subroutine swiftest_util_index_array(ind_arr, n) - !! author: David A. Minton - !! - !! Creates or resizes an index array of size n where ind_arr = [1, 2, ... n]. - !! This subroutine swiftest_assumes that if ind_arr is already allocated, it is a pre-existing index array of a different size. - implicit none - ! Arguments - integer(I4B), dimension(:), allocatable, intent(inout) :: ind_arr !! Index array. Input is a pre-existing index array where n /= size(ind_arr). Output is a new index array ind_arr = [1, 2, ... n] - integer(I4B), intent(in) :: n !! The new size of the index array - ! Internals - integer(I4B) :: nold, i - integer(I4B), dimension(:), allocatable :: itmp - - if (allocated(ind_arr)) then - nold = size(ind_arr) - if (nold == n) return ! Nothing to do, so go home - else - nold = 0 - end if - - allocate(itmp(n)) - if (n >= nold) then - if (nold > 0) itmp(1:nold) = ind_arr(1:nold) - itmp(nold+1:n) = [(i, i = nold + 1, n)] - call move_alloc(itmp, ind_arr) - else - itmp(1:n) = ind_arr(1:n) - call move_alloc(itmp, ind_arr) - end if - - return - end subroutine swiftest_util_index_array + end subroutine swiftest_util_get_potential_energy_triangular module subroutine swiftest_util_get_idvalues_system(self, idvals) @@ -1554,6 +1520,40 @@ module subroutine swiftest_util_get_vals_storage(self, idvals, tvals) end subroutine swiftest_util_get_vals_storage + module subroutine swiftest_util_index_array(ind_arr, n) + !! author: David A. Minton + !! + !! Creates or resizes an index array of size n where ind_arr = [1, 2, ... n]. + !! This subroutine swiftest_assumes that if ind_arr is already allocated, it is a pre-existing index array of a different size. + implicit none + ! Arguments + integer(I4B), dimension(:), allocatable, intent(inout) :: ind_arr !! Index array. Input is a pre-existing index array where n /= size(ind_arr). Output is a new index array ind_arr = [1, 2, ... n] + integer(I4B), intent(in) :: n !! The new size of the index array + ! Internals + integer(I4B) :: nold, i + integer(I4B), dimension(:), allocatable :: itmp + + if (allocated(ind_arr)) then + nold = size(ind_arr) + if (nold == n) return ! Nothing to do, so go home + else + nold = 0 + end if + + allocate(itmp(n)) + if (n >= nold) then + if (nold > 0) itmp(1:nold) = ind_arr(1:nold) + itmp(nold+1:n) = [(i, i = nold + 1, n)] + call move_alloc(itmp, ind_arr) + else + itmp(1:n) = ind_arr(1:n) + call move_alloc(itmp, ind_arr) + end if + + return + end subroutine swiftest_util_index_array + + module subroutine swiftest_util_index_map_storage(self) !! author: David A. Minton !! @@ -1740,7 +1740,6 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param) end select end select - ! Re-index the encounter list as the index values may have changed if (nenc_old > 0) then nencmin = min(nbody_system%plpl_encounter%nenc, plplenc_old%nenc)