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

Commit

Permalink
Merge branch 'OOPHelio' into OOPrestructure
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Jul 6, 2021
2 parents 0c30e50 + 27dbf18 commit 82cb999
Show file tree
Hide file tree
Showing 32 changed files with 1,577 additions and 1,714 deletions.
18 changes: 12 additions & 6 deletions Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@ lib:
ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
ln -s $(SWIFTEST_HOME)/Makefile .; \
make libdir
cd $(SWIFTEST_HOME)/src/driftkick; \
cd $(SWIFTEST_HOME)/src/drift; \
rm -f Makefile.Defines Makefile; \
ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
ln -s $(SWIFTEST_HOME)/Makefile .; \
Expand All @@ -106,6 +106,11 @@ lib:
ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
ln -s $(SWIFTEST_HOME)/Makefile .; \
make libdir
cd $(SWIFTEST_HOME)/src/kick; \
rm -f Makefile.Defines Makefile; \
ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
ln -s $(SWIFTEST_HOME)/Makefile .; \
make libdir
cd $(SWIFTEST_HOME)/src/obl; \
rm -f Makefile.Defines Makefile; \
ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
Expand Down Expand Up @@ -179,20 +184,21 @@ bin: *.f90
clean:
cd $(SWIFTEST_HOME)/src/modules; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/discard; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/driftkick; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/drift; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/eucl; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/gr; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/helio; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/io; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/kick; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/main; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/obl; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/operators; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/orbel; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/rmvs; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/setup; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/util; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/user; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/main; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/util; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/whm; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/rmvs; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/helio; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/bin; rm -f swiftest_*
cd $(SWIFTEST_HOME)/bin; rm -f tool_*
cd $(SWIFTEST_HOME)/lib; rm -f lib*
Expand Down
4 changes: 2 additions & 2 deletions Makefile.Defines
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,8 @@ GPAR = -fopenmp -ftree-parallelize-loops=4
GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak
GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporaries

#FFLAGS = $(IDEBUG) $(HEAPARR)
FFLAGS = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(SIMDVEC) $(PAR)
FFLAGS = $(IDEBUG) $(HEAPARR)
#FFLAGS = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(SIMDVEC) $(PAR)
FORTRAN = ifort
#AR = xiar

Expand Down
4 changes: 2 additions & 2 deletions examples/rmvs_swifter_comparison/9pl_18tp_encounters/param.swifter.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,8 @@ BIN_OUT bin.swifter.dat
OUT_TYPE REAL8
OUT_FORM XV
OUT_STAT UNKNOWN
J2 0.0
J4 0.0
J2 4.7535806948127355e-12
J4 -2.2473967953572827e-18
CHK_CLOSE yes
CHK_RMIN 0.004650467260962157
CHK_RMAX 1000.0
Expand Down
219 changes: 84 additions & 135 deletions examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb
View file

Large diffs are not rendered by default.

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 !! Input collection of parameters 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 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 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
0 src/driftkick/drift.f90 → src/drift/drift.f90
View file
File renamed without changes.
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 !! Input collection 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 !! Input collection 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 82cb999

Please sign in to comment.