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

Commit

Permalink
Fixing more array issues uncovered by gfortran
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Jan 13, 2023
1 parent 9a2399b commit 8e6425c
Show file tree
Hide file tree
Showing 2 changed files with 86 additions and 84 deletions.
149 changes: 77 additions & 72 deletions src/encounter/encounter_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,8 @@ module subroutine encounter_util_append_list(self, source, lsource_mask)
! Internals
integer(I4B) :: nold, nsrc

nold = self%nenc
nsrc = source%nenc
nold = int(self%nenc, kind=I4B)
nsrc = int(source%nenc, kind=I4B)
call swiftest_util_append(self%tcollision, source%tcollision, nold, nsrc, lsource_mask)
call swiftest_util_append(self%lclosest, source%lclosest, nold, nsrc, lsource_mask)
call swiftest_util_append(self%lvdotr, source%lvdotr, nold, nsrc, lsource_mask)
Expand Down Expand Up @@ -480,7 +480,8 @@ module subroutine encounter_util_snapshot(self, param, nbody_system, t, arg)
character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots)
! Arguments
class(encounter_snapshot), allocatable :: snapshot
integer(I4B) :: i, pi, pj, k, npl_snap, ntp_snap, iflag
integer(I4B) :: i, pii, pjj, npl_snap, ntp_snap, iflag
integer(I8B) :: k
real(DP), dimension(NDIM) :: rrel, vrel, rcom, vcom
real(DP) :: Gmtot, a, q, capm, tperi
real(DP), dimension(NDIM,2) :: rb,vb
Expand Down Expand Up @@ -597,79 +598,83 @@ module subroutine encounter_util_snapshot(self, param, nbody_system, t, arg)
call self%save(snapshot)
case("closest")
associate(plpl_encounter => nbody_system%plpl_encounter, pltp_encounter => nbody_system%pltp_encounter)
if (any(plpl_encounter%lclosest(:))) then
call pl_snap%setup(2, param)
do k = 1, plpl_encounter%nenc
if (plpl_encounter%lclosest(k)) then
pi = plpl_encounter%index1(k)
pj = plpl_encounter%index2(k)
select type(pl_snap)
class is (symba_pl)
select type(pl)
class is (symba_pl)
pl_snap%levelg(:) = pl%levelg([pi,pj])
end select
end select
pl_snap%id(:) = pl%id([pi,pj])
pl_snap%info(:) = pl%info([pi,pj])
pl_snap%Gmass(:) = pl%Gmass([pi,pj])
Gmtot = sum(pl_snap%Gmass(:))
if (param%lclose) pl_snap%radius(:) = pl%radius([pi,pj])
if (param%lrotation) then
do i = 1, NDIM
pl_snap%Ip(i,:) = pl%Ip(i,[pi,pj])
pl_snap%rot(i,:) = pl%rot(i,[pi,pj])
end do
if (plpl_encounter%nenc > 0) then
if (any(plpl_encounter%lclosest(:))) then
call pl_snap%setup(2, param)
do k = 1_I8B, plpl_encounter%nenc
if (plpl_encounter%lclosest(k)) then
pii = plpl_encounter%index1(k)
pjj = plpl_encounter%index2(k)
select type(pl_snap)
class is (symba_pl)
select type(pl)
class is (symba_pl)
pl_snap%levelg(:) = pl%levelg([pii,pjj])
end select
end select
pl_snap%id(:) = pl%id([pii,pjj])
pl_snap%info(:) = pl%info([pii,pjj])
pl_snap%Gmass(:) = pl%Gmass([pii,pjj])
Gmtot = sum(pl_snap%Gmass(:))
if (param%lclose) pl_snap%radius(:) = pl%radius([pii,pjj])
if (param%lrotation) then
do i = 1, NDIM
pl_snap%Ip(i,:) = pl%Ip(i,[pii,pjj])
pl_snap%rot(i,:) = pl%rot(i,[pii,pjj])
end do
end if

! Compute pericenter passage time to get the closest approach parameters
rrel(:) = plpl_encounter%r2(:,k) - plpl_encounter%r1(:,k)
vrel(:) = plpl_encounter%v2(:,k) - plpl_encounter%v1(:,k)
call swiftest_orbel_xv2aqt(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi)
snapshot%t = t + tperi
if ((snapshot%t < maxval(pl_snap%info(:)%origin_time)) .or. &
(snapshot%t > minval(pl_snap%info(:)%discard_time))) cycle

