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

Commit

Permalink
Merge branch 'debug'
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed May 17, 2023
2 parents 9af76f4 + 4ed6bf4 commit d0c8acb
Show file tree
Hide file tree
Showing 15 changed files with 372 additions and 85 deletions.
25 changes: 10 additions & 15 deletions cmake/Modules/SetFortranFlags.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,7 +158,7 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}"

# Strict model for floating-point calculations (precise and except)
SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}"
Fortran "-fp-model strict" # Intel
Fortran "-fp-model=strict" # Intel
)

# Enables floating-point invalid, divide-by-zero, and overflow exceptions
Expand Down Expand Up @@ -273,11 +273,6 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}"
Fortran "-align all" # Intel
)

# Assume all objects are contiguous in memory
SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}"
Fortran "-assume contiguous_assumed_shape" # Intel
)

# Generate an extended set of vector functions
SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}"
Fortran "-vecabi=cmdtarget" # Intel
Expand All @@ -289,17 +284,17 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}"
)

# Generate fused multiply-add instructions
SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}"
Fortran "-fma" # Intel
)
SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}"
Fortran "-fma" # Intel
)

# Generate fused multiply-add instructions
SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}"
Fortran "-qmkl=cluster" # Intel
Fortran "-qmkl" # Intel
Fortran "-mkl" # Old Intel
)

SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}"
Fortran "-qmkl=cluster" # Intel
Fortran "-qmkl" # Intel
Fortran "-mkl" # Old Intel
)
#####################
### MATH FLAGS ###
#####################
Expand Down
1 change: 1 addition & 0 deletions examples/rmvs_swifter_comparison/.gitignore
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
*
!.gitignore
!1pl_1tp_encounter
!8pl_16tp_encounters
5 changes: 5 additions & 0 deletions examples/rmvs_swifter_comparison/8pl_16tp_encounters/.gitignore
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
*
!.gitignore
!init_cond.py
!swiftest_vs_swifter.py
!swiftest_vs_swifter.ipynb
32 changes: 32 additions & 0 deletions examples/rmvs_swifter_comparison/8pl_16tp_encounters/init_cond.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
#!/usr/bin/env python3
import numpy as np
import swiftest

tstart = 0.0
dt = 1.0
tstop = 365.25e2
tstep_out = 100*dt

sim = swiftest.Simulation(simdir="swiftest_sim",init_cond_format="XV",output_format="XV",general_relativity=False, integrator="RMVS", rhill_present=True, MU="Msun", DU="AU", TU="d")
sim.clean()
sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"])

plname = sim.init_cond['name'].where(sim.init_cond['name'] != "Sun", drop=True)
pl = sim.init_cond.sel(name=plname)

for i,n in enumerate(pl.name):
pli = pl.sel(name=n)
rstart = 2 * pli['radius'].data[0] # Start the test particles at a multiple of the planet radius away
vstart = 1.5 * np.sqrt(2 * pli['Gmass'].data[0]) / rstart # Start the test particle velocities at a multiple of the escape speed
rstart_vec = np.array([rstart / np.sqrt(2.0), rstart / np.sqrt(2.0), 0.0])
vstart_vec = np.array([vstart, 0.0, 0.0])
rp = pli['rh'].data[0]
vp = pli['vh'].data[0]
sim.add_body(name=[f"TestParticle{100+i}",f"TestParticle{200+i}"],rh=[rp+rstart_vec, rp-rstart_vec],vh=[vp+vstart_vec, vp-vstart_vec])


sim.set_parameter(tstart=tstart, tstop=tstop, dt=dt, tstep_out=tstep_out, dump_cadence=0)
sim.save()

sim.set_parameter(simdir="swifter_sim",codename="Swifter",init_cond_file_type="ASCII",output_file_type="REAL8")
sim.save()
229 changes: 229 additions & 0 deletions examples/rmvs_swifter_comparison/8pl_16tp_encounters/swiftest_vs_swifter.ipynb
View file

Large diffs are not rendered by default.

9 changes: 9 additions & 0 deletions examples/rmvs_swifter_comparison/8pl_16tp_encounters/swiftest_vs_swifter.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
import swiftest
import numpy as np

