diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b2725e909..c2d95aaad 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -75,6 +75,7 @@ SET(FAST_MATH_FILES ${SRC}/rmvs/rmvs_util.f90 ${SRC}/swiftest/swiftest_discard.f90 ${SRC}/swiftest/swiftest_util.f90 + ${SRC}/swiftest/swiftest_driver.f90 ${SRC}/symba/symba_discard.f90 ${SRC}/symba/symba_encounter_check.f90 ${SRC}/symba/symba_util.f90 @@ -92,7 +93,7 @@ SET(COARRAY_FILES ${SRC}/rmvs/rmvs_coarray.f90 ) -SET(DRIVER_src ${SRC}/swiftest/swiftest_driver.f90) +SET(DRIVER_src ${SRC}/main/main.f90) # Combine all source files (includ) set(SWIFTEST_src ${FAST_MATH_FILES} ${STRICT_MATH_FILES}) diff --git a/src/main/main.f90 b/src/main/main.f90 new file mode 100644 index 000000000..250c05023 --- /dev/null +++ b/src/main/main.f90 @@ -0,0 +1,28 @@ +!! Copyright 2023 - David Minton +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +program main + !! author: David A. Minton + !! + !! This is used as a wrapper for the swiftest_driver subroutine so that the driver can be used as either + !! a function call or a standalone executable program. + use swiftest + implicit none + + 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 + character(len=:), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD" + + ! Read in command line arguments + call swiftest_io_get_args(integrator, param_file_name, display_style) + + ! Execute the driver + call swiftest_driver(integrator, param_file_name, display_style) + +end program main \ No newline at end of file diff --git a/src/swiftest/swiftest_driver.f90 b/src/swiftest/swiftest_driver.f90 index 3aea4aab7..b4dd7c997 100644 --- a/src/swiftest/swiftest_driver.f90 +++ b/src/swiftest/swiftest_driver.f90 @@ -6,186 +6,194 @@ !! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. !! You should have received a copy of the GNU General Public License along with Swiftest. !! If not, see: https://www.gnu.org/licenses. - -program swiftest_driver - !! author: David A. Minton - !! - !! Driver program for the Swiftest integrators. Unlike the earlier Swift and Swifter drivers, in Swiftest all integrators - !! are run from this single program. - !! - !! Adapted from Swifter by David E. Kaufmann's Swifter driver programs swifter_[bs,helio,ra15,rmvs,symba,tu4,whm].f90 - !! Adapted from Hal Levison and Martin Duncan's Swift driver programs - use netcdf - use swiftest - implicit none - - class(swiftest_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated - type(swiftest_parameters) :: param !! Run configuration parameters - class(swiftest_storage), allocatable :: system_history !! Stores the system history between output dumps - 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 - character(len=:), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD" - type(walltimer) :: integration_timer !! Object used for computing elapsed wall time - - call swiftest_io_get_args(integrator, param_file_name, display_style) - - !> Read in the user-defined parameters file and the initial conditions of the nbody_system - param%integrator = trim(adjustl(integrator)) - param%display_style = trim(adjustl(display_style)) - call param%read_in(param_file_name) +submodule(swiftest) s_swiftest_driver +contains + + subroutine swiftest_driver(integrator, param_file_name, display_style) + + !! author: David A. Minton + !! + !! Driver program for the Swiftest integrators. Unlike the earlier Swift and Swifter drivers, in Swiftest all integrators + !! are run from this single program. + !! + !! Adapted from Swifter by David E. Kaufmann's Swifter driver programs swifter_[bs,helio,ra15,rmvs,symba,tu4,whm].f90 + !! Adapted from Hal Levison and Martin Duncan's Swift driver programs + use netcdf + use swiftest + implicit none + + ! Arguments + character(len=:), intent(in), allocatable :: integrator !! Symbolic code of the requested integrator + character(len=:), intent(in), allocatable :: param_file_name !! Name of the input parameters file + character(len=:), intent(in), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD" + + ! Internals + class(swiftest_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated + type(swiftest_parameters) :: param !! Run configuration parameters + class(swiftest_storage), allocatable :: system_history !! Stores the system history between output dumps + type(walltimer) :: integration_timer !! Object used for computing elapsed wall time + + !> Read in the user-defined parameters file and the initial conditions of the nbody_system + param%integrator = trim(adjustl(integrator)) + param%display_style = trim(adjustl(display_style)) + call param%read_in(param_file_name) #ifdef COARRAY - if (.not.param%lcoarray .and. (this_image() /= 1)) stop ! Single image mode + if (.not.param%lcoarray .and. (this_image() /= 1)) stop ! Single image mode #endif - associate(t0 => param%t0, & - tstart => param%tstart, & - dt => param%dt, & - tstop => param%tstop, & - iloop => param%iloop, & - istart => param%istart, & - iout => param%iout, & - idump => param%idump, & - nout => param%nout, & - istep => param%istep, & - nloops => param%nloops, & - istep_out => param%istep_out, & - fstep_out => param%fstep_out, & - ltstretch => param%ltstretch, & - dump_cadence => param%dump_cadence, & - display_unit => param%display_unit) - - ! Set up loop and output cadence variables - nloops = ceiling((tstop - t0) / dt, kind=I8B) - istart = ceiling((tstart - t0) / dt + 1.0_DP, kind=I8B) - iloop = istart - 1 - iout = 0 - idump = 0 - if (ltstretch) then - nout = floor(log(1.0_DP + iloop * (fstep_out - 1.0_DP) / istep_out) / log(fstep_out)) - istep = floor(istep_out * fstep_out**nout, kind=I4B) - else - nout = 1 - istep = istep_out - end if - - ! Set up nbody_system storage for intermittent file dumps - if (dump_cadence == 0) dump_cadence = int(ceiling(nloops / (1.0_DP * istep_out), kind=I8B), kind=I4B) - - ! 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 + associate(t0 => param%t0, & + tstart => param%tstart, & + dt => param%dt, & + tstop => param%tstop, & + iloop => param%iloop, & + istart => param%istart, & + iout => param%iout, & + idump => param%idump, & + nout => param%nout, & + istep => param%istep, & + nloops => param%nloops, & + istep_out => param%istep_out, & + fstep_out => param%fstep_out, & + ltstretch => param%ltstretch, & + dump_cadence => param%dump_cadence, & + display_unit => param%display_unit) + + ! Set up loop and output cadence variables + nloops = ceiling((tstop - t0) / dt, kind=I8B) + istart = ceiling((tstart - t0) / dt + 1.0_DP, kind=I8B) + iloop = istart - 1 + iout = 0 + idump = 0 + if (ltstretch) then + nout = floor(log(1.0_DP + iloop * (fstep_out - 1.0_DP) / istep_out) / log(fstep_out)) + istep = floor(istep_out * fstep_out**nout, kind=I4B) + else + nout = 1 + istep = istep_out + end if + + ! Set up nbody_system storage for intermittent file dumps + if (dump_cadence == 0) dump_cadence = int(ceiling(nloops / (1.0_DP * istep_out), kind=I8B), kind=I4B) + + ! 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 #ifdef COARRAY - if (this_image() == 1 .or. param%log_output) then + if (this_image() == 1 .or. param%log_output) then #endif - !$ 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 + !$ 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 #ifdef COARRAY - if (param%lcoarray) then - write(param%display_unit,*) ' Coarray parameters:' - write(param%display_unit,*) ' -------------------' - write(param%display_unit,*) ' Number of images = ', num_images() - if (param%log_output .and. this_image() == 1) write(*,'(a,i3)') ' Coarray: Number of images = ',num_images() - else - write(param%display_unit,*) ' Coarrays disabled.' - if (param%log_output) write(*,*) ' Coarrays disabled.' + if (param%lcoarray) then + write(param%display_unit,*) ' Coarray parameters:' + write(param%display_unit,*) ' -------------------' + write(param%display_unit,*) ' Number of images = ', num_images() + if (param%log_output .and. this_image() == 1) write(*,'(a,i3)') ' Coarray: Number of images = ',num_images() + else + write(param%display_unit,*) ' Coarrays disabled.' + if (param%log_output) write(*,*) ' Coarrays disabled.' + end if end if - end if #endif - if (param%log_output) flush(param%display_unit) + if (param%log_output) flush(param%display_unit) #ifdef COARRAY - ! The following line lets us read in the input files one image at a time. Letting each image read the input in is faster than broadcasting all of the data - if (param%lcoarray .and. (this_image() /= 1)) sync images(this_image() - 1) + ! The following line lets us read in the input files one image at a time. Letting each image read the input in is faster than broadcasting all of the data + if (param%lcoarray .and. (this_image() /= 1)) sync images(this_image() - 1) #endif - call nbody_system%initialize(system_history, param) + call nbody_system%initialize(system_history, param) #ifdef COARRAY - if (param%lcoarray .and. (this_image() < num_images())) sync images(this_image() + 1) + if (param%lcoarray .and. (this_image() < num_images())) sync images(this_image() + 1) - ! Distribute test particles to the various images - if (param%lcoarray) call nbody_system%coarray_distribute(param) + ! Distribute test particles to the various images + if (param%lcoarray) call nbody_system%coarray_distribute(param) #endif - ! 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 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(system_history%nc, param) - else - call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum + if (param%lenergy) then + if (param%lrestart) then + call nbody_system%get_t0_values(system_history%nc, 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, system_history) - - 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, system_history) - if (idump == dump_cadence) then - idump = 0 - call nbody_system%dump(param, system_history) + 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, system_history) #ifdef COARRAY - if (param%lcoarray) call nbody_system%coarray_balance(param) + if (param%lcoarray) call nbody_system%coarray_balance(param) #endif - end if + end if #ifdef COARRAY - if (this_image() == 1 .or. param%log_output) then + if (this_image() == 1 .or. param%log_output) then #endif - call integration_timer%report(message="Integration steps:", unit=display_unit) + call integration_timer%report(message="Integration steps:", unit=display_unit) #ifdef COARRAY - end if !(this_image() == 1) + end if !(this_image() == 1) #endif - call nbody_system%display_run_information(param, integration_timer) - call integration_timer%reset() + call nbody_system%display_run_information(param, integration_timer) + call integration_timer%reset() #ifdef COARRAY - if (this_image() == 1 .or. param%log_output) then + if (this_image() == 1 .or. param%log_output) then #endif - if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) + if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) #ifdef COARRAY - end if ! (this_image() == 1) + end if ! (this_image() == 1) #endif + end if end if - end if - end do - ! Dump any remaining history if it exists - call nbody_system%dump(param, system_history) - 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, system_history) + call nbody_system%display_run_information(param, integration_timer, phase="last") + end associate #ifdef COARRAY - if (this_image() == 1) then + if (this_image() == 1) then #endif - call base_util_exit(SUCCESS,unit=param%display_unit) + call base_util_exit(SUCCESS,unit=param%display_unit) #ifdef COARRAY - end if ! (this_image() == 1) + end if ! (this_image() == 1) #endif -end program swiftest_driver + + return + end subroutine swiftest_driver + +end submodule s_swiftest_driver diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index b7bbd109c..ffb60a593 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -536,6 +536,13 @@ pure elemental module subroutine swiftest_drift_one(mu, rx, ry, rz, vx, vy, vz, integer(I4B), intent(out) :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR) end subroutine swiftest_drift_one + module subroutine swiftest_driver(integrator, param_file_name, display_style) + implicit none + character(len=:), intent(in), allocatable :: integrator !! Symbolic code of the requested integrator + character(len=:), intent(in), allocatable :: param_file_name !! Name of the input parameters file + character(len=:), intent(in), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD" + end subroutine swiftest_driver + pure module subroutine swiftest_gr_kick_getaccb_ns_body(self, nbody_system, param) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest generic body object