Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Hyperbolic and parabolic setups in set_binary #443

Merged
merged 5 commits into from
Jul 3, 2023
Merged

Hyperbolic and parabolic setups in set_binary #443

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
c33998b
(step) add fatal error if all particles accreted
danieljprice Jun 21, 2023
5612bb9
(timing) do not print blank timers
danieljprice Jul 3, 2023
2f7c71b
(set_flyby) typo in comment fixed
danieljprice Jul 3, 2023
16bbc6f
(set_binary; #302) handle parabolic and hyperbolic orbit setups
danieljprice Jul 3, 2023
7b79829
Merge branch 'master' into hyperbolic
danieljprice Jul 3, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 26 additions & 11 deletions src/main/step_leapfrog.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -753,7 +753,7 @@ end subroutine step_extern_sph_gr

subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,metrics,metricderivs,fext,time)
use dim, only:maxptmass,maxp,maxvxyzu
use io, only:iverbose,id,master,iprint,warning
use io, only:iverbose,id,master,iprint,warning,fatal
use externalforces, only:externalforce,accrete_particles,update_externalforce
use options, only:iexternalforce
use part, only:maxphase,isdead_or_accreted,iamboundary,igas,iphase,iamtype,&
Expand All @@ -769,7 +769,7 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me
real, intent(in) :: dtsph,time
real, intent(inout) :: dtextforce
real, intent(inout) :: xyzh(:,:),vxyzu(:,:),fext(:,:),pxyzu(:,:),dens(:),metrics(:,:,:,:),metricderivs(:,:,:,:)
integer :: i,itype,nsubsteps,naccreted,its,ierr
integer :: i,itype,nsubsteps,naccreted,its,ierr,nlive
real :: timei,t_end_step,hdt,pmassi
real :: dt,dtf,dtextforcenew,dtextforce_min
real :: pri,spsoundi,pondensi,tempi,gammai
Expand Down Expand Up @@ -866,7 +866,8 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me
tempi = eos_vars(itemp,i)
rhoi = rhoh(hi,massoftype(igas))

! Note: grforce needs derivatives of the metric, which do not change between pmom iterations
! Note: grforce needs derivatives of the metric,
! which do not change between pmom iterations
pmom_iterations: do while (its <= itsmax .and. .not. converged)
its = its + 1
pprev = pxyz
Expand All @@ -879,7 +880,8 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me
if (pmom_err < ptol) converged = .true.
fexti = fstar
enddo pmom_iterations
if (its > itsmax ) call warning('step_extern_gr','Reached max number of pmom iterations. pmom_err ',val=pmom_err)
if (its > itsmax ) call warning('step_extern_gr',&
'max # of pmom iterations',var='pmom_err',val=pmom_err)
pitsmax = max(its,pitsmax)
perrmax = max(pmom_err,perrmax)

Expand All @@ -892,9 +894,11 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me
its = 0
converged = .false.
vxyz_star = vxyz
! Note: since particle positions change between iterations the metric and its derivatives need to be updated.
! cons2prim does not require derivatives of the metric, so those can updated once the iterations
! are complete, in order to reduce the number of computations.
! Note: since particle positions change between iterations
! the metric and its derivatives need to be updated.
! cons2prim does not require derivatives of the metric,
! so those can updated once the iterations are complete
! in order to reduce the number of computations.
xyz_iterations: do while (its <= itsmax .and. .not. converged)
its = its+1
xyz_prev = xyz
Expand Down Expand Up @@ -942,6 +946,7 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me
!
accretedmass = 0.
naccreted = 0
nlive = 0
dtextforce_min = bignumber

!$omp parallel default(none) &
Expand All @@ -952,7 +957,7 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me
!$omp private(pri,pondensi,spsoundi,tempi,dtf) &
!$omp firstprivate(itype,pmassi) &
!$omp reduction(min:dtextforce_min) &
!$omp reduction(+:accretedmass,naccreted) &
!$omp reduction(+:accretedmass,naccreted,nlive) &
!$omp shared(idamp,damp_fac)
!$omp do
accreteloop: do i=1,npart
Expand Down Expand Up @@ -986,11 +991,16 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me
naccreted = naccreted + 1
endif
endif
nlive = nlive + 1
endif
enddo accreteloop
!$omp enddo
!$omp end parallel

if (npart > 2 .and. nlive < 2) then
call fatal('step','all particles accreted',var='nlive',ival=nlive)
endif

if (iverbose >= 2 .and. id==master .and. naccreted /= 0) write(iprint,"(a,es10.3,a,i4,a)") &
'Step: at time ',timei,', ',naccreted,' particles were accreted. Mass accreted = ',accretedmass

