From 639b53ee93b0827f22a12cc9d67e7a5ffab10fc2 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 14 Jun 2021 14:04:13 -0400 Subject: [PATCH] Switched failed fragmentation events to hit-and-run instead of merge --- src/modules/module_interfaces.f90 | 4 +- src/symba/symba_casedisruption.f90 | 11 ++-- src/symba/symba_casesupercatastrophic.f90 | 11 ++-- src/symba/symba_frag_pos.f90 | 62 ++++++++++++----------- 4 files changed, 46 insertions(+), 42 deletions(-) diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index db2d8fa21..de9ffac8d 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -923,7 +923,7 @@ END SUBROUTINE symba_energy INTERFACE subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lmerge) + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) use swiftest_globals USE swiftest_data_structures USE module_symba @@ -938,7 +938,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi real(DP), dimension(:,:), allocatable, intent(inout) :: Ip_frag real(DP), dimension(:,:), allocatable, intent(inout) :: xb_frag, vb_frag, rot_frag real(DP), intent(inout) :: Qloss - logical, intent(out) :: lmerge + logical, intent(out) :: lfailure end subroutine symba_frag_pos END INTERFACE diff --git a/src/symba/symba_casedisruption.f90 b/src/symba/symba_casedisruption.f90 index bd7107953..0a16b304d 100644 --- a/src/symba/symba_casedisruption.f90 +++ b/src/symba/symba_casedisruption.f90 @@ -29,7 +29,7 @@ function symba_casedisruption (symba_plA, family, nmergeadd, mergeadd_list, x, v real(DP), dimension(2) :: vol real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag real(DP), dimension(:), allocatable :: m_frag, rad_frag - logical :: lmerge + logical :: lfailure ! Collisional fragments will be uniformly distributed around the pre-impact barycenter nfrag = 12 ! This value is set for testing. This needs to be updated such that it is calculated or set by the user @@ -72,11 +72,12 @@ function symba_casedisruption (symba_plA, family, nmergeadd, mergeadd_list, x, v ! Put the fragments on the circle surrounding the center of mass of the system call symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lmerge) + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) - if (lmerge) then - write(*,*) 'Should have been a merge instead.' - status = symba_casemerge (symba_plA, family, nmergeadd, mergeadd_list, x, v, mass, radius, L_spin, Ip, param) + if (lfailure) then + write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' + status = ACTIVE + nfrag = 0 else ! Populate the list of new bodies write(*,'("Generating ",I2.0," fragments")') nfrag diff --git a/src/symba/symba_casesupercatastrophic.f90 b/src/symba/symba_casesupercatastrophic.f90 index 04d242ff0..98535e5b2 100644 --- a/src/symba/symba_casesupercatastrophic.f90 +++ b/src/symba/symba_casesupercatastrophic.f90 @@ -29,7 +29,7 @@ function symba_casesupercatastrophic (symba_plA, family, nmergeadd, mergeadd_lis real(DP), dimension(NDIM) :: Ip_new real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag real(DP), dimension(:), allocatable :: m_frag, rad_frag - logical :: lmerge + logical :: lfailure ! Collisional fragments will be uniformly distributed around the pre-impact barycenter nfrag = 20 ! This value is set for testing. This needs to be updated such that it is calculated or set by the user @@ -68,11 +68,12 @@ function symba_casesupercatastrophic (symba_plA, family, nmergeadd, mergeadd_lis ! Put the fragments on the circle surrounding the center of mass of the system call symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lmerge) + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) - if (lmerge) then - write(*,*) 'Should have been a merge instead.' - status = symba_casemerge(symba_plA, family, nmergeadd, mergeadd_list, x, v, mass, radius, L_spin, Ip, param) + if (lfailure) then + write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' + status = ACTIVE + nfrag = 0 else write(*,'("Generating ",I2.0," fragments")') nfrag status = SUPERCATASTROPHIC diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 98f53c146..c8aaa1452 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -1,5 +1,5 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lmerge) + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Places the collision fragments on a circle oriented with a plane defined @@ -23,7 +23,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi real(DP), dimension(:), allocatable, intent(inout) :: m_frag, rad_frag real(DP), dimension(:,:), allocatable, intent(inout) :: Ip_frag real(DP), dimension(:,:), allocatable, intent(inout) :: xb_frag, vb_frag, rot_frag - logical, intent(out) :: lmerge ! Answers the question: Should this have been a merger instead? + logical, intent(out) :: lfailure ! Answers the question: Should this have been a merger instead? real(DP), intent(inout) :: Qloss ! Internals real(DP) :: mscale, rscale, vscale, tscale, Lscale, Escale ! Scale factors that reduce quantities to O(~1) in the collisional system @@ -53,7 +53,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi if (nfrag < NFRAG_MIN) then write(*,*) "symba_frag_pos needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given." - lmerge = .true. + lfailure = .true. return end if @@ -87,50 +87,50 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi r_max_start = norm2(x(:,2) - x(:,1)) try = 1 - lmerge = .false. + lfailure = .false. ke_avg_deficit = 0.0_DP do while (try < MAXTRY) - lmerge = .false. + lfailure = .false. ke_avg_deficit_old = ke_avg_deficit ke_avg_deficit = 0.0_DP subtry = 1 do call set_fragment_position_vectors() - call set_fragment_tangential_velocities(lmerge) + call set_fragment_tangential_velocities(lfailure) ke_avg_deficit = ke_avg_deficit - ke_radial subtry = subtry + 1 - if (.not.lmerge .or. subtry == TANTRY) exit + if (.not.lfailure .or. subtry == TANTRY) exit write(*,*) 'Trying new arrangement' end do ke_avg_deficit = ke_avg_deficit / subtry - if (lmerge) write(*,*) 'Failed to find tangential velocities' + if (lfailure) write(*,*) 'Failed to find tangential velocities' - if (.not.lmerge) then - call set_fragment_radial_velocities(lmerge) - if (lmerge) write(*,*) 'Failed to find radial velocities' - if (.not.lmerge) then + if (.not.lfailure) then + call set_fragment_radial_velocities(lfailure) + if (lfailure) write(*,*) 'Failed to find radial velocities' + if (.not.lfailure) then call calculate_system_energy(linclude_fragments=.true.) write(*,*) 'Qloss : ',Qloss write(*,*) '-dEtot: ',-dEtot write(*,*) 'delta : ',abs((dEtot + Qloss)) if ((abs(dEtot + Qloss) > Etol) .or. (dEtot > 0.0_DP)) then write(*,*) 'Failed due to high energy error: ' - lmerge = .true. + lfailure = .true. else if (abs(dLmag) / Lmag_before > Ltol) then write(*,*) 'Failed due to high angular momentum error: ', dLmag / Lmag_before - lmerge = .true. + lfailure = .true. end if end if end if - if (.not.lmerge) exit + if (.not.lfailure) exit call restructure_failed_fragments() try = try + 1 end do write(*, "(' -------------------------------------------------------------------------------------')") write(*, "(' Final diagnostic')") write(*, "(' -------------------------------------------------------------------------------------')") - if (lmerge) then + if (lfailure) then write(*,*) "symba_frag_pos failed after: ",try," tries" do ii = 1, nfrag vb_frag(:, ii) = vcom(:) @@ -206,29 +206,32 @@ subroutine restore_scale_factors() integer(I4B) :: i call ieee_set_halting_mode(IEEE_ALL,.false.) - mtot = mscale + ! Restore scale factors + xcom(:) = xcom(:) * rscale + vcom(:) = vcom(:) * vscale + + mtot = mtot * mscale mass = mass * mscale radius = radius * rscale x = x * rscale v = v * vscale L_spin = L_spin * Lscale - rot = rot / tscale - - ! Avoid FP exceptions - x_frag(:,1:nfrag) = x_frag(:,1:nfrag) * rscale - v_frag(:,1:nfrag) = v_frag(:,1:nfrag) * vscale - m_frag(1:nfrag) = m_frag(1:nfrag) * mscale - rot_frag(:,1:nfrag) = rot_frag(:,1:nfrag) / tscale - rad_frag(1:nfrag) = rad_frag(1:nfrag) * rscale - ! Convert the final result to the system units - xcom(:) = xcom(:) * rscale - vcom(:) = vcom(:) * vscale + do i = 1, 2 + rot(:,i) = L_spin(:,i) * (mass(i) * radius(i)**2 * Ip(3, i)) + end do + + m_frag = m_frag * mscale + rad_frag = rad_frag * rscale + rot_frag = rot_frag / tscale + x_frag = x_frag * rscale + v_frag = v_frag * vscale + Qloss = Qloss * Escale + do i = 1, nfrag xb_frag(:, i) = x_frag(:, i) + xcom(:) vb_frag(:, i) = v_frag(:, i) + vcom(:) end do - Qloss = Qloss * Escale Etot_before = Etot_before * Escale pe_before = pe_before * Escale ke_spin_before = ke_spin_before * Escale @@ -242,7 +245,6 @@ subroutine restore_scale_factors() Ltot_after = Ltot_after * Lscale Lmag_after = Lmag_after * Lscale - mscale = 1.0_DP rscale = 1.0_DP vscale = 1.0_DP