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

Commit

Permalink
Major restructuring of classes and method interfaces. Trying to conso…
Browse files Browse the repository at this point in the history 
…lidate all cb, pl, tp arguments into system class
  • Loading branch information
@daminton
daminton committed Jul 6, 2021
1 parent 058889a commit 69acd29
Show file tree
Hide file tree
Showing 21 changed files with 1,120 additions and 1,214 deletions.
60 changes: 26 additions & 34 deletions src/discard/discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,37 +9,35 @@ module subroutine discard_system(self, param)
!! Adapted from David E. Kaufmann's Swifter routine: discard.f90
!! Adapted from Hal Levison's Swift routine discard.f
implicit none
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters

if (self%tp%nbody == 0) return
select type(self)
select type(system => self)
class is (whm_nbody_system)
associate(cb => self%cb, pl => self%pl, tp => self%tp, t => param%t, dt => param%dt, &
msys => self%msys, discards => self%tp_discards, &
ntp => self%tp%nbody, ldiscard => self%tp%ldiscard)
associate(cb => system%cb, pl => system%pl, npl => system%pl%nbody, tp => system%tp, ntp => system%tp%nbody)
if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. &
(param%rmaxu >= 0.0_DP) .or. ((param%qmin >= 0.0_DP) .and. (param%qmin_coord == "BARY"))) then
call pl%h2b(cb)
if (npl > 0) call pl%h2b(cb)
if (ntp > 0) call tp%h2b(cb)
end if
if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) then
if (ntp > 0) call tp%discard_sun(cb, param, t, msys)
if (ntp > 0) call tp%discard_sun(system, param)
end if
if (param%qmin >= 0.0_DP .and. ntp > 0) call tp%discard_peri(cb, pl, param, t, msys)
if (param%lclose .and. ntp > 0) call tp%discard_pl(pl, t, dt)
if (param%qmin >= 0.0_DP .and. ntp > 0) call tp%discard_peri(system, param)
if (param%lclose .and. ntp > 0) call tp%discard_pl(system, param)

if (any(tp%ldiscard(1:ntp))) then
! Spill the discards to the spill list
call tp%spill(discards, ldiscard)
call self%write_discard(param, discards)
call tp%spill(system%tp_discards, tp%ldiscard)
call self%write_discard(param, system%tp_discards)
end if
end associate
end select
return
end subroutine discard_system

module subroutine discard_sun_tp(self, cb, param, t, msys)
module subroutine discard_sun_tp(self, system, param)
!! author: David A. Minton
!!
!! Check to see if test particles should be discarded based on their positions relative to the Sun
Expand All @@ -49,16 +47,14 @@ module subroutine discard_sun_tp(self, cb, param, t, msys)
!! Adapted from Hal Levison's Swift routine discard_sun.f
implicit none
! Arguments
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object
class(swiftest_parameters), intent(in) :: param !! parameters
real(DP), intent(in) :: t !! Current simulation tim
real(DP), intent(in) :: msys !! Total system mass
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: i
real(DP) :: energy, vb2, rb2, rh2, rmin2, rmax2, rmaxu2

associate(tp => self, ntp => self%nbody)
associate(tp => self, ntp => self%nbody, cb => system%cb, t => param%t, msys => system%msys)
rmin2 = max(param%rmin * param%rmin, cb%radius * cb%radius)
rmax2 = param%rmax**2
rmaxu2 = param%rmaxu**2
Expand Down Expand Up @@ -90,7 +86,7 @@ module subroutine discard_sun_tp(self, cb, param, t, msys)
return
end subroutine discard_sun_tp

module subroutine discard_peri_tp(self, cb, pl, param, t, msys)
module subroutine discard_peri_tp(self, system, param)
!! author: David A. Minton
!!
!! Check to see if a test particle should be discarded because its perihelion distance becomes too small
Expand All @@ -99,19 +95,16 @@ module subroutine discard_peri_tp(self, cb, pl, param, t, msys)
!! Adapted from Hal Levison's Swift routine discard_peri.f
implicit none
! Arguments
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object
class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object
class(swiftest_parameters), intent(in) :: param !! parameters
real(DP), intent(in) :: t !! Current simulation tim
real(DP), intent(in) :: msys !! Total system mass
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameterss
! Internals
logical, save :: lfirst = .true.
integer(I4B) :: i, j, ih
real(DP) :: r2
real(DP), dimension(NDIM) :: dx