Expand Down Expand Up @@ -1096,7 +1106,7 @@ subroutine step_extern(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,fext,fxyzu,time,
real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:),fxyz_ptmass(:,:)
integer(kind=1), intent(in) :: nbinmax
integer(kind=1), intent(inout) :: ibin_wake(:)
integer :: i,itype,nsubsteps,naccreted,nfail,nfaili,merge_n
integer :: i,itype,nsubsteps,naccreted,nfail,nfaili,merge_n,nlive
integer :: merge_ij(nptmass)
integer(kind=1) :: ibin_wakei
real :: timei,hdt,fextx,fexty,fextz,fextxi,fextyi,fextzi,phii,pmassi
Expand Down Expand Up @@ -1375,6 +1385,7 @@ subroutine step_extern(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,fext,fxyzu,time,
accretedmass = 0.
nfail = 0
naccreted = 0
nlive = 0
ibin_wakei = 0
dptmass(:,:) = 0.

Expand All @@ -1387,7 +1398,7 @@ subroutine step_extern(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,fext,fxyzu,time,
!$omp reduction(+:dptmass) &
!$omp private(i,accreted,nfaili,fxi,fyi,fzi) &
!$omp firstprivate(itype,pmassi,ibin_wakei) &
!$omp reduction(+:accretedmass,nfail,naccreted)
!$omp reduction(+:accretedmass,nfail,naccreted,nlive)
!$omp do
accreteloop: do i=1,npart
if (.not.isdead_or_accreted(xyzh(4,i))) then
Expand Down Expand Up @@ -1430,12 +1441,16 @@ subroutine step_extern(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,fext,fxyzu,time,
endif
if (nfaili > 1) nfail = nfail + 1
endif
nlive = nlive + 1
endif

enddo accreteloop
!$omp enddo
!$omp end parallel

if (npart > 2 .and. nlive < 2) then
call fatal('step','all particles accreted',var='nlive',ival=nlive)
endif

!
! reduction of sink particle changes across MPI
!
Expand Down
30 changes: 28 additions & 2 deletions src/main/utils_binary.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,13 @@ real function get_E_from_mean_anomaly(M_ref,ecc) result(E)
M_guess = M_ref - 2.*tol

do while (abs(M_ref - M_guess) > tol)
M_guess = E_guess - ecc*sin(E_guess)
if (ecc < 1.) then ! eccentric
M_guess = E_guess - ecc*sin(E_guess)
elseif (ecc > 1.) then ! hyperbolic
M_guess = ecc*sinh(E_guess) - E_guess
else ! parabolic
M_guess = E_guess + 1./3.*E_guess**3
endif
if (M_guess > M_ref) then
E_right = E_guess
else
Expand All @@ -75,13 +81,33 @@ real function get_E_from_mean_anomaly(M_ref,ecc) result(E)

end function get_E_from_mean_anomaly

!---------------------------------------------------------------
!+
! Get eccentric (or parabolic/hyperbolic) anomaly from true anomaly
! https://space.stackexchange.com/questions/23128/design-of-an-elliptical-transfer-orbit/23130#23130
!+
!---------------------------------------------------------------
real function get_E_from_true_anomaly(theta,ecc) result(E)
real, intent(in) :: theta ! true anomaly in radians
real, intent(in) :: ecc ! eccentricity

if (ecc < 1.) then
E = atan2(sqrt(1. - ecc**2)*sin(theta),(ecc + cos(theta)))
elseif (ecc > 1.) then ! hyperbolic
!E = atanh(sqrt(ecc**2 - 1.)*sin(theta)/(ecc + cos(theta)))
E = 2.*atanh(sqrt((ecc - 1.)/(ecc + 1.))*tan(0.5*theta))
else ! parabolic
E = tan(0.5*theta)
endif

end function get_E_from_true_anomaly

!-----------------------------------------------------------------------
!+
! Calculate semi-major axis, ecc, ra and rp from radius(3), velocity(3)
! mass of central object and iexternalforce (for LT corrections)
!+
!-----------------------------------------------------------------------

subroutine get_orbit_bits(vel,rad,m1,iexternalforce,semia,ecc,ra,rp)
real, intent(in) :: m1, vel(3), rad(3)
integer, intent(in) :: iexternalforce
Expand Down
7 changes: 2 additions & 5 deletions src/main/utils_timing.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -293,11 +293,10 @@ subroutine print_timer(lu,itimer,time_total)
integer, intent(in) :: itimer
real(kind=4), intent(in) :: time_total

! Print label and tree structure
write(lu, "(1x,a)", advance="no") trim(timers(itimer)%treelabel)

! Print timings
if (timers(itimer)%wall > epsilon(0._4)) then
! Print label and tree structure
write(lu, "(1x,a)", advance="no") trim(timers(itimer)%treelabel)
if (time_total > 7200.0) then
write(lu,"(' ',f7.2,'h ',f8.2,'h ',f6.2,' ',f6.2,'%',' ',f6.2,'%')") &
timers(itimer)%wall/3600.,&
Expand All @@ -320,8 +319,6 @@ subroutine print_timer(lu,itimer,time_total)
timers(itimer)%loadbal*100.,&
timers(itimer)%wall/time_total*100.
endif
else
write(lu, "()")
endif

end subroutine print_timer
Expand Down
Loading