swiftersim = swiftest.Simulation(simdir="swifter_sim", read_data=True, codename="Swifter")
swiftestsim = swiftest.Simulation(simdir="swiftest_sim",read_data=True)
swiftestsim.data = swiftestsim.data.swap_dims({"name" : "id"})
swiftdiff = swiftestsim.data - swiftersim.data
swiftdiff['rh'].plot(x="time",hue="id",col="space")
swiftdiff['vh'].plot(x="time",hue="id",col="space")
52 changes: 26 additions & 26 deletions src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,17 +13,37 @@

# Add the source files
SET(STRICT_MATH_FILES
${SRC}/swiftest/swiftest_kick.f90
${SRC}/helio/helio_kick.f90
${SRC}/rmvs/rmvs_kick.f90
${SRC}/symba/symba_kick.f90
${SRC}/whm/whm_kick.f90
${SRC}/swiftest/swiftest_user.f90
${SRC}/collision/collision_generate.f90
${SRC}/fraggle/fraggle_generate.f90
${SRC}/fraggle/fraggle_util.f90
${SRC}/fraggle/fraggle_module.f90
${SRC}/helio/helio_drift.f90
${SRC}/helio/helio_gr.f90
${SRC}/helio/helio_kick.f90
${SRC}/helio/helio_step.f90
${SRC}/misc/lambda_function_module.f90
${SRC}/misc/solver_module.f90
${SRC}/operator/operator_module.f90
${SRC}/operator/operator_cross.f90
${SRC}/operator/operator_mag.f90
${SRC}/operator/operator_unit.f90
${SRC}/rmvs/rmvs_kick.f90
${SRC}/rmvs/rmvs_step.f90
${SRC}/swiftest/swiftest_drift.f90
${SRC}/swiftest/swiftest_gr.f90
${SRC}/swiftest/swiftest_kick.f90
${SRC}/swiftest/swiftest_user.f90
${SRC}/swiftest/swiftest_obl.f90
${SRC}/swiftest/swiftest_orbel.f90
${SRC}/symba/symba_drift.f90
${SRC}/symba/symba_gr.f90
${SRC}/symba/symba_kick.f90
${SRC}/symba/symba_step.f90
${SRC}/whm/whm_coord.f90
${SRC}/whm/whm_drift.f90
${SRC}/whm/whm_gr.f90
${SRC}/whm/whm_kick.f90
${SRC}/whm/whm_step.f90
)

