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

Commit

Permalink
Merge branch 'debug' into DocumentationTesting
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Sep 20, 2021
2 parents 76d235c + 41265b7 commit 7335412
Show file tree
Hide file tree
Showing 24 changed files with 493 additions and 5,524 deletions.
11 changes: 6 additions & 5 deletions Makefile.Defines
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,11 +53,11 @@ VTUNE_FLAGS = -g -O2 -qopt-report=5 -simd -shared-intel -qopenmp -debug inline-d

IDEBUG = -O0 -init=snan,arrays -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin
STRICTREAL = -fp-model strict -prec-div -prec-sqrt -assume protect-parens
SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -fp-model no-except
SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -fp-model no-except -fma
PAR = -qopenmp -parallel
HEAPARR = -heap-arrays 4194304
OPTREPORT = -qopt-report=5
IPRODUCTION = -no-wrap-margin -O3 $(PAR) $(SIMDVEC) #$(HEAPARR)
IPRODUCTION = -no-wrap-margin -O3 -ipo -qopt-prefetch=0 -sox $(PAR) $(SIMDVEC) #$(HEAPARR)

#gfortran flags
GDEBUG = -g -Og -fbacktrace -fbounds-check -ffree-line-length-none
Expand All @@ -68,10 +68,11 @@ GPRODUCTION = -O2 -ffree-line-length-none $(GPAR)


#FFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR)
FFLAGS = $(IPRODUCTION) $(STRICTREAL) $(OPTREPORT)
FFASTFLAGS = $(IPRODUCTION) -fp-model fast $(OPTREPORT)
#FFASTFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR)
FFLAGS = $(IPRODUCTION) $(STRICTREAL)
FFASTFLAGS = $(IPRODUCTION) -fp-model fast
FORTRAN = ifort
#AR = xiar
AR = xiar

#FORTRAN = gfortran
#FFLAGS = $(GDEBUG) $(GMEM) $(GPAR)
Expand Down
6 changes: 0 additions & 6 deletions examples/symba_mars_disk/Untitled.ipynb
View file

This file was deleted.

13 changes: 7 additions & 6 deletions examples/symba_mars_disk/param.in
View file
Original file line number Diff line number Diff line change
@@ -1,17 +1,17 @@
!Parameter file for the SyMBA-RINGMOONS test
T0 0.0
TSTOP 60000.0
TSTOP 12000.0
DT 600.0
CB_IN cb.in
PL_IN mars.in
TP_IN tp.in
IN_TYPE ASCII
ISTEP_OUT 100
ISTEP_DUMP 100
BIN_OUT bin.nc
ISTEP_OUT 1
ISTEP_DUMP 1
BIN_OUT bin.nc
PARTICLE_OUT particle.dat
OUT_TYPE REAL8
OUT_FORM XV
OUT_TYPE NETCDF_DOUBLE
OUT_FORM XVEL
OUT_STAT REPLACE
CHK_CLOSE yes
CHK_RMIN 3389500.0
Expand All @@ -31,3 +31,4 @@ MU2KG 1.0
DU2M 1.0
TU2S 1.0
SEED 2 3080983 2220830
FLATTEN_INTERACTIONS yes
33 changes: 33 additions & 0 deletions examples/symba_mars_disk/param.real8.in
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
!Parameter file for the SyMBA-RINGMOONS test
T0 0.0
TSTOP 6000.0
DT 600.0
CB_IN cb.in
PL_IN mars.in
TP_IN tp.in
IN_TYPE ASCII
ISTEP_OUT 1
ISTEP_DUMP 1
BIN_OUT bin.dat
PARTICLE_OUT particle.dat
OUT_TYPE REAL8
OUT_FORM EL
OUT_STAT REPLACE
CHK_CLOSE yes
CHK_RMIN 3389500.0
CHK_RMAX 3389500000.0
CHK_EJECT 3389500000.0
CHK_QMIN 3389500.0
CHK_QMIN_COORD HELIO
CHK_QMIN_RANGE 3389500.0 338950000000.0
EXTRA_FORCE no
RHILL_PRESENT yes
GMTINY 10000.0
MIN_GMFRAG 1000.0
ENERGY yes
FRAGMENTATION yes
ROTATION yes
MU2KG 1.0
DU2M 1.0
TU2S 1.0
SEED 2 3080983 2220830
5,414 changes: 73 additions & 5,341 deletions examples/symba_mars_disk/testnetcdf.ipynb
View file

Large diffs are not rendered by default.

6 changes: 3 additions & 3 deletions src/helio/helio_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ module subroutine helio_kick_getacch_pl(self, system, param, t, lbeg)
if (self%nbody == 0) return

