From d556d8a337b1982ae5b83a0f520bac48b53ecfa7 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 14 Dec 2022 12:44:47 -0500 Subject: [PATCH] Fixed a bunch of issues to get energy values computed in mergers --- src/fraggle/fraggle_set.f90 | 5 +++-- src/modules/fraggle_classes.f90 | 12 ++++++------ src/symba/symba_collision.f90 | 20 +++++++++++--------- 3 files changed, 20 insertions(+), 17 deletions(-) diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index dfa0f64e6..a8e61130d 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -287,7 +287,6 @@ module subroutine fraggle_set_original_scale_factors(self, colliders) do i = 1, 2 colliders%rot(:,i) = colliders%L_spin(:,i) * (colliders%mass(i) * colliders%radius(i)**2 * colliders%Ip(3, i)) end do - frag%Qloss = frag%Qloss * frag%Escale frag%mtot = frag%mtot * frag%mscale frag%mass = frag%mass * frag%mscale @@ -300,7 +299,9 @@ module subroutine fraggle_set_original_scale_factors(self, colliders) frag%rb(:, i) = frag%x_coll(:, i) + frag%rbcom(:) frag%vb(:, i) = frag%v_coll(:, i) + frag%vbcom(:) end do - + + frag%Qloss = frag%Qloss * frag%Escale + frag%Lorbit_before(:) = frag%Lorbit_before * frag%Lscale frag%Lspin_before(:) = frag%Lspin_before * frag%Lscale frag%Ltot_before(:) = frag%Ltot_before * frag%Lscale diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 00b791f71..8f6cd4310 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -86,12 +86,12 @@ module fraggle_classes real(DP) :: Etot_before, Etot_after !! Before/after total system energy ! 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 !! Distance dimension scale factor - real(DP) :: mscale !! Mass scale factor - real(DP) :: tscale !! Time scale factor - real(DP) :: vscale !! Velocity scale factor (a convenience unit that is derived from dscale and tscale) - real(DP) :: Escale !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale) - real(DP) :: Lscale !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale) + 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) contains procedure :: generate_fragments => fraggle_generate_fragments !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. procedure :: accel => fraggle_placeholder_accel !! Placeholder subroutine to fulfill requirement for an accel method diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index d9808f5cd..c6ad292c4 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -185,6 +185,10 @@ module function symba_collision_casemerge(system, param, t) result(status) class is (symba_pl) call fragments%set_mass_dist(colliders, param) + + ! Calculate the initial energy of the system without the collisional family + call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.true.) + ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) fragments%id(1) = pl%id(ibiggest) fragments%rb(:,1) = fragments%rbcom(:) @@ -197,18 +201,16 @@ module function symba_collision_casemerge(system, param, t) result(status) ! Assume prinicpal axis rotation on 3rd Ip axis fragments%rot(:,1) = L_spin_new(:) / (fragments%Ip(3,1) * fragments%mass(1) * fragments%radius(1)**2) else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable - param%Lescape(:) = param%Lescape(:) + colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + system%Lescape(:) = system%Lescape(:) + colliders%L_orbit(:,1) + colliders%L_orbit(:,2) end if ! Keep track of the component of potential energy due to the pre-impact colliders%idx for book-keeping - pe = 0.0_DP - do j = 1, colliders%ncoll - do i = j + 1, colliders%ncoll - pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%rb(:, i) - pl%rb(:, j)) - end do - end do - system%Ecollisions = system%Ecollisions + pe - system%Euntracked = system%Euntracked - pe + ! 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 + ! Update any encounter lists that have the removed bodies in them so that they instead point to the new do k = 1, system%plplenc_list%nenc