! Computer the center mass of the pair
rcom(:) = (plpl_encounter%r1(:,k) * pl_snap%Gmass(1) + plpl_encounter%r2(:,k) * pl_snap%Gmass(2)) / Gmtot
vcom(:) = (plpl_encounter%v1(:,k) * pl_snap%Gmass(1) + plpl_encounter%v2(:,k) * pl_snap%Gmass(2)) / Gmtot
rb(:,1) = plpl_encounter%r1(:,k) - rcom(:)
rb(:,2) = plpl_encounter%r2(:,k) - rcom(:)
vb(:,1) = plpl_encounter%v1(:,k) - vcom(:)
vb(:,2) = plpl_encounter%v2(:,k) - vcom(:)

! Drift the relative orbit to get the new relative position and velocity
call swiftest_drift_one(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), tperi, iflag)
if (iflag /= 0) write(*,*) "Danby error in encounter_util_snapshot_encounter. Closest approach positions and vectors may not be accurate."

! Get the new position and velocity vectors
rb(:,1) = -(pl_snap%Gmass(2) / Gmtot) * rrel(:)
rb(:,2) = (pl_snap%Gmass(1)) / Gmtot * rrel(:)

vb(:,1) = -(pl_snap%Gmass(2) / Gmtot) * vrel(:)
vb(:,2) = (pl_snap%Gmass(1)) / Gmtot * vrel(:)

! Move the CoM assuming constant velocity over the time it takes to reach periapsis
rcom(:) = rcom(:) + vcom(:) * tperi

! Compute the heliocentric position and velocity vector at periapsis
pl_snap%rh(:,1) = rb(:,1) + rcom(:)
pl_snap%rh(:,2) = rb(:,2) + rcom(:)
pl_snap%vh(:,1) = vb(:,1) + vcom(:)
pl_snap%vh(:,2) = vb(:,2) + vcom(:)

call pl_snap%sort("id", ascending=.true.)
call self%save(snapshot)
end if
end do

! Compute pericenter passage time to get the closest approach parameters
rrel(:) = plpl_encounter%r2(:,k) - plpl_encounter%r1(:,k)
vrel(:) = plpl_encounter%v2(:,k) - plpl_encounter%v1(:,k)
call swiftest_orbel_xv2aqt(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi)
snapshot%t = t + tperi
if ((snapshot%t < maxval(pl_snap%info(:)%origin_time)) .or. &
(snapshot%t > minval(pl_snap%info(:)%discard_time))) cycle

! Computer the center mass of the pair
rcom(:) = (plpl_encounter%r1(:,k) * pl_snap%Gmass(1) + plpl_encounter%r2(:,k) * pl_snap%Gmass(2)) / Gmtot
vcom(:) = (plpl_encounter%v1(:,k) * pl_snap%Gmass(1) + plpl_encounter%v2(:,k) * pl_snap%Gmass(2)) / Gmtot
rb(:,1) = plpl_encounter%r1(:,k) - rcom(:)
rb(:,2) = plpl_encounter%r2(:,k) - rcom(:)
vb(:,1) = plpl_encounter%v1(:,k) - vcom(:)
vb(:,2) = plpl_encounter%v2(:,k) - vcom(:)

! Drift the relative orbit to get the new relative position and velocity
call swiftest_drift_one(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), tperi, iflag)
if (iflag /= 0) write(*,*) "Danby error in encounter_util_snapshot_encounter. Closest approach positions and vectors may not be accurate."

! Get the new position and velocity vectors
rb(:,1) = -(pl_snap%Gmass(2) / Gmtot) * rrel(:)
rb(:,2) = (pl_snap%Gmass(1)) / Gmtot * rrel(:)

vb(:,1) = -(pl_snap%Gmass(2) / Gmtot) * vrel(:)
vb(:,2) = (pl_snap%Gmass(1)) / Gmtot * vrel(:)

! Move the CoM assuming constant velocity over the time it takes to reach periapsis
rcom(:) = rcom(:) + vcom(:) * tperi