associate(cb => system%cb, pl => self, npl => self%nbody)
call pl%accel_int()
call pl%accel_int(param)
if (param%loblatecb) then
call pl%accel_obl(system)
if (lbeg) then
Expand Down Expand Up @@ -65,9 +65,9 @@ module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg)
associate(tp => self, cb => system%cb, pl => system%pl, npl => system%pl%nbody)
system%lbeg = lbeg
if (system%lbeg) then
call tp%accel_int(pl%Gmass(1:npl), pl%xbeg(:,1:npl), npl)
call tp%accel_int(param, pl%Gmass(1:npl), pl%xbeg(:,1:npl), npl)
else
call tp%accel_int(pl%Gmass(1:npl), pl%xend(:,1:npl), npl)
call tp%accel_int(param, pl%Gmass(1:npl), pl%xend(:,1:npl), npl)
end if
if (param%loblatecb) call tp%accel_obl(system)
if (param%lextra_force) call tp%accel_user(system, param, t, lbeg)
Expand Down
9 changes: 7 additions & 2 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -601,6 +601,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
case ("TIDES")
call io_toupper(param_value)
if (param_value == "YES" .or. param_value == 'T') param%ltides = .true.
case ("FLATTEN_INTERACTIONS")
call io_toupper(param_value)
if (param_value == "NO" .or. param_value == 'F') param%lflatten_interactions = .false.
case ("FIRSTKICK")
call io_toupper(param_value)
if (param_value == "NO" .or. param_value == 'F') param%lfirstkick = .false.
Expand Down Expand Up @@ -755,9 +758,10 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
if ((param%qmin_alo > 0.0_DP) .or. (param%qmin_ahi > 0.0_DP)) write(*,*) "CHK_QMIN_RANGE = ",param%qmin_alo, param%qmin_ahi
write(*,*) "EXTRA_FORCE = ",param%lextra_force
write(*,*) "RHILL_PRESENT = ",param%lrhill_present
write(*,*) "ROTATION = ", param%lrotation
write(*,*) "TIDES = ", param%ltides
write(*,*) "ROTATION = ", param%lrotation
write(*,*) "TIDES = ", param%ltides
write(*,*) "ENERGY = ",param%lenergy
write(*,*) "FLATTEN_INTERACTIONS = ",param%lflatten_interactions
if (param%lenergy) write(*,*) "ENERGY_OUT = ",trim(adjustl(param%energy_out))
write(*,*) "MU2KG = ",param%MU2KG
write(*,*) "TU2S = ",param%TU2S
Expand Down Expand Up @@ -899,6 +903,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
write(param_name, *) "GR"; write(param_value, Lfmt) param%lgr; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value))
write(param_name, *) "ROTATION"; write(param_value, Lfmt) param%lrotation; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value))
write(param_name, *) "TIDES"; write(param_value, Lfmt) param%ltides; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value))
write(param_name, *) "FLATTEN_INTERACTIONS"; write(param_value, Lfmt) param%lflatten_interactions; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value))

if (param%lenergy) then
write(param_name, *) "FIRSTENERGY"; write(param_value, Lfmt) param%lfirstenergy; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value))
Expand Down
80 changes: 68 additions & 12 deletions src/kick/kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
use swiftest
contains

module subroutine kick_getacch_int_pl(self)
module subroutine kick_getacch_int_pl(self, param)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations of massive bodies
Expand All @@ -11,15 +11,20 @@ module subroutine kick_getacch_int_pl(self)
!! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f90
implicit none
! Arguments
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters

call kick_getacch_int_all_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)
if (param%lflatten_interactions) then
call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)
else
call kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%xh, self%Gmass, self%radius, self%ah)
end if

return
end subroutine kick_getacch_int_pl


module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
module subroutine kick_getacch_int_tp(self, param, GMpl, xhp, npl)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations of test particles by massive bodies
Expand All @@ -28,10 +33,11 @@ module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
!! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int_tp.f90
implicit none
! Arguments
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses
real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors
integer(I4B), intent(in) :: npl !! Number of active massive bodies
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters
real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses
real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors
integer(I4B), intent(in) :: npl !! Number of active massive bodies

if ((self%nbody == 0) .or. (npl == 0)) return

Expand All @@ -41,10 +47,11 @@ module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
end subroutine kick_getacch_int_tp


module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc)
module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization
!! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization.
!! This is the flattened (single loop) version that uses the k_plpl interaction pair index array
!!
!! Adapted from Hal Levison's Swift routine getacch_ah3.f
!! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f9
Expand Down Expand Up @@ -88,7 +95,56 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius,
end do

return
end subroutine kick_getacch_int_all_pl
end subroutine kick_getacch_int_all_flat_pl


module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization.
!! This is the upper triangular matrix (double loop) version.
!!
!! Adapted from Hal Levison's Swift routine getacch_ah3.f
!! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f9
implicit none
integer(I4B), intent(in) :: npl !! Total number of massive bodies
integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies
real(DP), dimension(:,:), intent(in) :: x !! Position vector array
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
! Internals
real(DP), dimension(NDIM,npl) :: ahi, ahj
integer(I4B) :: i, j
real(DP) :: rji2, rlim2
real(DP) :: xr, yr, zr