SET(FAST_MATH_FILES
Expand All @@ -33,52 +53,32 @@ SET(FAST_MATH_FILES
${SRC}/misc/io_progress_bar_module.f90
${SRC}/encounter/encounter_module.f90
${SRC}/collision/collision_module.f90
${SRC}/operator/operator_module.f90
${SRC}/walltime/walltime_module.f90
${SRC}/swiftest/swiftest_module.f90
${SRC}/whm/whm_module.f90
${SRC}/rmvs/rmvs_module.f90
${SRC}/helio/helio_module.f90
${SRC}/symba/symba_module.f90
${SRC}/collision/collision_check.f90
${SRC}/collision/collision_generate.f90
${SRC}/collision/collision_io.f90
${SRC}/collision/collision_regime.f90
${SRC}/collision/collision_resolve.f90
${SRC}/collision/collision_util.f90
${SRC}/encounter/encounter_check.f90
${SRC}/encounter/encounter_io.f90
${SRC}/encounter/encounter_util.f90
${SRC}/helio/helio_drift.f90
${SRC}/helio/helio_gr.f90
${SRC}/helio/helio_step.f90
${SRC}/helio/helio_util.f90
${SRC}/netcdf_io/netcdf_io_implementations.f90
${SRC}/operator/operator_cross.f90
${SRC}/operator/operator_mag.f90
${SRC}/operator/operator_unit.f90
${SRC}/rmvs/rmvs_discard.f90
${SRC}/rmvs/rmvs_encounter_check.f90
${SRC}/rmvs/rmvs_step.f90
${SRC}/rmvs/rmvs_util.f90
${SRC}/swiftest/swiftest_discard.f90
${SRC}/swiftest/swiftest_drift.f90
${SRC}/swiftest/swiftest_gr.f90
${SRC}/swiftest/swiftest_io.f90
${SRC}/swiftest/swiftest_obl.f90
${SRC}/swiftest/swiftest_orbel.f90
${SRC}/swiftest/swiftest_util.f90
${SRC}/symba/symba_discard.f90
${SRC}/symba/symba_drift.f90
${SRC}/symba/symba_encounter_check.f90
${SRC}/symba/symba_gr.f90
${SRC}/symba/symba_step.f90
${SRC}/symba/symba_util.f90
${SRC}/walltime/walltime_implementations.f90
${SRC}/whm/whm_coord.f90
${SRC}/whm/whm_drift.f90
${SRC}/whm/whm_gr.f90
${SRC}/whm/whm_step.f90
${SRC}/whm/whm_util.f90
${SRC}/swiftest/swiftest_driver.f90
)
Expand Down
4 changes: 4 additions & 0 deletions src/rmvs/rmvs_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,6 +202,10 @@ subroutine rmvs_step_out(cb, pl, tp, nbody_system, param, t, dt)
call tp%step(nbody_system, param, outer_time, dto)
tp%lfirst = lfirsttp
else
if (param%loblatecb) then
call swiftest_obl_acc(npl, cb%Gmass, cb%j2rp2, cb%j4rp4, pl%rbeg, pl%lmask, pl%outer(outer_index-1)%aobl, pl%Gmass, cb%aoblbeg)
call swiftest_obl_acc(npl, cb%Gmass, cb%j2rp2, cb%j4rp4, pl%rend, pl%lmask, pl%outer(outer_index)%aobl, pl%Gmass, cb%aoblend)
end if
call tp%step(nbody_system, param, outer_time, dto)
end if
do j = 1, npl
Expand Down
4 changes: 3 additions & 1 deletion src/rmvs/rmvs_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -298,13 +298,15 @@ module subroutine rmvs_util_setup_pl(self, n, param)
allocate(pl%outer(i)%v(NDIM, n))
pl%outer(i)%x(:,:) = 0.0_DP
pl%outer(i)%v(:,:) = 0.0_DP
allocate(pl%outer(i)%aobl(NDIM, n))
pl%outer(i)%aobl(:,:) = 0.0_DP
end do
do i = 0, NTPHENC
allocate(pl%inner(i)%x(NDIM, n))
allocate(pl%inner(i)%v(NDIM, n))
allocate(pl%inner(i)%aobl(NDIM, n))
pl%inner(i)%x(:,:) = 0.0_DP
pl%inner(i)%v(:,:) = 0.0_DP
allocate(pl%inner(i)%aobl(NDIM, n))
pl%inner(i)%aobl(:,:) = 0.0_DP
end do
! if (param%ltides) then
Expand Down
2 changes: 0 additions & 2 deletions src/swiftest/swiftest_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -983,7 +983,6 @@ module subroutine swiftest_io_netcdf_open(self, param, readonly)
! Optional variables The User Doesn't Need to Know About
status = nf90_inq_varid(nc%id, nc%mass_varname, nc%mass_varid)
status = nf90_inq_varid(nc%id, nc%rhill_varname, nc%rhill_varid)
param%lrhill_present = (status == NF90_NOERR)
status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid)
status = nf90_inq_varid(nc%id, nc%status_varname, nc%status_varid)
status = nf90_inq_varid(nc%id, nc%ntp_varname, nc%ntp_varid)
Expand Down Expand Up @@ -1232,7 +1231,6 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier
status = nf90_get_var(nc%id, nc%rhill_varid, rtemp, start=[1, tslot], count=[idmax,1])
if (status == NF90_NOERR) then
pl%rhill(:) = pack(rtemp, plmask)
param%lrhill_present = .true.
end if
end if

Expand Down
2 changes: 1 addition & 1 deletion src/swiftest/swiftest_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -430,7 +430,7 @@ pure module subroutine swiftest_kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl,
! Internals
real(DP) :: fac

