diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 8253fb12a..54b431287 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -72,7 +72,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa call frag%get_energy_and_momentum(colliders, system, param, lbefore=.true.) ! Start out the fragments close to the initial separation distance. This will be increased if there is any overlap or we fail to find a solution - r_max_start = 1 * norm2(colliders%rb(:,2) - colliders%rb(:,1)) + r_max_start = 2 * norm2(colliders%rb(:,2) - colliders%rb(:,1)) lfailure = .false. try = 1 do while (try < MAXTRY) @@ -177,8 +177,10 @@ subroutine fraggle_generate_pos_vec(frag, colliders, r_max_start) real(DP), intent(in) :: r_max_start !! Initial guess for the starting maximum radial distance of fragments ! Internals real(DP) :: dis, rad, r_max + real(DP), dimension(NDIM) :: runit logical, dimension(:), allocatable :: loverlap integer(I4B) :: i, j + logical :: lfixdir associate(nfrag => frag%nbody) allocate(loverlap(nfrag)) @@ -188,9 +190,13 @@ subroutine fraggle_generate_pos_vec(frag, colliders, r_max_start) r_max = r_max_start rad = sum(colliders%radius(:)) + lfixdir = (frag%regime /= COLLRESOLVE_REGIME_SUPERCATASTROPHIC) ! For this style of impact, make the fragments act like eject and point away from the impact point + + runit(:) = colliders%rb(:,2) - colliders%rb(:,1) + runit(:) = runit(:) / (.mag. runit(:)) + ! We will treat the first two fragments of the list as special cases. They get initialized the maximum distances apart along the original impactor distance vector. ! This is done because in a regular disruption, the first body is the largest, the second the second largest, and the rest are smaller equal-mass fragments. - call random_number(frag%x_coll(:,3:nfrag)) loverlap(:) = .true. do while (any(loverlap(3:nfrag))) @@ -201,6 +207,8 @@ subroutine fraggle_generate_pos_vec(frag, colliders, r_max_start) if (loverlap(i)) then call random_number(frag%x_coll(:,i)) frag%x_coll(:, i) = 2 * (frag%x_coll(:, i) - 0.5_DP) * r_max + frag%x_coll(:, i) = frag%x_coll(:, i) + (frag%rbimp(:) - frag%rbcom(:)) ! Shift the center of the fragment cloud to the impact point rather than the CoM + if (lfixdir .and. dot_product(frag%x_coll(:,i), runit(:)) < 0.0_DP) frag%x_coll(:, i) = -frag%x_coll(:, i) ! Make sure the fragment cloud points away from the impact point end if end do loverlap(:) = .false. @@ -214,12 +222,12 @@ subroutine fraggle_generate_pos_vec(frag, colliders, r_max_start) call fraggle_util_shift_vector_to_origin(frag%mass, frag%x_coll) call frag%set_coordinate_system(colliders) - do i = 1, nfrag + do concurrent(i = 1:nfrag) frag%rb(:,i) = frag%x_coll(:,i) + frag%rbcom(:) end do frag%rbcom(:) = 0.0_DP - do i = 1, nfrag + do concurrent(i = 1:nfrag) frag%rbcom(:) = frag%rbcom(:) + frag%mass(i) * frag%rb(:,i) end do frag%rbcom(:) = frag%rbcom(:) / frag%mtot diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 index 7b3191149..7962e6c25 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/fraggle/fraggle_regime.f90 @@ -27,7 +27,7 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) integer(I4B) :: jtarg, jproj real(DP), dimension(2) :: radius_si, mass_si, density_si real(DP) :: min_mfrag_si, Mcb_si - real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si + real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si, runit real(DP) :: mlr, mslr, mtot, dentot associate(colliders => self) @@ -71,6 +71,11 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) frag%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / frag%mtot frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot + ! Find the point of impact between the two bodies + runit(:) = colliders%rb(:,2) - colliders%rb(:,1) + runit(:) = runit(:) / (.mag. runit(:)) + frag%rbimp(:) = colliders%rb(:,1) + colliders%radius(1) * runit(:) + ! Convert quantities back to the system units and save them into the fragment system frag%mass_dist(:) = (frag%mass_dist(:) / param%MU2KG) frag%Qloss = frag%Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index a8e61130d..45cf41a92 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -185,10 +185,11 @@ module subroutine fraggle_set_coordinate_system(self, colliders) Ltot = colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + colliders%L_spin(:,1) + colliders%L_spin(:,2) frag%y_coll_unit(:) = delta_r(:) / r_col_norm L_mag = .mag.Ltot(:) - if (L_mag > tiny(L_mag)) then + if (L_mag > sqrt(tiny(L_mag))) then frag%z_coll_unit(:) = Ltot(:) / L_mag else - frag%z_coll_unit(:) = 0.0_DP + call random_number(frag%z_coll_unit(:)) + frag%z_coll_unit(:) = frag%z_coll_unit(:) / (.mag.frag%z_coll_unit(:)) end if ! The cross product of the y- by z-axis will give us the x-axis frag%x_coll_unit(:) = frag%y_coll_unit(:) .cross. frag%z_coll_unit(:) @@ -198,6 +199,7 @@ module subroutine fraggle_set_coordinate_system(self, colliders) call random_number(L_sigma(:,:)) ! Randomize the tangential velocity direction. This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, ! otherwise we can get an ill-conditioned system + do concurrent(i = 1:nfrag, frag%rmag(i) > 0.0_DP) frag%v_r_unit(:, i) = frag%x_coll(:, i) / frag%rmag(i) frag%v_n_unit(:, i) = frag%z_coll_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP) @@ -205,6 +207,7 @@ module subroutine fraggle_set_coordinate_system(self, colliders) frag%v_t_unit(:, i) = frag%v_n_unit(:, i) .cross. frag%v_r_unit(:, i) frag%v_t_unit(:, i) = frag%v_t_unit(:, i) / (.mag. frag%v_t_unit(:, i)) end do + end associate return @@ -236,6 +239,7 @@ module subroutine fraggle_set_natural_scale_factors(self, colliders) ! Scale all dimensioned quantities of colliders and fragments frag%rbcom(:) = frag%rbcom(:) / frag%dscale frag%vbcom(:) = frag%vbcom(:) / frag%vscale + frag%rbimp(:) = frag%rbimp(:) / frag%dscale colliders%rb(:,:) = colliders%rb(:,:) / frag%dscale colliders%vb(:,:) = colliders%vb(:,:) / frag%vscale colliders%mass(:) = colliders%mass(:) / frag%mscale @@ -278,6 +282,7 @@ module subroutine fraggle_set_original_scale_factors(self, colliders) ! Restore scale factors frag%rbcom(:) = frag%rbcom(:) * frag%dscale frag%vbcom(:) = frag%vbcom(:) * frag%vscale + frag%rbimp(:) = frag%rbimp(:) * frag%dscale colliders%mass = colliders%mass * frag%mscale colliders%radius = colliders%radius * frag%dscale diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 038b3c1a5..3ed18f32d 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -241,7 +241,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para frag%ke_orbit_before = tmpsys%ke_orbit frag%ke_spin_before = tmpsys%ke_spin frag%pe_before = tmpsys%pe - frag%Etot_before = tmpsys%te + frag%Etot_before = tmpsys%te else frag%Lorbit_after(:) = tmpsys%Lorbit(:) frag%Lspin_after(:) = tmpsys%Lspin(:) @@ -249,7 +249,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para frag%ke_orbit_after = tmpsys%ke_orbit frag%ke_spin_after = tmpsys%ke_spin frag%pe_after = tmpsys%pe - frag%Etot_after = tmpsys%te + frag%Etot_after = tmpsys%te - (frag%pe_after - frag%pe_before) ! Gotta be careful with PE when number of bodies changes. end if end associate diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 8f6cd4310..ebd5212cf 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -54,6 +54,7 @@ module fraggle_classes ! Values in a coordinate frame centered on the collider barycenter and collisional system unit vectors (these are used internally by the fragment generation subroutine) real(DP), dimension(NDIM) :: rbcom !! Center of mass position vector of the collider system in system barycentric coordinates real(DP), dimension(NDIM) :: vbcom !! Velocity vector of the center of mass of the collider system in system barycentric coordinates + real(DP), dimension(NDIM) :: rbimp !! Impact point position vector of the collider system in system barycentric coordinates real(DP), dimension(NDIM) :: x_coll_unit !! x-direction unit vector of collisional system real(DP), dimension(NDIM) :: y_coll_unit !! y-direction unit vector of collisional system real(DP), dimension(NDIM) :: z_coll_unit !! z-direction unit vector of collisional system diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index d7542d854..b5ed877c3 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -28,6 +28,7 @@ module function symba_collision_casedisruption(system, param, t) result(status) integer(I4B) :: i, ibiggest, nfrag logical :: lfailure character(len=STRMAX) :: message + real(DP) :: dpe associate(colliders => system%colliders, fragments => system%fragments) @@ -46,6 +47,10 @@ module function symba_collision_casedisruption(system, param, t) result(status) ! Generate the position and velocity distributions of the fragments call fragments%generate_fragments(colliders, system, param, lfailure) + dpe = fragments%pe_after - fragments%pe_before + system%Ecollisions = system%Ecollisions - dpe + system%Euntracked = system%Euntracked + dpe + if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, "No fragment solution found, so treat as a pure hit-and-run") status = ACTIVE @@ -98,6 +103,7 @@ module function symba_collision_casehitandrun(system, param, t) result(status) integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj logical :: lpure character(len=STRMAX) :: message + real(DP) :: dpe associate(colliders => system%colliders, fragments => system%fragments) message = "Hit and run between" @@ -123,6 +129,10 @@ module function symba_collision_casehitandrun(system, param, t) result(status) ! Generate the position and velocity distributions of the fragments call fragments%generate_fragments(colliders, system, param, lpure) + dpe = fragments%pe_after - fragments%pe_before + system%Ecollisions = system%Ecollisions - dpe + system%Euntracked = system%Euntracked + dpe + if (lpure) then call io_log_one_message(FRAGGLE_LOG_OUT, "Should have been a pure hit and run instead") nfrag = 0 @@ -172,8 +182,8 @@ module function symba_collision_casemerge(system, param, t) result(status) integer(I4B) :: status !! Status flag assigned to this outcome ! Internals integer(I4B) :: i, j, k, ibiggest - real(DP) :: pe real(DP), dimension(NDIM) :: L_spin_new + real(DP) :: dpe character(len=STRMAX) :: message associate(colliders => system%colliders, fragments => system%fragments) @@ -207,9 +217,9 @@ module function symba_collision_casemerge(system, param, t) result(status) ! Keep track of the component of potential energy due to the pre-impact colliders%idx for book-keeping ! Get the energy of the system after the collision call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.false.) - pe = fragments%pe_after - fragments%pe_before - system%Ecollisions = system%Ecollisions - pe - system%Euntracked = system%Euntracked + pe + dpe = fragments%pe_after - fragments%pe_before + system%Ecollisions = system%Ecollisions - dpe + system%Euntracked = system%Euntracked + dpe ! Update any encounter lists that have the removed bodies in them so that they instead point to the new