! Compute the heliocentric position and velocity vector at periapsis
pl_snap%rh(:,1) = rb(:,1) + rcom(:)
pl_snap%rh(:,2) = rb(:,2) + rcom(:)
pl_snap%vh(:,1) = vb(:,1) + vcom(:)
pl_snap%vh(:,2) = vb(:,2) + vcom(:)

call pl_snap%sort("id", ascending=.true.)
call self%save(snapshot)
end if
end do

plpl_encounter%lclosest(:) = .false.
plpl_encounter%lclosest(:) = .false.
end if
end if

if (any(pltp_encounter%lclosest(:))) then
do k = 1, pltp_encounter%nenc
end do
pltp_encounter%lclosest(:) = .false.
if (pltp_encounter%nenc > 0) then
if (any(pltp_encounter%lclosest(:))) then
do k = 1_I8B, pltp_encounter%nenc
end do
pltp_encounter%lclosest(:) = .false.
end if
end if
end associate
case default
Expand Down
21 changes: 9 additions & 12 deletions src/swiftest/swiftest_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1828,21 +1828,19 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i
!! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
!!
!! Read in parameters for the integration
!! Currently this procedure does not work in user-defined derived-type input mode
!! e.g. read(unit,'(DT)') param
!! as the newline characters are ignored in the input file when compiled in ifort.
!!
!! Adapted from David E. Kaufmann's Swifter routine io_init_param.f90
!! Adapted from Martin Duncan's Swift routine io_init_param.f
implicit none
! Arguments
class(swiftest_parameters), intent(inout) :: self !! Collection of parameters
integer, intent(in) :: unit !! File unit number
character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT.
integer(I4B), intent(in) :: unit !! File unit number
character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT.
!! If you do not include a char-literal-constant, the iotype argument contains only DT.
character(len=*), intent(in) :: v_list(:) !! The first element passes the integrator code to the reader
integer, intent(out) :: iostat !! IO status code
character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0
character(len=*), intent(in) :: v_list(:) !! The first element passes the integrator code to the reader
integer(I4B), intent(out) :: iostat !! IO status code
character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0
! Internals
logical :: tstart_set = .false. !! Is the final time set in the input file?
logical :: tstop_set = .false. !! Is the final time set in the input file?
Expand Down Expand Up @@ -2168,16 +2166,15 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i
param%lenc_save_trajectory = (param%encounter_save == "TRAJECTORY") .or. (param%encounter_save == "BOTH")
param%lenc_save_closest = (param%encounter_save == "CLOSEST") .or. (param%encounter_save == "BOTH")

integrator = v_list(1)
if ((integrator == INT_RMVS) .or. (integrator == INT_SYMBA)) then
if ((param%integrator == INT_RMVS) .or. (param%integrator == INT_SYMBA)) then
if (.not.param%lclose) then
write(iomsg,*) 'This integrator requires CHK_CLOSE to be enabled.'
iostat = -1
return
end if
end if

param%lmtiny_pl = (integrator == INT_SYMBA)
param%lmtiny_pl = (param%integrator == INT_SYMBA)

if (param%lmtiny_pl .and. param%GMTINY < 0.0_DP) then
write(iomsg,*) "GMTINY invalid or not set: ", param%GMTINY
Expand All @@ -2204,7 +2201,7 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i
end if

! Determine if the GR flag is set correctly for this integrator
select case(integrator)
select case(param%integrator)
case(INT_WHM, INT_RMVS, INT_HELIO, INT_SYMBA)
case default
if (param%lgr) write(iomsg, *) 'GR is not yet implemented for this integrator. This parameter will be ignored.'
Expand Down Expand Up @@ -2867,7 +2864,7 @@ module subroutine swiftest_io_read_in_param(self, param_file_name)
!! as the newline characters are ignored in the input file when compiled in ifort.

!read(LUN,'(DT)', iostat= ierr, iomsg = errmsg) self
call self%reader(LUN, iotype= "none", v_list = [self%integrator], iostat = ierr, iomsg = errmsg)
call self%reader(LUN, iotype= "none", v_list=[""], iostat = ierr, iomsg = errmsg)
if (ierr == 0) return

667 continue
Expand Down

0 comments on commit 8e6425c

Please sign in to comment.