From 8e6425c9b2b9263b4681b42b47504f11225ce948 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 12 Jan 2023 23:11:48 -0500 Subject: [PATCH] Fixing more array issues uncovered by gfortran --- src/encounter/encounter_util.f90 | 149 ++++++++++++++++--------------- src/swiftest/swiftest_io.f90 | 21 ++--- 2 files changed, 86 insertions(+), 84 deletions(-) diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 94ea7d942..d489ab2ff 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -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) @@ -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 @@ -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 diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index b58e3c20a..da0693b95 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -1828,8 +1828,6 @@ 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 @@ -1837,12 +1835,12 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i 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? @@ -2168,8 +2166,7 @@ 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 @@ -2177,7 +2174,7 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i 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 @@ -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.' @@ -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