associate(tp => self, ntp => self%nbody, npl => pl%nbody, qmin_coord => param%qmin_coord)
associate(cb => system%cb, tp => self, ntp => self%nbody, pl => system%pl, npl => system%pl%nbody, qmin_coord => param%qmin_coord, t => param%t, msys => system%msys)
if (lfirst) then
call util_hills(npl, pl)
call util_peri(lfirst, ntp, tp, cb%Gmass, msys, param%qmin_coord)
Expand Down Expand Up @@ -145,7 +138,7 @@ module subroutine discard_peri_tp(self, cb, pl, param, t, msys)

end subroutine discard_peri_tp

module subroutine discard_pl_tp(self, pl, t, dt)
module subroutine discard_pl_tp(self, system, param)
!! author: David A. Minton
!!
!! Check to see if test particles should be discarded based on their positions relative to the massive bodies
Expand All @@ -154,16 +147,15 @@ module subroutine discard_pl_tp(self, pl, t, dt)
!! Adapted from Hal Levison's Swift routine discard_pl.f
implicit none
! Arguments
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object
real(DP), intent(in) :: t !! Current simulation tim
real(DP), intent(in) :: dt !! Stepsize
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: i, j, isp
real(DP) :: r2min, radius
real(DP), dimension(NDIM) :: dx, dv

associate(tp => self, ntp => self%nbody, npl => pl%nbody)
associate(tp => self, ntp => self%nbody, pl => system%pl, npl => system%pl%nbody, t => param%t, dt => param%dt)
do i = 1, ntp
if (tp%status(i) == ACTIVE) then
do j = 1, npl
Expand All @@ -187,7 +179,7 @@ module subroutine discard_pl_tp(self, pl, t, dt)

end subroutine discard_pl_tp

module subroutine discard_pl_close(dx, dv, dt, r2crit, iflag, r2min)
subroutine discard_pl_close(dx, dv, dt, r2crit, iflag, r2min)
!! author: David A. Minton
!!
!! Check to see if a test particle and massive body are having, or will have within the next time step, an encounter such
Expand Down
14 changes: 10 additions & 4 deletions src/gr/gr.f90
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,22 @@
submodule(swiftest_classes) s_gr
use swiftest
contains
module procedure gr_getaccb_ns_body
subroutine gr_getaccb_ns_body(self, cb, param, agr, agr0)

!! author: David A. Minton
!!
!! Add relativistic correction acceleration for non-symplectic integrators
!! Based on Quinn, Tremaine, & Duncan 1998
!!
!! Adapted from David A. Minton's Swifter routine routine gr_getaccb_ns.f90
implicit none

! Arguments
class(swiftest_body), intent(inout) :: self
class(swiftest_cb), intent(inout) :: cb
class(swiftest_parameters), intent(in) :: param
real(DP), dimension(:, :), intent(inout) :: agr
real(DP), dimension(NDIM), intent(out) :: agr0
! Internals
real(DP), dimension(NDIM) :: xh, vh
real(DP) :: rmag, rdotv, vmag2
integer(I4B) :: i
Expand All @@ -29,7 +36,6 @@
agr0 = 0.0_DP
select type(self)
class is (swiftest_pl)
!do concurrent(i = 1:NDIM)
do i = 1, NDIM
agr0(i) = -sum(self%Gmass(1:n) * agr(1:n, i) / msun)
end do
Expand All @@ -39,6 +45,6 @@

return

end procedure gr_getaccb_ns_body
end subroutine gr_getaccb_ns_body

end submodule s_gr
68 changes: 30 additions & 38 deletions src/helio/helio_drift.f90
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
submodule (helio_classes) s_helio_drift
use swiftest
contains
module subroutine helio_drift_pl(self, cb, param, dt)
module subroutine helio_drift_pl(self, system, param, dt)

!! author: David A. Minton
!!
!! Loop through massive bodies and call Danby drift routine
Expand All @@ -11,22 +12,17 @@ module subroutine helio_drift_pl(self, cb, param, dt)
!! Adapted from Hal Levison's Swift routine drift.f
implicit none
! Arguments
class(helio_pl), intent(inout) :: self !! Helio test particle data structure
class(swiftest_cb), intent(inout) :: cb !! Helio central body particle data structuree
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of
real(DP), intent(in) :: dt !! Stepsize
class(helio_pl), intent(inout) :: self !! Helio massive body object
class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of
real(DP), intent(in) :: dt !! Stepsize)
! Internals
integer(I4B) :: i !! Loop counter
real(DP) :: rmag, vmag2, energy
integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag
real(DP), dimension(:), allocatable :: dtp

associate(npl => self%nbody, &
xh => self%xh, &
vb => self%vb, &
status => self%status,&
mu => self%mu)

associate(pl => self, npl => self%nbody)
if (npl == 0) return

allocate(iflag(npl))
Expand All @@ -35,23 +31,23 @@ module subroutine helio_drift_pl(self, cb, param, dt)

