From 622502d1dc3b30f1b343697ef607a07aae54b248 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 21 Mar 2023 14:04:43 -0400 Subject: [PATCH] Reverted coarray branch source code in order to take a different approach --- src/collision/collision_regime.f90 | 2 - src/misc/solver_module.f90 | 2 +- src/rmvs/rmvs_module.f90 | 2 +- src/rmvs/rmvs_step.f90 | 186 ++++++++++++++--------------- src/rmvs/rmvs_util.f90 | 67 +++++------ src/swiftest/swiftest_driver.f90 | 133 ++++++++++----------- src/swiftest/swiftest_module.f90 | 29 +++-- src/swiftest/swiftest_util.f90 | 127 ++++++++++---------- 8 files changed, 265 insertions(+), 283 deletions(-) diff --git a/src/collision/collision_regime.f90 b/src/collision/collision_regime.f90 index 5bb77099a..7480c6a0d 100644 --- a/src/collision/collision_regime.f90 +++ b/src/collision/collision_regime.f90 @@ -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, & diff --git a/src/misc/solver_module.f90 b/src/misc/solver_module.f90 index 504e9215e..410116b2a 100644 --- a/src/misc/solver_module.f90 +++ b/src/misc/solver_module.f90 @@ -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 diff --git a/src/rmvs/rmvs_module.f90 b/src/rmvs/rmvs_module.f90 index 74910d7d3..733a81e60 100644 --- a/src/rmvs/rmvs_module.f90 +++ b/src/rmvs/rmvs_module.f90 @@ -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 diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 18f9bd644..bcb38797a 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -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 @@ -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 diff --git a/src/rmvs/rmvs_util.f90 b/src/rmvs/rmvs_util.f90 index d4d486da7..b9462840b 100644 --- a/src/rmvs/rmvs_util.f90 +++ b/src/rmvs/rmvs_util.f90 @@ -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 diff --git a/src/swiftest/swiftest_driver.f90 b/src/swiftest/swiftest_driver.f90 index efc50e7c0..bfa950222 100644 --- a/src/swiftest/swiftest_driver.f90 +++ b/src/swiftest/swiftest_driver.f90 @@ -19,7 +19,7 @@ program swiftest_driver use swiftest implicit none - class(base_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated + class(swiftest_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated class(swiftest_parameters), allocatable :: param !! Run configuration parameters character(len=:), allocatable :: integrator !! Integrator type code (see globals for symbolic names) character(len=:), allocatable :: param_file_name !! Name of the file containing user-defined parameters @@ -70,78 +70,75 @@ program swiftest_driver ! Construct the main n-body nbody_system using the user-input integrator to choose the type of nbody_system call swiftest_util_setup_construct_system(nbody_system, param) - select type(nbody_system) - class is (swiftest_nbody_system) - - !> Define the maximum number of threads - nthreads = 1 ! In the *serial* case - !$ nthreads = omp_get_max_threads() ! In the *parallel* case - !$ write(param%display_unit,'(a)') ' OpenMP parameters:' - !$ write(param%display_unit,'(a)') ' ------------------' - !$ write(param%display_unit,'(a,i3,/)') ' Number of threads = ', nthreads - !$ if (param%log_output) write(*,'(a,i3)') ' OpenMP: Number of threads = ',nthreads - - call nbody_system%initialize(param) - - associate (system_history => nbody_system%system_history) - ! If this is a new run, compute energy initial conditions (if energy tracking is turned on) and write the initial conditions to file. - call nbody_system%display_run_information(param, integration_timer, phase="first") - if (param%lenergy) then - if (param%lrestart) then - call nbody_system%get_t0_values(param) - else - call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum - end if - call nbody_system%conservation_report(param, lterminal=.true.) + + !> Define the maximum number of threads + nthreads = 1 ! In the *serial* case + !$ nthreads = omp_get_max_threads() ! In the *parallel* case + !$ write(param%display_unit,'(a)') ' OpenMP parameters:' + !$ write(param%display_unit,'(a)') ' ------------------' + !$ write(param%display_unit,'(a,i3,/)') ' Number of threads = ', nthreads + !$ if (param%log_output) write(*,'(a,i3)') ' OpenMP: Number of threads = ',nthreads + + call nbody_system%initialize(param) + + associate (system_history => nbody_system%system_history) + ! If this is a new run, compute energy initial conditions (if energy tracking is turned on) and write the initial conditions to file. + call nbody_system%display_run_information(param, integration_timer, phase="first") + if (param%lenergy) then + if (param%lrestart) then + call nbody_system%get_t0_values(param) + else + call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum end if - call system_history%take_snapshot(param,nbody_system) - call nbody_system%dump(param) - - do iloop = istart, nloops - !> Step the nbody_system forward in time - call integration_timer%start() - call nbody_system%step(param, nbody_system%t, dt) - call integration_timer%stop() - - nbody_system%t = t0 + iloop * dt - - !> Evaluate any discards or collisional outcomes - call nbody_system%discard(param) - - !> If the loop counter is at the output cadence value, append the data file with a single frame - if (istep_out > 0) then - iout = iout + 1 - if ((iout == istep) .or. (iloop == nloops)) then - iout = 0 - idump = idump + 1 - if (ltstretch) then - nout = nout + 1 - istep = floor(istep_out * fstep_out**nout, kind=I4B) - end if - - call system_history%take_snapshot(param,nbody_system) - - if (idump == dump_cadence) then - idump = 0 - call nbody_system%dump(param) - end if - - call integration_timer%report(message="Integration steps:", unit=display_unit) - call nbody_system%display_run_information(param, integration_timer) - call integration_timer%reset() - if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) + call nbody_system%conservation_report(param, lterminal=.true.) + end if + call system_history%take_snapshot(param,nbody_system) + call nbody_system%dump(param) + + do iloop = istart, nloops + !> Step the nbody_system forward in time + call integration_timer%start() + call nbody_system%step(param, nbody_system%t, dt) + call integration_timer%stop() + + nbody_system%t = t0 + iloop * dt + + !> Evaluate any discards or collisional outcomes + call nbody_system%discard(param) + + !> If the loop counter is at the output cadence value, append the data file with a single frame + if (istep_out > 0) then + iout = iout + 1 + if ((iout == istep) .or. (iloop == nloops)) then + iout = 0 + idump = idump + 1 + if (ltstretch) then + nout = nout + 1 + istep = floor(istep_out * fstep_out**nout, kind=I4B) + end if + + call system_history%take_snapshot(param,nbody_system) + if (idump == dump_cadence) then + idump = 0 + call nbody_system%dump(param) end if + + call integration_timer%report(message="Integration steps:", unit=display_unit) + call nbody_system%display_run_information(param, integration_timer) + call integration_timer%reset() + if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) + end if + end if - end do - ! Dump any remaining history if it exists - call nbody_system%dump(param) - call system_history%dump(param) - call nbody_system%display_run_information(param, integration_timer, phase="last") - - end associate - end select + end do + ! Dump any remaining history if it exists + call nbody_system%dump(param) + call system_history%dump(param) + call nbody_system%display_run_information(param, integration_timer, phase="last") + + end associate end associate call base_util_exit(SUCCESS) diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 00c5ba228..4ac6cb3c3 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -304,13 +304,14 @@ module swiftest !> An abstract class for a basic Swiftest nbody system - type, extends(base_nbody_system) :: swiftest_nbody_system + type, abstract, extends(base_nbody_system) :: swiftest_nbody_system !! This superclass contains a minimial nbody_system of a set of test particles (tp), massive bodies (pl), and a central body (cb) !! The full swiftest_nbody_system type that is used as the parent class of all integrators is defined in collision class(swiftest_cb), allocatable :: cb !! Central body data structure class(swiftest_pl), allocatable :: pl !! Massive body data structure class(swiftest_tp), allocatable :: tp !! Test particle data structure + class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure class(swiftest_pl), allocatable :: pl_discards !! Discarded massive body particle data structure class(swiftest_pl), allocatable :: pl_adds !! List of added bodies in mergers or collisions @@ -373,7 +374,7 @@ module swiftest !! the test particles contains !> Each integrator will have its own version of the step - + procedure(abstract_step_system), deferred :: step ! Concrete classes that are common to the basic integrator (only test particles considered for discard) procedure :: discard => swiftest_discard_system !! Perform a discard step on the nbody_system @@ -390,8 +391,7 @@ module swiftest procedure :: read_in => swiftest_io_read_in_system !! Reads the initial conditions for an nbody system procedure :: write_frame_system => swiftest_io_write_frame_system !! Write a frame of input data from file procedure :: obl_pot => swiftest_obl_pot_system !! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body - procedure :: step => swiftest_step_system !! Placeholder step method. Each integrator will override - procedure :: dealloc => swiftest_util_dealloc_system !! Deallocates all allocatables and resets all values to defaults. Acts as a base for a finalizer + procedure :: dealloc => swiftest_util_dealloc_system !! Deallocates all allocatables and resets all values to defaults. Acts as a base for a finalizer procedure :: get_energy_and_momentum => swiftest_util_get_energy_and_momentum_system !! Calculates the total nbody_system energy and momentum procedure :: get_idvals => swiftest_util_get_idvalues_system !! Returns an array of all id values in use in the nbody_system procedure :: rescale => swiftest_util_rescale_system !! Rescales the nbody_system into a new set of units @@ -448,6 +448,15 @@ subroutine abstract_step_body(self, nbody_system, param, t, dt) real(DP), intent(in) :: t !! Simulation time real(DP), intent(in) :: dt !! Current stepsize end subroutine abstract_step_body + + subroutine abstract_step_system(self, param, t, dt) + import DP, swiftest_nbody_system, swiftest_parameters + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody_system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Simulation time + real(DP), intent(in) :: dt !! Current stepsize + end subroutine abstract_step_system end interface @@ -1049,7 +1058,7 @@ end subroutine swiftest_util_setup_body module subroutine swiftest_util_setup_construct_system(nbody_system, param) implicit none - class(base_nbody_system), allocatable, intent(inout) :: nbody_system !! Swiftest nbody_system object + class(swiftest_nbody_system), allocatable, intent(inout) :: nbody_system !! Swiftest nbody_system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine swiftest_util_setup_construct_system @@ -1943,16 +1952,6 @@ subroutine swiftest_final_kin(self) end subroutine swiftest_final_kin - subroutine swiftest_step_system(self, param, t, dt) - implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody_system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Simulation time - real(DP), intent(in) :: dt !! Current stepsize - return - end subroutine swiftest_step_system - - subroutine swiftest_final_storage(self) !! author: David A. Minton !! diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index bf6a30087..fcf130f83 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -2703,8 +2703,8 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) !! implicit none ! Arguments - class(base_nbody_system), allocatable, intent(inout) :: nbody_system !! Swiftest nbody_system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_nbody_system), allocatable, intent(inout) :: nbody_system !! Swiftest nbody_system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters select case(param%integrator) case (INT_BS) @@ -2767,11 +2767,11 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) call base_util_exit(FAILURE) end select - select type(nbody_system) - class is (swiftest_nbody_system) - allocate(swiftest_particle_info :: nbody_system%cb%info) - nbody_system%t = param%tstart - end select + allocate(swiftest_particle_info :: nbody_system%cb%info) + + + nbody_system%t = param%tstart + return end subroutine swiftest_util_setup_construct_system @@ -3055,7 +3055,7 @@ module subroutine swiftest_util_snapshot_system(self, param, nbody_system, t, ar real(DP), intent(in), optional :: t !! Time of snapshot if different from nbody_system time character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) ! Internals - class(base_nbody_system), allocatable :: snapshot + class(swiftest_nbody_system), allocatable :: snapshot ! To allow for runs to be restarted in a bit-identical way, we'll need to run the same coordinate conversion routines we would run upon restarting select type(pl => nbody_system%pl) @@ -3077,63 +3077,60 @@ module subroutine swiftest_util_snapshot_system(self, param, nbody_system, t, ar ! Take a minimal snapshot wihout all of the extra storage objects allocate(snapshot, mold=nbody_system) - select type(snapshot) - class is (swiftest_nbody_system) - allocate(snapshot%cb, source=nbody_system%cb ) - allocate(snapshot%pl, source=nbody_system%pl ) - allocate(snapshot%tp, source=nbody_system%tp ) - allocate(snapshot%system_history) - allocate(snapshot%system_history%nc, source=nbody_system%system_history%nc) - snapshot%system_history%nc%lfile_is_open = .true. - - snapshot%t = nbody_system%t - snapshot%GMtot = nbody_system%GMtot - snapshot%ke_orbit = nbody_system%ke_orbit - snapshot%ke_spin = nbody_system%ke_spin - snapshot%pe = nbody_system%pe - snapshot%be = nbody_system%be - snapshot%te = nbody_system%te - snapshot%oblpot = nbody_system%oblpot - snapshot%L_orbit = nbody_system%L_orbit - snapshot%L_spin = nbody_system%L_spin - snapshot%L_total = nbody_system%L_total - snapshot%ke_orbit_orig = nbody_system%ke_orbit_orig - snapshot%ke_spin_orig = nbody_system%ke_spin_orig - snapshot%pe_orig = nbody_system%pe_orig - snapshot%be_orig = nbody_system%be_orig - snapshot%E_orbit_orig = nbody_system%E_orbit_orig - snapshot%GMtot_orig = nbody_system%GMtot_orig - snapshot%L_total_orig = nbody_system%L_total_orig - snapshot%L_orbit_orig = nbody_system%L_orbit_orig - snapshot%L_spin_orig = nbody_system%L_spin_orig - snapshot%L_escape = nbody_system%L_escape - snapshot%GMescape = nbody_system%GMescape - snapshot%E_collisions = nbody_system%E_collisions - snapshot%E_untracked = nbody_system%E_untracked - snapshot%ke_orbit_error = nbody_system%ke_orbit_error - snapshot%ke_spin_error = nbody_system%ke_spin_error - snapshot%pe_error = nbody_system%pe_error - snapshot%be_error = nbody_system%be_error - snapshot%E_orbit_error = nbody_system%E_orbit_error - snapshot%Ecoll_error = nbody_system%Ecoll_error - snapshot%E_untracked_error = nbody_system%E_untracked_error - snapshot%te_error = nbody_system%te_error - snapshot%L_orbit_error = nbody_system%L_orbit_error - snapshot%L_spin_error = nbody_system%L_spin_error - snapshot%L_escape_error = nbody_system%L_escape_error - snapshot%L_total_error = nbody_system%L_total_error - snapshot%Mtot_error = nbody_system%Mtot_error - snapshot%Mescape_error = nbody_system%Mescape_error - snapshot%lbeg = nbody_system%lbeg - - - ! Store a snapshot of the nbody_system for posterity - call base_util_snapshot_save(self, snapshot) - self%nt = self%iframe - self%nid = self%nid + 1 ! Central body - if (allocated(nbody_system%pl)) self%nid = self%nid + nbody_system%pl%nbody - if (allocated(nbody_system%tp)) self%nid = self%nid + nbody_system%tp%nbody - end select + allocate(snapshot%cb, source=nbody_system%cb ) + allocate(snapshot%pl, source=nbody_system%pl ) + allocate(snapshot%tp, source=nbody_system%tp ) + allocate(snapshot%system_history) + allocate(snapshot%system_history%nc, source=nbody_system%system_history%nc) + snapshot%system_history%nc%lfile_is_open = .true. + + snapshot%t = nbody_system%t + snapshot%GMtot = nbody_system%GMtot + snapshot%ke_orbit = nbody_system%ke_orbit + snapshot%ke_spin = nbody_system%ke_spin + snapshot%pe = nbody_system%pe + snapshot%be = nbody_system%be + snapshot%te = nbody_system%te + snapshot%oblpot = nbody_system%oblpot + snapshot%L_orbit = nbody_system%L_orbit + snapshot%L_spin = nbody_system%L_spin + snapshot%L_total = nbody_system%L_total + snapshot%ke_orbit_orig = nbody_system%ke_orbit_orig + snapshot%ke_spin_orig = nbody_system%ke_spin_orig + snapshot%pe_orig = nbody_system%pe_orig + snapshot%be_orig = nbody_system%be_orig + snapshot%E_orbit_orig = nbody_system%E_orbit_orig + snapshot%GMtot_orig = nbody_system%GMtot_orig + snapshot%L_total_orig = nbody_system%L_total_orig + snapshot%L_orbit_orig = nbody_system%L_orbit_orig + snapshot%L_spin_orig = nbody_system%L_spin_orig + snapshot%L_escape = nbody_system%L_escape + snapshot%GMescape = nbody_system%GMescape + snapshot%E_collisions = nbody_system%E_collisions + snapshot%E_untracked = nbody_system%E_untracked + snapshot%ke_orbit_error = nbody_system%ke_orbit_error + snapshot%ke_spin_error = nbody_system%ke_spin_error + snapshot%pe_error = nbody_system%pe_error + snapshot%be_error = nbody_system%be_error + snapshot%E_orbit_error = nbody_system%E_orbit_error + snapshot%Ecoll_error = nbody_system%Ecoll_error + snapshot%E_untracked_error = nbody_system%E_untracked_error + snapshot%te_error = nbody_system%te_error + snapshot%L_orbit_error = nbody_system%L_orbit_error + snapshot%L_spin_error = nbody_system%L_spin_error + snapshot%L_escape_error = nbody_system%L_escape_error + snapshot%L_total_error = nbody_system%L_total_error + snapshot%Mtot_error = nbody_system%Mtot_error + snapshot%Mescape_error = nbody_system%Mescape_error + snapshot%lbeg = nbody_system%lbeg + + + ! Store a snapshot of the nbody_system for posterity + call base_util_snapshot_save(self, snapshot) + self%nt = self%iframe + self%nid = self%nid + 1 ! Central body + if (allocated(nbody_system%pl)) self%nid = self%nid + nbody_system%pl%nbody + if (allocated(nbody_system%tp)) self%nid = self%nid + nbody_system%tp%nbody return end subroutine swiftest_util_snapshot_system