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

Commit

Permalink
Fixed scaling issue in radial velocity step of fragmentation initiali…
Browse files Browse the repository at this point in the history 
…zation
  • Loading branch information
@daminton
daminton committed Aug 11, 2021
1 parent 8228d43 commit 46f333e
Show file tree
Hide file tree
Showing 3 changed files with 37 additions and 32 deletions.
63 changes: 33 additions & 30 deletions src/fragmentation/fragmentation.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,24 +133,24 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
call restore_scale_factors()
call calculate_system_energy(linclude_fragments=.true.)

write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' Final diagnostic')")
write(*, "(' -------------------------------------------------------------------------------------')")
if (lfailure) then
write(*,*) "symba_frag_pos failed after: ",try," tries"
do ii = 1, nfrag
vb_frag(:, ii) = vcom(:)
end do
else
write(*,*) "symba_frag_pos succeeded after: ",try," tries"
write(*, "(' dL_tot should be very small' )")
write(*,fmtlabel) ' dL_tot |', dLmag / Lmag_before
write(*, "(' dE_tot should be negative and equal to Qloss' )")
write(*,fmtlabel) ' dE_tot |', dEtot / abs(Etot_before)
write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before)
write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before)
end if
write(*, "(' -------------------------------------------------------------------------------------')")
! write(*, "(' -------------------------------------------------------------------------------------')")
! write(*, "(' Final diagnostic')")
! write(*, "(' -------------------------------------------------------------------------------------')")
! if (lfailure) then
! write(*,*) "symba_frag_pos failed after: ",try," tries"
! do ii = 1, nfrag
! vb_frag(:, ii) = vcom(:)
! end do
! else
! write(*,*) "symba_frag_pos succeeded after: ",try," tries"
! write(*, "(' dL_tot should be very small' )")
! write(*,fmtlabel) ' dL_tot |', dLmag / Lmag_before
! write(*, "(' dE_tot should be negative and equal to Qloss' )")
! write(*,fmtlabel) ' dE_tot |', dEtot / abs(Etot_before)
! write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before)
! write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before)
! end if
! write(*, "(' -------------------------------------------------------------------------------------')")

call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily

Expand Down Expand Up @@ -660,12 +660,12 @@ 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(*,*) 'Failure? ',lerr
write(*,*) 'ke_frag_budget: ',ke_frag_budget
write(*,*) 'ke_frag_spin : ',ke_frag_spin
write(*,*) 'ke_tangential : ',ke_frag_orbit
write(*,*) 'ke_remainder : ',ke_radial
! write(*,*) 'Tangential'
! write(*,*) 'Failure? ',lerr
! write(*,*) 'ke_frag_budget: ',ke_frag_budget
! write(*,*) 'ke_frag_spin : ',ke_frag_spin
! write(*,*) 'ke_tangential : ',ke_frag_orbit
! write(*,*) 'ke_remainder : ',ke_radial

return
end subroutine set_fragment_tan_vel
Expand Down Expand Up @@ -782,6 +782,9 @@ subroutine set_fragment_radial_velocities(lerr)
v_r_mag = util_minimize_bfgs(objective_function, nfrag, v_r_initial, TOL, lerr)
! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin)
vb_frag(:,1:nfrag) = vmag_to_vb(v_r_mag(1:nfrag), v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:))
do i = 1, nfrag
v_frag(:, i) = vb_frag(:, i) - vcom(:)
end do
call add_fragments_to_tmpsys()

do concurrent(i = 1:nfrag)
Expand All @@ -790,12 +793,12 @@ subroutine set_fragment_radial_velocities(lerr)
end do
ke_frag_orbit = 0.5_DP * sum(kearr(:))
ke_frag_spin = 0.5_DP * sum(kespinarr(:))
write(*,*) 'Radial'
write(*,*) 'Failure? ',lerr
write(*,*) 'ke_frag_budget: ',ke_frag_budget
write(*,*) 'ke_frag_spin : ',ke_frag_spin
write(*,*) 'ke_frag_orbit : ',ke_frag_orbit
write(*,*) 'ke_remainder : ',ke_frag_budget - (ke_frag_orbit + ke_frag_spin)
! write(*,*) 'Radial'
! write(*,*) 'Failure? ',lerr
! write(*,*) 'ke_frag_budget: ',ke_frag_budget
! write(*,*) 'ke_frag_spin : ',ke_frag_spin
! write(*,*) 'ke_frag_orbit : ',ke_frag_orbit
! write(*,*) 'ke_remainder : ',ke_frag_budget - (ke_frag_orbit + ke_frag_spin)
lerr = .false.

return
Expand Down
5 changes: 3 additions & 2 deletions src/symba/symba_fragmentation.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -390,11 +390,12 @@ subroutine symba_fragmentation_mergeaddsub(system, param, family, id_frag, Ip_fr
nfrag = size(m_frag(:))
lmask(:) = .false.
lmask(family(:)) = .true.
pl%status(family(:)) = INACTIVE

pl%status(family(:)) = MERGED
nstart = pl_discards%nbody + 1
nend = pl_discards%nbody + nfamily
call pl_discards%append(pl, lmask)
pl%ldiscard(family(:)) = .true.
pl%lcollision(family(:)) = .true.

! Record how many bodies were subtracted in this event
pl_discards%ncomp(nstart:nend) = nfamily
Expand Down
1 change: 1 addition & 0 deletions src/util/util_rescale.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ module subroutine util_rescale_system(self, param, mscale, dscale, tscale)
cb%radius = cb%radius / dscale
cb%xb(:) = cb%xb(:) / dscale
cb%vb(:) = cb%vb(:) / vscale
cb%rot(:) = cb%rot(:) * tscale
pl%mass(1:npl) = pl%mass(1:npl) / mscale
pl%Gmass(1:npl) = param%GU * pl%mass(1:npl)
pl%radius(1:npl) = pl%radius(1:npl) / dscale
Expand Down

0 comments on commit 46f333e

Please sign in to comment.