fac = GMpl * sqrt(1.0_DP / (rji2*rji2*rji2))
fac = GMpl / (rji2*sqrt(rji2))
ax = ax - fac * xr
ay = ay - fac * yr
az = az - fac * zr
Expand Down
18 changes: 12 additions & 6 deletions src/swiftest/swiftest_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,7 +146,6 @@ module swiftest
procedure :: read_in => swiftest_io_read_in_body !! Read in body initial conditions from an ascii file
procedure :: write_frame => swiftest_io_netcdf_write_frame_body !! I/O routine for writing out a single frame of time-series data for all bodies in the nbody_system in NetCDF format
procedure :: write_info => swiftest_io_netcdf_write_info_body !! Dump contents of particle information metadata to file
procedure :: accel_obl => swiftest_obl_acc_body !! Compute the barycentric accelerations of bodies due to the oblateness of the central body
procedure :: el2xv => swiftest_orbel_el2xv_vec !! Convert orbital elements to position and velocity vectors
procedure :: xv2el => swiftest_orbel_xv2el_vec !! Convert position and velocity vectors to orbital elements
procedure :: setup => swiftest_util_setup_body !! A constructor that sets the number of bodies and allocates all allocatable arrays
Expand Down Expand Up @@ -995,11 +994,18 @@ pure module subroutine swiftest_kick_getacch_int_one_tp(rji2, xr, yr, zr, Gmpl,
real(DP), intent(inout) :: ax, ay, az !! Acceleration vector components of test particle
end subroutine swiftest_kick_getacch_int_one_tp

module subroutine swiftest_obl_acc_body(self, nbody_system)
implicit none
class(swiftest_body), intent(inout) :: self !! Swiftest body object
class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
end subroutine swiftest_obl_acc_body
module subroutine swiftest_obl_acc(n, GMcb, j2rp2, j4rp4, rh, lmask, aobl, GMpl, aoblcb)
implicit none
integer(I4B), intent(in) :: n !! Number of bodies
real(DP), intent(in) :: GMcb !! Central body G*Mass
real(DP), intent(in) :: j2rp2 !! J2 * R**2 for the central body
real(DP), intent(in) :: j4rp4 !! J4 * R**4 for the central body
real(DP), dimension(:,:), intent(in) :: rh !! Heliocentric positions of bodies
logical, dimension(:), intent(in) :: lmask !! Logical mask of bodies to compute aobl
real(DP), dimension(:,:), intent(out) :: aobl !! Barycentric acceleration of bodies due to central body oblateness
real(DP), dimension(:), intent(in), optional :: GMpl !! Masses of input bodies if they are not test particles
real(DP), dimension(:), intent(out), optional :: aoblcb !! Barycentric acceleration of central body (only needed if input bodies are massive)
end subroutine swiftest_obl_acc

module subroutine swiftest_obl_acc_pl(self, nbody_system)
implicit none
Expand Down
63 changes: 36 additions & 27 deletions src/swiftest/swiftest_obl.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@

submodule (swiftest) s_swiftest_obl
contains
module subroutine swiftest_obl_acc_body(self, nbody_system)
module subroutine swiftest_obl_acc(n, GMcb, j2rp2, j4rp4, rh, lmask, aobl, GMpl, aoblcb)
!! author: David A. Minton
!!
!! Compute the barycentric accelerations of bodies due to the oblateness of the central body
Expand All @@ -19,33 +19,46 @@ module subroutine swiftest_obl_acc_body(self, nbody_system)
!! Adapted from Hal Levison's Swift routine obl_acc.f and obl_acc_tp.f
implicit none
! Arguments
class(swiftest_body), intent(inout) :: self !! Swiftest body object
class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
integer(I4B), intent(in) :: n !! Number of bodies
real(DP), intent(in) :: GMcb !! Central body G*Mass
real(DP), intent(in) :: j2rp2 !! J2 * R**2 for the central body
real(DP), intent(in) :: j4rp4 !! J4 * R**4 for the central body
real(DP), dimension(:,:), intent(in) :: rh !! Heliocentric positions of bodies
logical, dimension(:), intent(in) :: lmask !! Logical mask of bodies to compute aobl
real(DP), dimension(:,:), intent(out) :: aobl !! Barycentric acceleration of bodies due to central body oblateness
real(DP), dimension(:), intent(in), optional :: GMpl !! Masses of input bodies if they are not test particles
real(DP), dimension(:), intent(out), optional :: aoblcb !! Barycentric acceleration of central body (only needed if input bodies are massive)
! Internals
integer(I4B) :: i
real(DP) :: r2, irh, rinv2, t0, t1, t2, t3, fac1, fac2

if (self%nbody == 0) return

associate(n => self%nbody, cb => nbody_system%cb)
self%aobl(:,:) = 0.0_DP
do concurrent(i = 1:n, self%lmask(i))
r2 = dot_product(self%rh(:, i), self%rh(:, i))
irh = 1.0_DP / sqrt(r2)
rinv2 = irh**2
t0 = -cb%Gmass * rinv2 * rinv2 * irh
t1 = 1.5_DP * cb%j2rp2
t2 = self%rh(3, i) * self%rh(3, i) * rinv2
t3 = 1.875_DP * cb%j4rp4 * rinv2
fac1 = t0 * (t1 - t3 - (5 * t1 - (14.0_DP - 21.0_DP * t2) * t3) * t2)
fac2 = 2 * t0 * (t1 - (2.0_DP - (14.0_DP * t2 / 3.0_DP)) * t3)
self%aobl(:, i) = fac1 * self%rh(:, i)
self%aobl(3, i) = fac2 * self%rh(3, i) + self%aobl(3, i)
if (n == 0) return

aobl(:,:) = 0.0_DP
do concurrent(i = 1:n, lmask(i))
r2 = dot_product(rh(:, i), rh(:, i))
irh = 1.0_DP / sqrt(r2)
rinv2 = irh**2
t0 = -GMcb * rinv2 * rinv2 * irh
t1 = 1.5_DP * j2rp2
t2 = rh(3, i) * rh(3, i) * rinv2
t3 = 1.875_DP * j4rp4 * rinv2
fac1 = t0 * (t1 - t3 - (5 * t1 - (14.0_DP - 21.0_DP * t2) * t3) * t2)
fac2 = 2 * t0 * (t1 - (2.0_DP - (14.0_DP * t2 / 3.0_DP)) * t3)
aobl(:, i) = fac1 * rh(:, i)
aobl(3, i) = fac2 * rh(3, i) + aobl(3, i)
end do

if (present(GMpl) .and. present(aoblcb)) then
aoblcb(:) = 0.0_DP
do i = n, 1, -1
if (lmask(i)) aoblcb(:) = aoblcb(:) - GMpl(i) * aobl(:, i) / GMcb
end do
end associate
end if

return

end subroutine swiftest_obl_acc_body
end subroutine swiftest_obl_acc


module subroutine swiftest_obl_acc_pl(self, nbody_system)
Expand All @@ -65,11 +78,7 @@ module subroutine swiftest_obl_acc_pl(self, nbody_system)
if (self%nbody == 0) return

associate(pl => self, npl => self%nbody, cb => nbody_system%cb)
call swiftest_obl_acc_body(pl, nbody_system)
cb%aobl(:) = 0.0_DP
do i = npl, 1, -1
if (pl%lmask(i)) cb%aobl(:) = cb%aobl(:) - pl%Gmass(i) * pl%aobl(:, i) / cb%Gmass
end do
call swiftest_obl_acc(npl, cb%Gmass, cb%j2rp2, cb%j4rp4, pl%rh, pl%lmask, pl%aobl, pl%Gmass, cb%aobl)

do concurrent(i = 1:npl, pl%lmask(i))
pl%ah(:, i) = pl%ah(:, i) + pl%aobl(:, i) - cb%aobl(:)
Expand Down Expand Up @@ -99,7 +108,7 @@ module subroutine swiftest_obl_acc_tp(self, nbody_system)
if (self%nbody == 0) return

associate(tp => self, ntp => self%nbody, cb => nbody_system%cb)
call swiftest_obl_acc_body(tp, nbody_system)
call swiftest_obl_acc(ntp, cb%Gmass, cb%j2rp2, cb%j4rp4, tp%rh, tp%lmask, tp%aobl)
if (nbody_system%lbeg) then
aoblcb = cb%aoblbeg
else
Expand Down
Loading

0 comments on commit d0c8acb

Please sign in to comment.