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

Commit

Permalink
Put together basic SyMBA infrastructure
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Jul 9, 2021
1 parent 8b32ff2 commit e186156
Show file tree
Hide file tree
Showing 17 changed files with 904 additions and 557 deletions.
8 changes: 4 additions & 4 deletions Makefile.Defines
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,13 +65,13 @@ 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 = $(IDEBUG) $(HEAPARR)
#FFLAGS = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(SIMDVEC) $(PAR)
FORTRAN = ifort
#FORTRAN = ifort
#AR = xiar

#FORTRAN = gfortran
#FFLAGS = -ffree-line-length-none $(GDEBUG) #$(GMEM)
FORTRAN = gfortran
FFLAGS = -ffree-line-length-none $(GDEBUG) #$(GMEM)
AR = ar

# DO NOT include in CFLAGS the "-c" option to compile object only
Expand Down
10 changes: 5 additions & 5 deletions src/discard/discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,19 +88,19 @@ subroutine discard_sun_tp(tp, system, param)
rh2 = dot_product(tp%xh(:, i), tp%xh(:, i))
if ((param%rmax >= 0.0_DP) .and. (rh2 > rmax2)) then
tp%status(i) = DISCARDED_RMAX
write(*, *) "Particle ", tp%name(i), " too far from sun at t = ", t
write(*, *) "Particle ", tp%id(i), " too far from sun at t = ", t
tp%ldiscard(i) = .true.
else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then
tp%status(i) = DISCARDED_RMIN
write(*, *) "Particle ", tp%name(i), " too close to sun at t = ", t
write(*, *) "Particle ", tp%id(i), " too close to sun at t = ", t
tp%ldiscard(i) = .true.
else if (param%rmaxu >= 0.0_DP) then
rb2 = dot_product(tp%xb(:, i), tp%xb(:, i))
vb2 = dot_product(tp%vb(:, i), tp%vb(:, i))
energy = 0.5_DP * vb2 - msys / sqrt(rb2)
if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then
tp%status(i) = DISCARDED_RMAXU
write(*, *) "Particle ", tp%name(i), " is unbound and too far from barycenter at t = ", t
write(*, *) "Particle ", tp%id(i), " is unbound and too far from barycenter at t = ", t
tp%ldiscard(i) = .true.
end if
end if
Expand Down Expand Up @@ -150,7 +150,7 @@ subroutine discard_peri_tp(tp, system, param)
(tp%atp(i) <= param%qmin_ahi) .and. &
(tp%peri(i) <= param%qmin)) then
tp%status(i) = DISCARDED_PERI
write(*, *) "Particle ", tp%name(i), " perihelion distance too small at t = ", t
write(*, *) "Particle ", tp%id(i), " perihelion distance too small at t = ", t
tp%ldiscard(i) = .true.
end if
end if
Expand Down Expand Up @@ -191,7 +191,7 @@ subroutine discard_pl_tp(tp, system, param)
if (isp /= 0) then
tp%status(i) = DISCARDED_PLR
pl%ldiscard(j) = .true.
write(*, *) "Particle ", tp%name(i), " too close to massive body ", pl%name(j), " at t = ", t
write(*, *) "Particle ", tp%id(i), " too close to massive body ", pl%id(j), " at t = ", t
tp%ldiscard(i) = .true.
exit
end if
Expand Down
4 changes: 2 additions & 2 deletions src/helio/helio_drift.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ module subroutine helio_drift_pl(self, system, param, dt)
dtp(1:npl), iflag(1:npl))
if (any(iflag(1:npl) /= 0)) then
do i = 1, npl
write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!"
write(*, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!"
write(*, *) pl%xh(:,i)
write(*, *) pl%vb(:,i)
write(*, *) " stopping "
Expand Down Expand Up @@ -136,7 +136,7 @@ module subroutine helio_drift_tp(self, system, param, dt)
if (any(iflag(1:ntp) /= 0)) then
tp%status = DISCARDED_DRIFTERR
do i = 1, ntp
if (iflag(i) /= 0) write(*, *) "Particle ", tp%name(i), " lost due to error in Danby drift"
if (iflag(i) /= 0) write(*, *) "Particle ", tp%id(i), " lost due to error in Danby drift"
end do
end if
end associate
Expand Down
287 changes: 81 additions & 206 deletions src/io/io.f90
View file

Large diffs are not rendered by default.

4 changes: 2 additions & 2 deletions src/modules/helio_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ module helio_classes
!! Definition of classes and methods specific to the Democratic Heliocentric Method
!! Adapted from David E. Kaufmann's Swifter routine: helio.f90
use swiftest_globals
use swiftest_classes, only : swiftest_cb, swiftest_pl, swiftest_tp
use swiftest_classes, only : swiftest_cb, swiftest_pl, swiftest_tp, swiftest_nbody_system
use whm_classes, only : whm_nbody_system
implicit none

Expand Down Expand Up @@ -171,7 +171,7 @@ module subroutine helio_step_pl(self, system, param, t, dt)
end subroutine helio_step_pl

module subroutine helio_step_tp(self, system, param, t, dt)
use swiftest_classes, only : swiftest_cb, swiftest_parameters
use swiftest_classes, only : swiftest_cb, swiftest_parameters, swiftest_nbody_system
implicit none
class(helio_tp), intent(inout) :: self !! Helio test particle object
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
Expand Down
2 changes: 1 addition & 1 deletion src/modules/rmvs_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ module rmvs_classes
!*******************************************************************************************************************************

!> RMVS massive body particle class
type, private, extends(whm_pl) :: rmvs_pl
type, public, extends(whm_pl) :: rmvs_pl
integer(I4B), dimension(:), allocatable :: nenc !! number of test particles encountering planet this full rmvs time step
integer(I4B), dimension(:), allocatable :: tpenc1P !! index of first test particle encountering planet
integer(I4B), dimension(:), allocatable :: plind ! Connects the planetocentric indices back to the heliocentric planet list
Expand Down
Loading

0 comments on commit e186156

Please sign in to comment.