diff --git a/CMakeLists.txt b/CMakeLists.txt index 9ce5af9cf..eb709cbe0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,7 +29,7 @@ ENDIF(NOT CMAKE_Fortran_COMPILER_SUPPORTS_F90) # Set some options the user may choose # Uncomment the below if you want the user to choose a parallelization library -OPTION(USE_MPI "Use the MPI library for parallelization" ON) +OPTION(USE_CAF "Use Coarray Fortran for parallelization" ON) OPTION(USE_OPENMP "Use OpenMP for parallelization" ON) diff --git a/cmake/Modules/SetParallelizationLibrary.cmake b/cmake/Modules/SetParallelizationLibrary.cmake index d809cc00b..3f3153f4d 100644 --- a/cmake/Modules/SetParallelizationLibrary.cmake +++ b/cmake/Modules/SetParallelizationLibrary.cmake @@ -7,11 +7,6 @@ # You should have received a copy of the GNU General Public License along with Swiftest. # If not, see: https://www.gnu.org/licenses. -# Turns on either OpenMP or MPI -# If both are requested, the other is disabled -# When one is turned on, the other is turned off -# If both are off, we explicitly disable them just in case - IF (USE_OPENMP) # Find OpenMP IF (NOT OpenMP_Fortran_FLAGS) @@ -20,28 +15,23 @@ IF (USE_OPENMP) MESSAGE (FATAL_ERROR "Fortran compiler does not support OpenMP") ENDIF (NOT OpenMP_Fortran_FLAGS) ENDIF (NOT OpenMP_Fortran_FLAGS) - # Turn of MPI ENDIF (USE_OPENMP) -IF (USE_MPI) +IF (USE_CAF) # Find MPI - IF (NOT MPI_Fortran_FOUND) - FIND_PACKAGE (MPI REQUIRED) + IF (NOT Coarray_Fortran_FLAGS) FIND_PACKAGE (Coarray_Fortran) IF (NOT Coarray_Fortran_FLAGS) MESSAGE (FATAL_ERROR "Fortran compiler does not support Coarrays") ENDIF (NOT Coarray_Fortran_FLAGS) - ENDIF (NOT MPI_Fortran_FOUND) -ENDIF (USE_MPI) + ENDIF (NOT Coarray_Fortran_FLAGS) +ENDIF (USE_CAF) -IF (NOT USE_OPENMP AND NOT USE_MPI) - # Turn off both OpenMP and MPI +IF (NOT USE_OPENMP AND NOT USE_CAF) + # Turn off both OpenMP and CAF SET (OMP_NUM_PROCS 0 CACHE STRING "Number of processors OpenMP may use" FORCE) UNSET (OpenMP_Fortran_FLAGS CACHE) UNSET (Coarray_Fortran_FLAGS CACHE) UNSET (GOMP_Fortran_LINK_FLAGS CACHE) - UNSET (MPI_FOUND CACHE) - UNSET (MPI_COMPILER CACHE) - UNSET (MPI_LIBRARY CACHE) -ENDIF (NOT USE_OPENMP AND NOT USE_MPI) +ENDIF (NOT USE_OPENMP AND NOT USE_CAF) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c41433ede..b8cc5be2c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -101,13 +101,11 @@ IF(USE_OPENMP) LINK_FLAGS "${OpenMP_Fortran_FLAGS}") ENDIF(USE_OPENMP) -IF(USE_MPI) +IF(USE_CAF) SET_TARGET_PROPERTIES(${SWIFTEST_DRIVER} PROPERTIES - COMPILE_FLAGS "${MPI_Fortran_COMPILE_FLAGS} ${Coarray_Fortran_FLAGS}" - LINK_FLAGS "${MPI_Fortran_LINK_FLAGS} ${Coarray_Fortran_FLAGS}") - INCLUDE_DIRECTORIES(${MPI_Fortran_INCLUDE_PATH}) - TARGET_LINK_LIBRARIES(${SWIFTEST_DRIVER} ${MPI_Fortran_LIBRARIES}) -ENDIF(USE_MPI) + COMPILE_FLAGS "${Coarray_Fortran_FLAGS}" + LINK_FLAGS "${Coarray_Fortran_FLAGS}") +ENDIF(USE_CAF) ##################################### diff --git a/src/rmvs/rmvs_module.f90 b/src/rmvs/rmvs_module.f90 index 3a950bc8d..74910d7d3 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(swiftest_nbody_system), dimension(:), allocatable :: planetocentric !! Planetocentric version of the massive body objects (one for each massive body) + class(base_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 bcb38797a..18f9bd644 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -451,64 +451,67 @@ subroutine rmvs_make_planetocentric(param, cb, pl, tp) if (allocated(encmask)) deallocate(encmask) allocate(encmask(ntp)) encmask(1:ntp) = tp%plencP(1:ntp) == i - 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) + 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 end do - end do - call tpenci%set_mu(cbenci) - end associate + call tpenci%set_mu(cbenci) + end associate + end select end select end select end select @@ -611,39 +614,42 @@ 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(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 + 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 end select end select end select diff --git a/src/swiftest/swiftest_driver.f90 b/src/swiftest/swiftest_driver.f90 index 2da264bea..336578eb2 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(swiftest_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated + class(base_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,75 +70,78 @@ 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) - - !> 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 + 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.) end if - 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) + 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, nsubsteps=istep_out) + call nbody_system%display_run_information(param, integration_timer) + call integration_timer%reset() + if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) - if (idump == dump_cadence) then - idump = 0 - call nbody_system%dump(param) end if - - call integration_timer%report(message="Integration steps:", unit=display_unit, nsubsteps=istep_out) - 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 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 associate call base_util_exit(SUCCESS) diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 9f8b7e0c1..92e1a4006 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -304,7 +304,7 @@ module swiftest !> An abstract class for a basic Swiftest nbody system - type, abstract, extends(base_nbody_system) :: swiftest_nbody_system + type, 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 @@ -373,7 +373,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,7 +390,8 @@ 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 :: dealloc => swiftest_util_dealloc_system !! Deallocates all allocatables and resets all values to defaults. Acts as a base for a finalizer + 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 :: 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 @@ -447,15 +448,6 @@ 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 @@ -1057,7 +1049,7 @@ end subroutine swiftest_util_setup_body module subroutine swiftest_util_setup_construct_system(nbody_system, param) implicit none - class(swiftest_nbody_system), allocatable, intent(inout) :: nbody_system !! Swiftest nbody_system object + class(base_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 @@ -1951,6 +1943,16 @@ 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 cc1b9e398..c68dd1ac3 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(swiftest_nbody_system), allocatable, intent(inout) :: nbody_system !! Swiftest nbody_system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(base_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 - allocate(swiftest_particle_info :: nbody_system%cb%info) - - - nbody_system%t = param%tstart - + select type(nbody_system) + class is (swiftest_nbody_system) + allocate(swiftest_particle_info :: nbody_system%cb%info) + nbody_system%t = param%tstart + end select 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(swiftest_nbody_system), allocatable :: snapshot + class(base_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,60 +3077,63 @@ 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) - 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 + 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 return end subroutine swiftest_util_snapshot_system diff --git a/src/whm/whm_module.f90 b/src/whm/whm_module.f90 index fc97eee3b..d2d10f971 100644 --- a/src/whm/whm_module.f90 +++ b/src/whm/whm_module.f90 @@ -70,7 +70,6 @@ module whm !> An abstract class for the WHM integrator nbody system type, extends(swiftest_nbody_system) :: whm_nbody_system - class(whm_tp), dimension(:),codimension[:], allocatable :: cotp !! Co-array test particle data structure contains !> Replace the abstract procedures with concrete ones procedure :: initialize => whm_util_setup_initialize_system !! Performs WHM-specific initilization steps, like calculating the Jacobi masses