if (param%lgr) then
do i = 1,npl
rmag = norm2(xh(:, i))
vmag2 = dot_product(vb(:, i), vb(:, i))
energy = 0.5_DP * vmag2 - mu(i) / rmag
rmag = norm2(pl%xh(:, i))
vmag2 = dot_product(pl%vb(:, i), pl%vb(:, i))
energy = 0.5_DP * vmag2 - pl%mu(i) / rmag
dtp(i) = dt * (1.0_DP + 3 * param%inv_c2 * energy)
end do
else
dtp(:) = dt
end if

call drift_one(mu(1:npl), xh(1,1:npl), xh(2,1:npl), xh(3,1:npl), &
vb(1,1:npl), vb(2,1:npl), vb(3,1:npl), &
dtp(1:npl), iflag(1:npl))
call drift_one(pl%mu(1:npl), pl%xh(1,1:npl), pl%xh(2,1:npl), pl%xh(3,1:npl), &
pl%vb(1,1:npl), pl%vb(2,1:npl), pl%vb(3,1:npl), &
dtp(1:npl), iflag(1:npl))
if (any(iflag(1:npl) /= 0)) then
do i = 1, npl
write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!"
write(*, *) xh(:,i)
write(*, *) vb(:,i)
write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!"
write(*, *) pl%xh(:,i)
write(*, *) pl%vb(:,i)
write(*, *) " stopping "
call util_exit(FAILURE)
end do
Expand Down Expand Up @@ -94,7 +90,7 @@ module subroutine helio_drift_linear_pl(self, cb, dt, pt)

end subroutine helio_drift_linear_pl

module subroutine helio_drift_tp(self, cb, param, dt)
module subroutine helio_drift_tp(self, system, param, dt)
!! author: David A. Minton
!!
!! Loop through test particles and call Danby drift routine
Expand All @@ -103,21 +99,17 @@ module subroutine helio_drift_tp(self, cb, param, dt)
!! Adapted from Hal Levison's Swift routine drift_tp.f
implicit none
! Arguments
class(helio_tp), intent(inout) :: self !! Helio test particle data structure
class(swiftest_cb), intent(inout) :: cb !! Helio central body particle data structuree
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of
real(DP), intent(in) :: dt !! Stepsize
class(helio_tp), intent(inout) :: self !! Helio test particle object
class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of
real(DP), intent(in) :: dt !! Stepsiz
! Internals
integer(I4B) :: i !! Loop counter
real(DP) :: rmag, vmag2, energy
real(DP), dimension(:), allocatable :: dtp
integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag

associate(ntp => self%nbody, &
xh => self%xh, &
vh => self%vh, &
status => self%status,&
mu => self%mu)
associate(tp => self, ntp => self%nbody)
if (ntp == 0) return
allocate(iflag(ntp))
allocate(dtp(ntp))
Expand All @@ -128,21 +120,21 @@ module subroutine helio_drift_tp(self, cb, param, dt)

if (param%lgr) then
do i = 1,ntp
rmag = norm2(xh(:, i))
vmag2 = dot_product(vh(:, i), vh(:, i))
energy = 0.5_DP * vmag2 - mu(i) / rmag
rmag = norm2(tp%xh(:, i))
vmag2 = dot_product(tp%vh(:, i), tp%vh(:, i))
energy = 0.5_DP * vmag2 - tp%mu(i) / rmag
dtp(i) = dt * (1.0_DP + 3 * param%inv_c2 * energy)
end do
else
dtp(:) = dt
end if
call drift_one(mu(1:ntp), xh(1,1:ntp), xh(2,1:ntp), xh(3,1:ntp), &
vh(1,1:ntp), vh(2,1:ntp), vh(3,1:ntp), &
dtp(1:ntp), iflag(1:ntp))
call drift_one(tp%mu(1:ntp), tp%xh(1,1:ntp), tp%xh(2,1:ntp), tp%xh(3,1:ntp), &
tp%vh(1,1:ntp), tp%vh(2,1:ntp), tp%vh(3,1:ntp), &
dtp(1:ntp), iflag(1:ntp))
if (any(iflag(1:ntp) /= 0)) then
status = DISCARDED_DRIFTERR
tp%status = DISCARDED_DRIFTERR
do i = 1, ntp
if (iflag(i) /= 0) write(*, *) "Particle ", self%name(i), " lost due to error in Danby drift"
if (iflag(i) /= 0) write(*, *) "Particle ", tp%name(i), " lost due to error in Danby drift"
end do
end if
end associate
Expand Down
Loading

0 comments on commit 69acd29

Please sign in to comment.