diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 23a0d20b4..fed0191b4 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -122,10 +122,11 @@ module collision 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_kinetic_energy => collision_util_get_kinetic_energy !! Calcualtes the current kinetic energy of the fragments - final :: collision_final_fragments !! Finalizer deallocates all allocatables + 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_kinetic_energy => collision_util_get_kinetic_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 @@ -231,7 +232,6 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) real(DP), intent(in) :: t !! The time of the collision end subroutine collision_generate_bounce - module subroutine collision_generate_hitandrun(self, nbody_system, param, t) implicit none class(collision_basic), intent(inout) :: self @@ -247,7 +247,6 @@ module subroutine collision_generate_merge(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 end subroutine collision_generate_merge - module subroutine collision_io_collider_message(pl, collidx, collider_message) implicit none @@ -406,6 +405,11 @@ module subroutine collision_util_set_coordinate_collider(self) class(collision_basic), intent(inout) :: self !! collisional system end subroutine collision_util_set_coordinate_collider + module subroutine collision_util_set_coordinate_fragments(self) + implicit none + class(collision_fragments(*)), intent(inout) :: self !! Collisional nbody_system + end subroutine collision_util_set_coordinate_fragments + module subroutine collision_util_set_coordinate_impactors(self) implicit none class(collision_impactors), intent(inout) :: self !! collisional system diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 2878f95ec..3659bbc63 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -449,22 +449,39 @@ module subroutine collision_util_set_coordinate_collider(self) call impactors%set_coordinate_system() if (.not.allocated(self%fragments)) return + call fragments%set_coordinate_system() + + + end associate + + return + end subroutine collision_util_set_coordinate_collider + + + module subroutine collision_util_set_coordinate_fragments(self) + !! author: David A. Minton + !! + !! Defines the collisional coordinate nbody_system, including the unit vectors of both the nbody_system and individual fragments. + implicit none + ! Arguments + class(collision_fragments(*)), intent(inout) :: self !! Collisional nbody_system + + associate(fragments => self, nfrag => self%nbody) if ((nfrag == 0) .or. (.not.any(fragments%rc(:,:) > 0.0_DP))) return fragments%rmag(:) = .mag. fragments%rc(:,:) fragments%vmag(:) = .mag. fragments%vc(:,:) fragments%rotmag(:) = .mag. fragments%rot(:,:) - + ! Define the radial, normal, and tangential unit vectors for each individual fragment fragments%r_unit(:,:) = .unit. fragments%rc(:,:) fragments%v_unit(:,:) = .unit. fragments%vc(:,:) fragments%n_unit(:,:) = .unit. (fragments%rc(:,:) .cross. fragments%vc(:,:)) fragments%t_unit(:,:) = -.unit. (fragments%r_unit(:,:) .cross. fragments%n_unit(:,:)) - end associate return - end subroutine collision_util_set_coordinate_collider + end subroutine collision_util_set_coordinate_fragments module subroutine collision_util_set_coordinate_impactors(self) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index f00080f6b..925b04620 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -98,7 +98,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur logical :: lk_plpl, lfailure_local logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes character(len=STRMAX) :: message - real(DP), parameter :: fail_scale_initial = 1.001_DP + real(DP), parameter :: fail_scale_initial = 1.01_DP ! The minimization and linear solvers can sometimes lead to floating point exceptions. Rather than halting the code entirely if this occurs, we ! can simply fail the attempt and try again. So we need to turn off any floating point exception halting modes temporarily @@ -136,23 +136,16 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur ! Compute the "before" energy/momentum and the budgets call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) call self%set_budgets() - do try = 1, MAXTRY - self%fail_scale = (fail_scale_initial)**try - write(message,*) try - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle try " // trim(adjustl(message))) - call fraggle_generate_pos_vec(self) - call fraggle_generate_rot_vec(self) - call fraggle_generate_vel_vec(self,lfailure_local) - if (.not.lfailure_local) exit - end do + 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%set_original_scale() if (lfailure_local) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle fragment generation failed after " // & - trim(adjustl(message)) // " tries. This collision may have gained energy.") + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "This collision may have gained energy.") else - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle fragment generation succeeded after " // & - trim(adjustl(message)) // " tries") + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle fragment generation succeeded.") end if @@ -253,8 +246,8 @@ module subroutine fraggle_generate_pos_vec(collider) logical, dimension(collider%fragments%nbody) :: loverlap integer(I4B) :: i, j, loop logical :: lcat, lhitandrun - integer(I4B), parameter :: MAXLOOP = 10000 - real(DP), parameter :: rdistance_scale_factor = 0.01_DP ! Scale factor to apply to distance scaling in the event of a fail + integer(I4B), parameter :: MAXLOOP = 20000 + real(DP), parameter :: rdistance_scale_factor = 0.1_DP ! Scale factor to apply to distance scaling in the event of a fail associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) @@ -269,6 +262,8 @@ module subroutine fraggle_generate_pos_vec(collider) if (.not.any(loverlap(:))) exit fragment_cloud_center(:,1) = impactors%rc(:,1) - rdistance(1) * impactors%bounce_unit(:) fragment_cloud_center(:,2) = impactors%rc(:,2) + rdistance(2) * impactors%bounce_unit(:) + fragment_cloud_radius(1) = .mag.(fragment_cloud_center(:,1) - impactors%rbimp(:)) + fragment_cloud_radius(2) = .mag.(fragment_cloud_center(:,2) - impactors%rbimp(:)) do concurrent(i = 1:nfrag, loverlap(i)) if (i < 3) then fragments%rc(:,i) = fragment_cloud_center(:,i) @@ -298,7 +293,6 @@ module subroutine fraggle_generate_pos_vec(collider) end do end do rdistance(:) = rdistance(:) * collider%fail_scale - fragment_cloud_radius(:) = fragment_cloud_radius(:) * collider%fail_scale end do call collision_util_shift_vector_to_origin(fragments%mass, fragments%rc) @@ -375,16 +369,19 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) integer(I4B) :: i, j, loop, istart, n, ndof logical :: lhitandrun, lnoncat real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual, vshear, vunit - real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof, ke_tot + real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof, ke_tot, ke_residual_min integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale, ke_avail - integer(I4B), parameter :: MAXLOOP = 100 + integer(I4B), parameter :: MAXLOOP = 2000 real(DP), parameter :: TOL = 1e-2 + class(collision_fragments(:)), allocatable :: fragments - associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) + associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) lnoncat = (impactors%regime /= COLLRESOLVE_REGIME_SUPERCATASTROPHIC) ! For non-catastrophic impacts, make the fragments act like ejecta and point away from the impact point + 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) @@ -414,6 +411,7 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) 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(:)) + ! Set the velocities of all fragments using all of the scale factors determined above do concurrent(i = 1:nfrag) j = fragments%origin_body(i) @@ -439,10 +437,17 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) else istart = 1 end if + ke_residual_min = huge(1.0_DP) do loop = 1, MAXLOOP - call collider%set_coordinate_system() + call fragments%set_coordinate_system() call fragments%get_kinetic_energy() ke_residual = fragments%ke_budget - (fragments%ke_orbit_tot + fragments%ke_spin_tot) + if (abs(ke_residual) < abs(ke_residual_min)) then ! This is our best case so far. Save it for posterity + if (allocated(collider%fragments)) deallocate(collider%fragments) + allocate(collider%fragments, source=fragments) + ke_residual_min = ke_residual + end if + ke_residual_min = max(ke_residual_min,ke_residual) if (ke_residual > 0.0_DP) exit ! 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(:) @@ -475,7 +480,7 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) end do ! Check for any residual angular momentum, and if there is any, put it into spin - call collider%set_coordinate_system() + call fragments%set_coordinate_system() call fragments%get_angular_momentum() Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) do concurrent(i = istart:nfrag)