Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Fixed bad indexing on potential energy calculation. Set mass scaling …
Browse files Browse the repository at this point in the history 
…to be total mass in fragmentation
  • Loading branch information
@daminton
daminton committed Aug 8, 2021
1 parent fb38873 commit 6b04d41
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 27 deletions.
50 changes: 26 additions & 24 deletions src/fragmentation/fragmentation.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,7 +168,7 @@ subroutine set_scale_factors()
! Set scale factors
Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2)))
dscale = sum(radius(:))
mscale = sqrt(Escale * dscale / param%GU)
mscale = mtot
vscale = sqrt(Escale / mscale)
tscale = dscale / vscale
Lscale = mscale * dscale * vscale
Expand Down Expand Up @@ -338,6 +338,7 @@ subroutine calculate_system_energy(linclude_fragments)
call tmpsys%tp%setup(0, param)
deallocate(tmpsys%cb)
allocate(tmpsys%cb, source=cb)
if (allocated(tmpparam)) deallocate(tmpparam)
allocate(tmpparam, source=param)
allocate(ltmp(npl_new))
ltmp(:) = .false.
Expand All @@ -347,13 +348,14 @@ subroutine calculate_system_energy(linclude_fragments)
call tmpsys%rescale(tmpparam, mscale, dscale, tscale)

if (linclude_fragments) then ! Append the fragments if they are included
! Energy calculation requires the fragments to be in the system barcyentric frame, s
! Energy calculation requires the fragments to be in the system barcyentric frame

tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag)
tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * param%GU * dscale**3 * mscale / tscale**2
tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * tmpparam%GU
tmpsys%pl%radius(npl+1:npl_new) = rad_frag(1:nfrag)
tmpsys%pl%xb(:,npl+1:npl_new) = xb_frag(:,1:nfrag)
tmpsys%pl%vb(:,npl+1:npl_new) = vb_frag(:,1:nfrag)
tmpsys%pl%status(npl+1:npl_new) = COLLISION
tmpsys%pl%status(npl+1:npl_new) = ACTIVE
if (param%lrotation) then
tmpsys%pl%Ip(:,npl+1:npl_new) = Ip_frag(:,1:nfrag)
tmpsys%pl%rot(:,npl+1:npl_new) = rot_frag(:,1:nfrag)
Expand Down Expand Up @@ -529,21 +531,21 @@ subroutine set_fragment_tan_vel(lerr)

call calculate_system_energy(linclude_fragments=.true.)
ke_frag_budget = -dEtot - Qloss
write(*,*) '***************************************************'
write(*,*) 'Original dis : ',norm2(x(:,2) - x(:,1))
write(*,*) 'r_max : ',r_max
write(*,*) 'f_spin : ',f_spin
write(*,*) '***************************************************'
write(*,*) 'Energy balance so far: '
write(*,*) 'ke_frag_budget : ',ke_frag_budget
write(*,*) 'ke_orbit_before: ',ke_orbit_before
write(*,*) 'ke_orbit_after : ',ke_orbit_after
write(*,*) 'ke_spin_before : ',ke_spin_before
write(*,*) 'ke_spin_after : ',ke_spin_after
write(*,*) 'pe_before : ',pe_before
write(*,*) 'pe_after : ',pe_after
write(*,*) 'Qloss : ',Qloss
write(*,*) '***************************************************'
!write(*,*) '***************************************************'
!write(*,*) 'Original dis : ',norm2(x(:,2) - x(:,1))
!write(*,*) 'r_max : ',r_max
!write(*,*) 'f_spin : ',f_spin
!write(*,*) '***************************************************'
!write(*,*) 'Energy balance so far: '
!write(*,*) 'ke_frag_budget : ',ke_frag_budget
!write(*,*) 'ke_orbit_before: ',ke_orbit_before
!write(*,*) 'ke_orbit_after : ',ke_orbit_after
!write(*,*) 'ke_spin_before : ',ke_spin_before
!write(*,*) 'ke_spin_after : ',ke_spin_after
!write(*,*) 'pe_before : ',pe_before
!write(*,*) 'pe_after : ',pe_after
!write(*,*) 'Qloss : ',Qloss
!write(*,*) '***************************************************'
if (ke_frag_budget < 0.0_DP) then
write(*,*) 'Negative ke_frag_budget: ',ke_frag_budget
r_max_start = r_max_start / 2
Expand Down Expand Up @@ -608,11 +610,11 @@ subroutine set_fragment_tan_vel(lerr)

! If we are over the energy budget, flag this as a failure so we can try again
lerr = (ke_radial < 0.0_DP)
write(*,*) 'Tangential'
write(*,*) 'ke_frag_budget: ',ke_frag_budget
write(*,*) 'ke_frag_orbit : ',ke_frag_orbit
write(*,*) 'ke_frag_spin : ',ke_frag_spin
write(*,*) 'ke_radial : ',ke_radial
!write(*,*) 'Tangential'
!write(*,*) 'ke_frag_budget: ',ke_frag_budget
!write(*,*) 'ke_frag_orbit : ',ke_frag_orbit
!write(*,*) 'ke_frag_spin : ',ke_frag_spin
!write(*,*) 'ke_radial : ',ke_radial

return

Expand Down
6 changes: 3 additions & 3 deletions src/util/util_get_energy_momentum.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,23 +75,23 @@ module subroutine util_get_energy_momentum_system(self, param)
! Do the central body potential energy component first
!$omp simd
do i = 1, npl
associate(px => pl%xh(1,i), py => pl%xh(2,i), pz => pl%xh(3,i))
associate(px => pl%xb(1,i), py => pl%xb(2,i), pz => pl%xb(3,i))
pecb(i) = -cb%Gmass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2)
end associate
end do

! Do the potential energy between pairs of massive bodies
do k = 1, pl%nplpl
associate(ik => pl%k_plpl(1, k), jk => pl%k_plpl(2, k))
pepl(k) = -pl%Gmass(ik) * pl%mass(jk) / norm2(pl%xh(:, jk) - pl%xh(:, ik))
pepl(k) = -pl%Gmass(ik) * pl%mass(jk) / norm2(pl%xb(:, jk) - pl%xb(:, ik))
lstatpl(k) = (lstatus(ik) .and. lstatus(jk))
end associate
end do

system%ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:))
if (param%lrotation) system%ke_spin = 0.5_DP * sum(kespinpl(1:npl), lstatus(:))

system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(2:npl), lstatus(2:npl))
system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), lstatus(1:npl))

! Potential energy from the oblateness term
if (param%loblatecb) then
Expand Down

0 comments on commit 6b04d41

Please sign in to comment.