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

Commit

Permalink
Merge branch 'OOPSymba' into OOPTides
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Jul 24, 2021
2 parents f808a7c + 9a201c4 commit 7878958
Show file tree
Hide file tree
Showing 3 changed files with 74 additions and 15 deletions.
9 changes: 5 additions & 4 deletions src/modules/symba_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -175,21 +175,22 @@ module subroutine symba_discard_tp(self, system, param)
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
end subroutine symba_discard_tp

module function symba_encounter_check_pl(self, system, dt, irec) result(lencounter)
module function symba_encounter_check_pl(self, system, dt, irec) result(lany_encounter)
implicit none
class(symba_pl), intent(inout) :: self !! SyMBA test particle object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
logical :: lencounter !! Returns true if there is at least one close encounter
integer(I4B), intent(in) :: irec !! Current recursion level
logical :: lany_encounter !! Returns true if there is at least one close encounter
end function symba_encounter_check_pl

module function symba_encounter_check_tp(self, system, dt) result(lencounter)
module function symba_encounter_check_tp(self, system, dt, irec) result(lany_encounter)
implicit none
class(symba_tp), intent(inout) :: self !! SyMBA test particle object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
logical :: lencounter !! Returns true if there is at least one close encounter
integer(I4B), intent(in) :: irec !! Current recursion level
logical :: lany_encounter !! Returns true if there is at least one close encounter
end function symba_encounter_check_tp

module subroutine symba_io_dump_particle_info(self, param, msg)
Expand Down
69 changes: 64 additions & 5 deletions src/symba/symba_encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,10 +12,9 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc
logical :: lany_encounter !! Returns true if there is at least one close encounter
! Internals
real(DP) :: r2crit, vdotr, r2, v2, tmin, r2min, term2
integer(I4B) :: j, nenc_old
integer(I4B) :: nenc_old
integer(I8B) :: k
real(DP), dimension(NDIM) :: xr, vr
integer(I4B), dimension(:,:), allocatable :: ind
logical, dimension(:), allocatable :: lencounter, loc_lvdotr

associate(pl => self, npl => self%nbody, nplpl => self%nplpl)
Expand Down Expand Up @@ -71,16 +70,76 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc
return
end function symba_encounter_check_pl

module function symba_encounter_check_tp(self, system, dt) result(lencounter)
module function symba_encounter_check_tp(self, system, dt, irec) result(lany_encounter)
implicit none
! Arguments
class(symba_tp), intent(inout) :: self !! SyMBA test particle object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
! Result
logical :: lencounter !! Returns true if there is at least one close encounter
logical :: lany_encounter !! Returns true if there is at least one close encounter
! Internals
real(DP) :: r2crit, vdotr, r2, v2, tmin, r2min, term2
integer(I4B) :: i, j, nenc_old
real(DP), dimension(NDIM) :: xr, vr
logical, dimension(:,:), allocatable :: lencounter, loc_lvdotr

associate(tp => self, ntp => self%nbody, pl => system%pl, npl => system%pl%nbody)
allocate(lencounter(npl, ntp), loc_lvdotr(npl, ntp))
lencounter(:,:) = .false.

term2 = RHSCALE * (RSHELL**irec)

do j = 1, ntp
do i = 1, npl
xr(:) = tp%xh(:, j) - pl%xh(:, i)
r2 = dot_product(xr(:), xr(:))
r2crit = (pl%rhill(i) * term2)**2
vr(:) = tp%vh(:, j) - pl%vh(:, i)
vdotr = dot_product(vr(:), xr(:))
if (r2 < r2crit) then
lencounter(i,j) = .true.
loc_lvdotr(i,j) = (vdotr < 0.0_DP)
else
if (vdotr < 0.0_DP) then
v2 = dot_product(vr(:), vr(:))
tmin = -vdotr / v2
if (tmin < dt) then
r2min = r2 - vdotr * vdotr / v2
else
r2min = r2 + 2 * vdotr * dt + v2 * dt * dt
end if
r2min = min(r2min, r2)
if (r2min <= r2crit) then
lencounter(i,j) = .true.
loc_lvdotr(i,j) = (vdotr < 0.0_DP)
end if
end if
end if
end do
end do

lencounter = .false.
lany_encounter = any(lencounter(:,:))
if (lany_encounter) then
associate(pltpenc_list => system%pltpenc_list, nenc => system%pltpenc_list%nenc)
nenc_old = nenc
call pltpenc_list%resize(nenc_old + count(lencounter(:,:)))
pltpenc_list%status(nenc_old+1:nenc) = ACTIVE
pltpenc_list%level(nenc_old+1:nenc) = irec
pltpenc_list%lvdotr(nenc_old+1:nenc) = pack(loc_lvdotr(:,:), lencounter(:,:))
!*********************************************************************************************************
! This needs to be tested
pltpenc_list%index1(nenc_old+1:nenc) = pack(spread([(i, i = 1, npl)], dim=2, ncopies=ntp), lencounter(:,:))
pltpenc_list%index2(nenc_old+1:nenc) = pack(spread([(j, j = 1, ntp)], dim=1, ncopies=npl), lencounter(:,:))
!*********************************************************************************************************
select type(pl)
class is (symba_pl)
pl%lencounter(pltpenc_list%index1(nenc_old+1:nenc)) = .true.
end select
end associate
end if
end associate
return
end function symba_encounter_check_tp

Expand Down
11 changes: 5 additions & 6 deletions src/symba/symba_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,14 +16,14 @@ module subroutine symba_step_system(self, param, t, dt)
real(DP), intent(in) :: t !! Simulation time
real(DP), intent(in) :: dt !! Current stepsize
! Internals
logical :: lencounter_pl, lencounter_tp, lencounter
logical :: lencounter

call self%reset()
select type(pl => self%pl)
class is (symba_pl)
select type(tp => self%tp)
class is (symba_tp)
lencounter = pl%encounter_check(self, dt, 0) .or. tp%encounter_check(self, dt)
lencounter = pl%encounter_check(self, dt, 0) .or. tp%encounter_check(self, dt, 0)
if (lencounter) then
call self%interp(param, t, dt)
else
Expand Down Expand Up @@ -120,15 +120,14 @@ module recursive subroutine symba_step_recur_system(self, param, t, dt, ireci)
dtl = param%dt / (NTENC**ireci)
dth = 0.5_DP * dtl
IF (dtl / param%dt < VSMALL) THEN
WRITE(*, *) "SWIFTEST Warning:"
WRITE(*, *) " In symba_step_recur_system, local time step is too small"
WRITE(*, *) " Roundoff error will be important!"
write(*, *) "SWIFTEST Warning:"
write(*, *) " In symba_step_recur_system, local time step is too small"
write(*, *) " Roundoff error will be important!"
call util_exit(FAILURE)
END IF
irecp = ireci + 1
if (ireci == 0) then
icflg = 0

end if
end associate

Expand Down

0 comments on commit 7878958

Please sign in to comment.