ahi(:,:) = 0.0_DP
ahj(:,:) = 0.0_DP

!$omp parallel do default(private) schedule(static)&
!$omp shared(npl, nplm, x, Gmass, radius) &
!$omp lastprivate(rji2, rlim2, xr, yr, zr) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do i = 1, nplm
do concurrent(j = i+1:npl)
xr = x(1, j) - x(1, i)
yr = x(2, j) - x(2, i)
zr = x(3, j) - x(3, i)
rji2 = xr**2 + yr**2 + zr**2
rlim2 = (radius(i) + radius(j))**2
if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
end do
end do
!$omp end parallel do

do concurrent(i = 1:npl)
acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i)
end do

return
end subroutine kick_getacch_int_all_triangular_pl


module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc)
Expand All @@ -112,7 +168,7 @@ module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc)
integer(I4B) :: i, j

!$omp parallel do default(private) schedule(static)&
!$omp shared(npl, ntp, lmask, xtp, xpl) &
!$omp shared(npl, ntp, lmask, xtp, xpl, GMpl) &
!$omp reduction(-:acc)
do i = 1, ntp
if (lmask(i)) then
Expand Down
11 changes: 6 additions & 5 deletions src/modules/rmvs_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,12 +122,13 @@ module subroutine rmvs_discard_tp(self, system, param)
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
end subroutine rmvs_discard_tp

module function rmvs_encounter_check_tp(self, system, dt) result(lencounter)
module function rmvs_encounter_check_tp(self, param, system, dt) result(lencounter)
implicit none
class(rmvs_tp), intent(inout) :: self !! RMVS test particle object
class(rmvs_nbody_system), intent(inout) :: system !! RMVS nbody system object
real(DP), intent(in) :: dt !! step size
logical :: lencounter !! Returns true if there is at least one close encounter
class(rmvs_tp), intent(inout) :: self !! RMVS test particle object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
class(rmvs_nbody_system), intent(inout) :: system !! RMVS nbody system object
real(DP), intent(in) :: dt !! step size
logical :: lencounter !! Returns true if there is at least one close encounter
end function rmvs_encounter_check_tp

module subroutine rmvs_io_write_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2, enc_out)
Expand Down
32 changes: 22 additions & 10 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,7 @@ module swiftest_classes
logical :: loblatecb = .false. !! Calculate acceleration from oblate central body (automatically turns true if nonzero J2 is input)
logical :: lrotation = .false. !! Include rotation states of big bodies
logical :: ltides = .false. !! Include tidal dissipation
logical :: lflatten_interactions = .true. !! Use the flattened upper triangular matrix for pl-pl interactions (turning this on improves the speed but uses more memory)

! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values)
real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy
Expand Down Expand Up @@ -432,7 +433,6 @@ module swiftest_classes
integer(I4B) :: nenc !! Total number of encounters
logical, dimension(:), allocatable :: lvdotr !! relative vdotr flag
integer(I4B), dimension(:), allocatable :: status !! status of the interaction
integer(I8B), dimension(:), allocatable :: kidx !! index value of the encounter from the master k_plpl encounter list
integer(I4B), dimension(:), allocatable :: index1 !! position of the first body in the encounter
integer(I4B), dimension(:), allocatable :: index2 !! position of the second body in the encounter
integer(I4B), dimension(:), allocatable :: id1 !! id of the first body in the encounter
Expand Down Expand Up @@ -840,20 +840,22 @@ module subroutine io_write_hdr_system(self, iu, param)
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
end subroutine io_write_hdr_system

module subroutine kick_getacch_int_pl(self)
module subroutine kick_getacch_int_pl(self, param)
implicit none
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters
end subroutine kick_getacch_int_pl

module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
module subroutine kick_getacch_int_tp(self, param, GMpl, xhp, npl)
implicit none
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle
real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses
real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors
integer(I4B), intent(in) :: npl !! Number of active massive bodies
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters
real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses
real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors
integer(I4B), intent(in) :: npl !! Number of active massive bodies
end subroutine kick_getacch_int_tp

module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc)
module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc)
implicit none
integer(I4B), intent(in) :: npl !! Number of massive bodies
integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute
Expand All @@ -862,7 +864,17 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius,
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
end subroutine kick_getacch_int_all_pl
end subroutine kick_getacch_int_all_flat_pl

module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc)
implicit none
integer(I4B), intent(in) :: npl !! Total number of massive bodies
integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies
real(DP), dimension(:,:), intent(in) :: x !! Position vector array
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
end subroutine kick_getacch_int_all_triangular_pl

module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc)
implicit none
Expand Down
Loading

0 comments on commit 7335412

Please sign in to comment.