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

Commit

Permalink
Reverted coarray branch source code in order to take a different appr…
Browse files Browse the repository at this point in the history 
…oach
  • Loading branch information
@daminton
daminton committed Mar 21, 2023
1 parent 1963fa7 commit 622502d
Show file tree
Hide file tree
Showing 8 changed files with 265 additions and 283 deletions.
2 changes: 0 additions & 2 deletions src/collision/collision_regime.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -249,8 +249,6 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2,
ke = 0.5_DP * Vimp**2
pe = - GC * m1 * m2 / (Mtot * norm2(rh2 - rh1))

Qmerge = Gc * m1 * m2 / (.mag.(rh2 - rh1)) + 3*m1 / (5*rad1) + 3*m2 / (5*rad1) - U_binding ! Change in energy due to a pure merger

if ((m1 < min_mfrag).or.(m2 < min_mfrag)) then
regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime
call swiftest_io_log_one_message(COLLISION_LOG_OUT, &
Expand Down
2 changes: 1 addition & 1 deletion src/misc/solver_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ module solver
use lambda_function
use, intrinsic :: ieee_exceptions
private
public :: solve_linear_system, solve_roots
public :: solve_linear_system, solve_rkf45, solve_roots

interface solve_linear_system
module procedure solve_linear_system_dp
Expand Down
2 changes: 1 addition & 1 deletion src/rmvs/rmvs_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ module rmvs
integer(I4B), dimension(:), allocatable :: plind !! Connects the planetocentric indices back to the heliocentric planet list
type(rmvs_interp), dimension(:), allocatable :: outer !! interpolated heliocentric central body position for outer encounters
type(rmvs_interp), dimension(:), allocatable :: inner !! interpolated heliocentric central body position for inner encounters
class(base_nbody_system), dimension(:), allocatable :: planetocentric !! Planetocentric version of the massive body objects (one for each massive body)
class(rmvs_nbody_system), dimension(:), allocatable :: planetocentric !! Planetocentric version of the massive body objects (one for each massive body)
logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations
contains
procedure :: setup => rmvs_util_setup_pl !! Constructor method - Allocates space for the input number of bodiess
Expand Down
186 changes: 90 additions & 96 deletions src/rmvs/rmvs_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -451,67 +451,64 @@ subroutine rmvs_make_planetocentric(param, cb, pl, tp)
if (allocated(encmask)) deallocate(encmask)
allocate(encmask(ntp))
encmask(1:ntp) = tp%plencP(1:ntp) == i
select type (planetocentric => pl%planetocentric)
class is(rmvs_nbody_system)
allocate(rmvs_tp :: planetocentric(i)%tp)
! Create encountering test particle structure
select type(cbenci => planetocentric(i)%cb)
class is (rmvs_cb)
select type(plenci => planetocentric(i)%pl)
class is (rmvs_pl)
select type(tpenci => planetocentric(i)%tp)
class is (rmvs_tp)
tpenci%lplanetocentric = .true.
associate(nenci => pl%nenc(i))
call tpenci%setup(nenci, param)
tpenci%cb_heliocentric = cb
tpenci%ipleP = i
tpenci%lmask(1:nenci) = .true.
tpenci%status(1:nenci) = ACTIVE
! Grab all the encountering test particles and convert them to a planetocentric frame
tpenci%id(1:nenci) = pack(tp%id(1:ntp), encmask(1:ntp))
do j = 1, NDIM
tpenci%rheliocentric(j, 1:nenci) = pack(tp%rh(j,1:ntp), encmask(:))
tpenci%rh(j, 1:nenci) = tpenci%rheliocentric(j, 1:nenci) - pl%inner(0)%x(j, i)
tpenci%vh(j, 1:nenci) = pack(tp%vh(j, 1:ntp), encmask(1:ntp)) - pl%inner(0)%v(j, i)
end do
tpenci%lperi(1:nenci) = pack(tp%lperi(1:ntp), encmask(1:ntp))
tpenci%plperP(1:nenci) = pack(tp%plperP(1:ntp), encmask(1:ntp))
! Make sure that the test particles get the planetocentric value of mu
allocate(cbenci%inner(0:NTPHENC))
do inner_index = 0, NTPHENC
allocate(plenci%inner(inner_index)%x, mold=pl%inner(inner_index)%x)
allocate(plenci%inner(inner_index)%v, mold=pl%inner(inner_index)%x)
allocate(cbenci%inner(inner_index)%x(NDIM,1))
allocate(cbenci%inner(inner_index)%v(NDIM,1))
cbenci%inner(inner_index)%x(:,1) = pl%inner(inner_index)%x(:, i)
cbenci%inner(inner_index)%v(:,1) = pl%inner(inner_index)%v(:, i)
plenci%inner(inner_index)%x(:,1) = -cbenci%inner(inner_index)%x(:,1)
plenci%inner(inner_index)%v(:,1) = -cbenci%inner(inner_index)%v(:,1)

if (param%loblatecb) then
allocate(plenci%inner(inner_index)%aobl, mold=pl%inner(inner_index)%aobl)
allocate(cbenci%inner(inner_index)%aobl(NDIM,1))
cbenci%inner(inner_index)%aobl(:,1) = pl%inner(inner_index)%aobl(:, i)
end if

if (param%ltides) then
allocate(plenci%inner(inner_index)%atide, mold=pl%inner(inner_index)%atide)
allocate(cbenci%inner(inner_index)%atide(NDIM,1))
cbenci%inner(inner_index)%atide(:,1) = pl%inner(inner_index)%atide(:, i)
end if

do j = 2, npl
ipc2hc = plenci%plind(j)
plenci%inner(inner_index)%x(:,j) = pl%inner(inner_index)%x(:, ipc2hc) &
- cbenci%inner(inner_index)%x(:,1)
plenci%inner(inner_index)%v(:,j) = pl%inner(inner_index)%v(:, ipc2hc) &
- cbenci%inner(inner_index)%v(:,1)
end do
allocate(rmvs_tp :: pl%planetocentric(i)%tp)
! Create encountering test particle structure
select type(cbenci => pl%planetocentric(i)%cb)
class is (rmvs_cb)
select type(plenci => pl%planetocentric(i)%pl)
class is (rmvs_pl)
select type(tpenci => pl%planetocentric(i)%tp)
class is (rmvs_tp)
tpenci%lplanetocentric = .true.
associate(nenci => pl%nenc(i))
call tpenci%setup(nenci, param)
tpenci%cb_heliocentric = cb
tpenci%ipleP = i
tpenci%lmask(1:nenci) = .true.
tpenci%status(1:nenci) = ACTIVE
! Grab all the encountering test particles and convert them to a planetocentric frame
tpenci%id(1:nenci) = pack(tp%id(1:ntp), encmask(1:ntp))
do j = 1, NDIM
tpenci%rheliocentric(j, 1:nenci) = pack(tp%rh(j,1:ntp), encmask(:))
tpenci%rh(j, 1:nenci) = tpenci%rheliocentric(j, 1:nenci) - pl%inner(0)%x(j, i)
tpenci%vh(j, 1:nenci) = pack(tp%vh(j, 1:ntp), encmask(1:ntp)) - pl%inner(0)%v(j, i)
end do
tpenci%lperi(1:nenci) = pack(tp%lperi(1:ntp), encmask(1:ntp))
tpenci%plperP(1:nenci) = pack(tp%plperP(1:ntp), encmask(1:ntp))
! Make sure that the test particles get the planetocentric value of mu
allocate(cbenci%inner(0:NTPHENC))
do inner_index = 0, NTPHENC
allocate(plenci%inner(inner_index)%x, mold=pl%inner(inner_index)%x)
allocate(plenci%inner(inner_index)%v, mold=pl%inner(inner_index)%x)
allocate(cbenci%inner(inner_index)%x(NDIM,1))
allocate(cbenci%inner(inner_index)%v(NDIM,1))
cbenci%inner(inner_index)%x(:,1) = pl%inner(inner_index)%x(:, i)
cbenci%inner(inner_index)%v(:,1) = pl%inner(inner_index)%v(:, i)
plenci%inner(inner_index)%x(:,1) = -cbenci%inner(inner_index)%x(:,1)
plenci%inner(inner_index)%v(:,1) = -cbenci%inner(inner_index)%v(:,1)

if (param%loblatecb) then
allocate(plenci%inner(inner_index)%aobl, mold=pl%inner(inner_index)%aobl)
allocate(cbenci%inner(inner_index)%aobl(NDIM,1))
cbenci%inner(inner_index)%aobl(:,1) = pl%inner(inner_index)%aobl(:, i)
end if

if (param%ltides) then
allocate(plenci%inner(inner_index)%atide, mold=pl%inner(inner_index)%atide)
allocate(cbenci%inner(inner_index)%atide(NDIM,1))
cbenci%inner(inner_index)%atide(:,1) = pl%inner(inner_index)%atide(:, i)
end if

do j = 2, npl
ipc2hc = plenci%plind(j)
plenci%inner(inner_index)%x(:,j) = pl%inner(inner_index)%x(:, ipc2hc) &
- cbenci%inner(inner_index)%x(:,1)
plenci%inner(inner_index)%v(:,j) = pl%inner(inner_index)%v(:, ipc2hc) &
- cbenci%inner(inner_index)%v(:,1)
end do
call tpenci%set_mu(cbenci)
end associate
end select
end do
call tpenci%set_mu(cbenci)
end associate
end select
end select
end select
Expand Down Expand Up @@ -614,42 +611,39 @@ subroutine rmvs_end_planetocentric(pl, tp)
associate (npl => pl%nbody, ntp => tp%nbody)
do i = 1, npl
if (pl%nenc(i) == 0) cycle
select type(planetocentric => pl%planetocentric)
class is(rmvs_nbody_system)
select type(cbenci => planetocentric(i)%cb)
class is (rmvs_cb)
select type(plenci => planetocentric(i)%pl)
class is (rmvs_pl)
select type(tpenci => planetocentric(i)%tp)
class is (rmvs_tp)
associate(nenci => pl%nenc(i))
if (allocated(tpind)) deallocate(tpind)
allocate(tpind(nenci))
! Index array of encountering test particles
if (allocated(encmask)) deallocate(encmask)
allocate(encmask(ntp))
encmask(1:ntp) = tp%plencP(1:ntp) == i
tpind(1:nenci) = pack([(j,j=1,ntp)], encmask(1:ntp))

! Copy the results of the integration back over and shift back to heliocentric reference
tp%status(tpind(1:nenci)) = tpenci%status(1:nenci)
tp%lmask(tpind(1:nenci)) = tpenci%lmask(1:nenci)
do j = 1, NDIM
tp%rh(j, tpind(1:nenci)) = tpenci%rh(j,1:nenci) + pl%inner(NTPHENC)%x(j, i)
tp%vh(j, tpind(1:nenci)) = tpenci%vh(j,1:nenci) + pl%inner(NTPHENC)%v(j, i)
end do
tp%lperi(tpind(1:nenci)) = tpenci%lperi(1:nenci)
tp%plperP(tpind(1:nenci)) = tpenci%plperP(1:nenci)
deallocate(planetocentric(i)%tp)
deallocate(cbenci%inner)
do inner_index = 0, NTPHENC
deallocate(plenci%inner(inner_index)%x)
deallocate(plenci%inner(inner_index)%v)
if (allocated(plenci%inner(inner_index)%aobl)) deallocate(plenci%inner(inner_index)%aobl)
if (allocated(plenci%inner(inner_index)%atide)) deallocate(plenci%inner(inner_index)%atide)
end do
end associate
end select
select type(cbenci => pl%planetocentric(i)%cb)
class is (rmvs_cb)
select type(plenci => pl%planetocentric(i)%pl)
class is (rmvs_pl)
select type(tpenci => pl%planetocentric(i)%tp)
class is (rmvs_tp)
associate(nenci => pl%nenc(i))
if (allocated(tpind)) deallocate(tpind)
allocate(tpind(nenci))
! Index array of encountering test particles
if (allocated(encmask)) deallocate(encmask)
allocate(encmask(ntp))
encmask(1:ntp) = tp%plencP(1:ntp) == i
tpind(1:nenci) = pack([(j,j=1,ntp)], encmask(1:ntp))

! Copy the results of the integration back over and shift back to heliocentric reference
tp%status(tpind(1:nenci)) = tpenci%status(1:nenci)
tp%lmask(tpind(1:nenci)) = tpenci%lmask(1:nenci)
do j = 1, NDIM
tp%rh(j, tpind(1:nenci)) = tpenci%rh(j,1:nenci) + pl%inner(NTPHENC)%x(j, i)
tp%vh(j, tpind(1:nenci)) = tpenci%vh(j,1:nenci) + pl%inner(NTPHENC)%v(j, i)
end do
tp%lperi(tpind(1:nenci)) = tpenci%lperi(1:nenci)
tp%plperP(tpind(1:nenci)) = tpenci%plperP(1:nenci)
deallocate(pl%planetocentric(i)%tp)
deallocate(cbenci%inner)
do inner_index = 0, NTPHENC
deallocate(plenci%inner(inner_index)%x)
deallocate(plenci%inner(inner_index)%v)
if (allocated(plenci%inner(inner_index)%aobl)) deallocate(plenci%inner(inner_index)%aobl)
if (allocated(plenci%inner(inner_index)%atide)) deallocate(plenci%inner(inner_index)%atide)
end do
end associate
end select
end select
end select
Expand Down
67 changes: 32 additions & 35 deletions src/rmvs/rmvs_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -356,42 +356,39 @@ module subroutine rmvs_util_setup_initialize_system(self, param)
tp%lplanetocentric = .false.
cb%lplanetocentric = .false.
associate(npl => pl%nbody)
allocate(rmvs_nbody_system :: pl%planetocentric(npl))
select type(planetocentric => pl%planetocentric)
class is (rmvs_nbody_system)
planetocentric(:)%lplanetocentric = .true.
do i = 1, npl
allocate(planetocentric(i)%cb, source=cb)
allocate(rmvs_pl :: planetocentric(i)%pl)
select type(cbenci => planetocentric(i)%cb)
class is (rmvs_cb)
select type(plenci => planetocentric(i)%pl)
class is (rmvs_pl)
cbenci%lplanetocentric = .true.
plenci%lplanetocentric = .true.
call plenci%setup(npl, param)
plenci%status(:) = ACTIVE
plenci%lmask(:) = .true.
! plind stores the heliocentric index value of a planetocentric planet
! e.g. Consider an encounter with planet 3.
! Then the following will be the values of plind:
! pl%planetocentric(3)%pl%plind(1) = 0 (central body - never used)
! pl%planetocentric(3)%pl%plind(2) = 1
! pl%planetocentric(3)%pl%plind(3) = 2
! pl%planetocentric(3)%pl%plind(4) = 4
! pl%planetocentric(3)%pl%plind(5) = 5
! etc.
allocate(plenci%plind(npl))
plenci%plind(1:npl) = [(j,j=1,npl)]
plenci%plind(2:npl) = pack(plenci%plind(1:npl), plenci%plind(1:npl) /= i)
plenci%plind(1) = 0
plenci%Gmass(1) = cb%Gmass
plenci%Gmass(2:npl) = pl%Gmass(plenci%plind(2:npl))
cbenci%Gmass = pl%Gmass(i)
end select
allocate(pl%planetocentric(npl))
pl%planetocentric(:)%lplanetocentric = .true.
do i = 1, npl
allocate(pl%planetocentric(i)%cb, source=cb)
allocate(rmvs_pl :: pl%planetocentric(i)%pl)
select type(cbenci => pl%planetocentric(i)%cb)
class is (rmvs_cb)
select type(plenci => pl%planetocentric(i)%pl)
class is (rmvs_pl)
cbenci%lplanetocentric = .true.
plenci%lplanetocentric = .true.
call plenci%setup(npl, param)
plenci%status(:) = ACTIVE
plenci%lmask(:) = .true.
! plind stores the heliocentric index value of a planetocentric planet
! e.g. Consider an encounter with planet 3.
! Then the following will be the values of plind:
! pl%planetocentric(3)%pl%plind(1) = 0 (central body - never used)
! pl%planetocentric(3)%pl%plind(2) = 1
! pl%planetocentric(3)%pl%plind(3) = 2
! pl%planetocentric(3)%pl%plind(4) = 4
! pl%planetocentric(3)%pl%plind(5) = 5
! etc.
allocate(plenci%plind(npl))
plenci%plind(1:npl) = [(j,j=1,npl)]
plenci%plind(2:npl) = pack(plenci%plind(1:npl), plenci%plind(1:npl) /= i)
plenci%plind(1) = 0
plenci%Gmass(1) = cb%Gmass
plenci%Gmass(2:npl) = pl%Gmass(plenci%plind(2:npl))
cbenci%Gmass = pl%Gmass(i)
end select
end do
end select
end select
end do
end associate
end select
end select
Expand Down
Loading

0 comments on commit 622502d

Please sign in to comment.