From 8f4e883353e6d6fd5e3a00c39ba12d8dcf0af5e8 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 2 Jul 2021 17:19:47 -0400 Subject: [PATCH] Restructured the swiftest and whm classes to accept system variabiles in place of individual cb, pl, and tp for most cases --- src/main/swiftest_driver.f90 | 8 +- src/modules/rmvs_classes.f90 | 21 +-- src/modules/swiftest_classes.f90 | 218 ++++++++++++++----------------- src/modules/whm_classes.f90 | 112 ++++++++-------- src/rmvs/rmvs_step.f90 | 46 +++---- src/whm/whm_drift.f90 | 38 +++--- src/whm/whm_getacch.f90 | 46 +++---- src/whm/whm_gr.f90 | 24 ++-- src/whm/whm_step.f90 | 51 ++++---- 9 files changed, 263 insertions(+), 301 deletions(-) diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index e88118d12..3c9f3c8ac 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -48,7 +48,7 @@ program swiftest_driver iloop = 0 iout = istep_out idump = istep_dump - if (istep_out > 0) call nbody_system%write_frame(iu, param, t, dt) + if (istep_out > 0) call nbody_system%write_frame(iu, param) !> Define the maximum number of threads nthreads = 1 ! In the *serial* case !$ nthreads = omp_get_max_threads() ! In the *parallel* case @@ -60,7 +60,7 @@ program swiftest_driver ntp = nbody_system%tp%nbody npl = nbody_system%pl%nbody !> Step the system forward in time - call nbody_system%step(param) + call nbody_system%step(param, t, dt) t = t0 + iloop * dt if (t > tstop) exit @@ -72,7 +72,7 @@ program swiftest_driver if (istep_out > 0) then iout = iout - 1 if (iout == 0) then - call nbody_system%write_frame(iu, param, t, dt) + call nbody_system%write_frame(iu, param) iout = istep_out end if end if @@ -81,7 +81,7 @@ program swiftest_driver if (istep_dump > 0) then idump = idump - 1 if (idump == 0) then - call nbody_system%dump(param, t, dt, statusfmt) + call nbody_system%dump(param, statusfmt) idump = istep_dump end if end if diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 4dea2e59f..b67c354b0 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -117,16 +117,15 @@ module subroutine rmvs_discard_pl_tp(self, pl, t, dt) real(DP), intent(in) :: dt !! Stepsize end subroutine rmvs_discard_pl_tp - module subroutine rmvs_getacch_tp(self, cb, pl, param, t, xh) + module subroutine rmvs_getacch_tp(self, system, param, t, xh) use swiftest_classes, only : swiftest_cb, swiftest_parameters - use whm_classes, only : whm_pl + use whm_classes, only : whm_nbody_system implicit none - class(rmvs_tp), intent(inout) :: self !! RMVS test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree - class(whm_pl), intent(inout) :: pl !! WHM massive body particle data structure. + class(rmvs_tp), intent(inout) :: self !! RMVS test particle data structure + class(whm_nbody_system), intent(inout) :: system !! Swiftest central body particle data structuree class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameter - real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets + real(DP), intent(in) :: t !! Current time + real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets end subroutine rmvs_getacch_tp module subroutine rmvs_setup_system(self, param) @@ -148,11 +147,13 @@ module subroutine rmvs_setup_set_beg_end(self, xbeg, xend, vbeg) real(DP), dimension(:,:), intent(in), optional :: xbeg, xend, vbeg end subroutine rmvs_setup_set_beg_end - module subroutine rmvs_step_system(self, param) + module subroutine rmvs_step_system(self, param, t, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters + class(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters + real(DP), intent(in) :: t !! Simulation time + real(DP), intent(in) :: dt !! Current stepsize end subroutine rmvs_step_system module subroutine rmvs_setup_tp(self,n) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 5dd39f74d..62ea97d6f 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -142,10 +142,8 @@ module swiftest_classes contains private procedure(abstract_set_mu), public, deferred :: set_mu - procedure(abstract_gr_getacch), public, deferred :: gr_getacch procedure(abstract_step_body), public, deferred :: step ! These are concrete because the implementation is the same for all types of particles - procedure, public :: gr_getaccb => gr_getaccb_ns_body !! Add relativistic correction acceleration for non-symplectic integrators procedure, public :: initialize => io_read_body_in !! Read in body initial conditions from a file procedure, public :: read_frame => io_read_frame_body !! I/O routine for writing out a single frame of time-series data for the central body procedure, public :: write_frame => io_write_frame_body !! I/O routine for writing out a single frame of time-series data for the central body @@ -161,7 +159,6 @@ module swiftest_classes procedure, public :: fill => util_fill_body !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) procedure, public :: spill => util_spill_body !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) procedure, public :: reverse_status => util_reverse_status !! Reverses the active/inactive status of all particles in a structure - procedure, public :: step end type swiftest_body !******************************************************************************************************************************** @@ -270,13 +267,6 @@ subroutine abstract_copy(self, src, mask) logical, dimension(:), intent(in) :: mask end subroutine abstract_copy - subroutine abstract_gr_getacch(self, cb, param) - import swiftest_body, swiftest_cb, swiftest_parameters - class(swiftest_body), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameter - end subroutine abstract_gr_getacch - subroutine abstract_initialize(self, param) import swiftest_base, swiftest_parameters class(swiftest_base), intent(inout) :: self !! Swiftest base object @@ -299,48 +289,50 @@ subroutine abstract_set_mu(self, cb) class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object end subroutine abstract_set_mu - subroutine abstract_step_body(self, system, param, dt) - import swiftest_nbody_system, swiftest_parameters + subroutine abstract_step_body(self, system, param, t, dt) + import DP, swiftest_body, swiftest_nbody_system, swiftest_parameters implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest particle object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters + class(swiftest_body), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters + 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) - import swiftest_nbody_system, swiftest_parameters + 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 system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters + real(DP), intent(in) :: t !! Simulation time + real(DP), intent(in) :: dt !! Current stepsize end subroutine abstract_step_system - subroutine abstract_write_frame(self, iu, param, t, dt) + subroutine abstract_write_frame(self, iu, param) import DP, I4B, swiftest_base, swiftest_parameters - class(swiftest_base), intent(in) :: self !! Swiftest base object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_base), intent(in) :: self !! Swiftest base object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Step size end subroutine abstract_write_frame end interface interface module subroutine discard_peri_tp(self, cb, pl, param, t, msys) implicit none - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object class(swiftest_parameters), intent(in) :: param !! parameters - real(DP), intent(in) :: t !! Current simulation tim - real(DP), intent(in) :: msys !! Total system mass + real(DP), intent(in) :: t !! Current simulation tim + real(DP), intent(in) :: msys !! Total system mass end subroutine discard_peri_tp module subroutine discard_pl_close(dx, dv, dt, r2crit, iflag, r2min) implicit none - real(DP), dimension(:), intent(in) :: dx, dv - real(DP), intent(in) :: dt, r2crit - integer(I4B), intent(out) :: iflag - real(DP), intent(out) :: r2min + real(DP), dimension(:), intent(in) :: dx, dv + real(DP), intent(in) :: dt, r2crit + integer(I4B), intent(out) :: iflag + real(DP), intent(out) :: r2min end subroutine discard_pl_close module subroutine discard_pl_tp(self, pl, t, dt) @@ -351,19 +343,18 @@ module subroutine discard_pl_tp(self, pl, t, dt) real(DP), intent(in) :: dt !! Stepsize end subroutine discard_pl_tp - module subroutine discard_sun_tp(self, cb, param, t, msys) + module subroutine discard_sun_tp(self, cb, param, msys) implicit none - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object class(swiftest_parameters), intent(in) :: param !! parameters - real(DP), intent(in) :: t !! Current simulation tim - real(DP), intent(in) :: msys !! Total system mass + real(DP), intent(in) :: msys !! Total system mass end subroutine discard_sun_tp module subroutine discard_system(self, param) implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters end subroutine discard_system module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) @@ -376,33 +367,24 @@ end subroutine drift_one module subroutine eucl_dist_index_plpl(self) implicit none - class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object end subroutine module subroutine eucl_dist_index_pltp(self, pl) implicit none - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object end subroutine module subroutine eucl_irij3_plpl(self) implicit none - class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object end subroutine eucl_irij3_plpl - module subroutine gr_getaccb_ns_body(self, cb, param, agr, agr0) - implicit none - class(swiftest_body), intent(inout) :: self - class(swiftest_cb), intent(inout) :: cb - class(swiftest_parameters), intent(in) :: param - real(DP), dimension(:, :), intent(inout) :: agr - real(DP), dimension(NDIM), intent(out) :: agr0 - end subroutine gr_getaccb_ns_body - module subroutine kick_vb_body(self, dt) implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest generic body object - real(DP), intent(in) :: dt !! Stepsize + class(swiftest_body), intent(inout) :: self !! Swiftest generic body object + real(DP), intent(in) :: dt !! Stepsize end subroutine kick_vb_body module subroutine kick_vh_body(self, dt) @@ -412,47 +394,44 @@ module subroutine kick_vh_body(self, dt) end subroutine kick_vh_body module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) + implicit none 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. - !! If you do not include a char-literal-constant, the iotype argument contains only DT. - integer, 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 + 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. + integer(I4B), 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 end subroutine io_param_reader module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) - class(swiftest_parameters),intent(in) :: 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. - !! If you do not include a char-literal-constant, the iotype argument contains only DT. - integer, intent(in) :: v_list(:) !! Not used in this procedure - integer, intent(out) :: iostat !! IO status code - character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 + implicit none + class(swiftest_parameters), intent(in) :: self !! Collection of parameters + 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. + integer(I4B), intent(in) :: v_list(:) !! Not used in this procedure + integer(I4B), intent(out) :: iostat !! IO status code + character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 end subroutine io_param_writer - module subroutine io_dump_param(self, param_file_name, t, dt) - class(swiftest_parameters),intent(in) :: self !! Output collection of parameters - character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) - real(DP),intent(in) :: t !! Current simulation time - real(DP),intent(in) :: dt !! Step size + module subroutine io_dump_param(self, param_file_name) + implicit none + class(swiftest_parameters),intent(in) :: self !! Output collection of parameters + character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) end subroutine io_dump_param - module subroutine io_dump_swiftest(self, param, t, dt, msg) + module subroutine io_dump_swiftest(self, param, msg) implicit none class(swiftest_base), intent(inout) :: self !! Swiftest base object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Stepsize + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters character(*), optional, intent(in) :: msg !! Message to display with dump operation end subroutine io_dump_swiftest - module subroutine io_dump_system(self, param, t, dt, msg) + module subroutine io_dump_system(self, param, msg) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Stepsize + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters character(*), optional, intent(in) :: msg !! Message to display with dump operation end subroutine io_dump_system @@ -473,59 +452,60 @@ end function io_get_token module subroutine io_read_body_in(self, param) implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest particle object + class(swiftest_body), intent(inout) :: self !! Swiftest particle object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters end subroutine io_read_body_in module subroutine io_read_cb_in(self, param) implicit none - class(swiftest_cb), intent(inout) :: self + class(swiftest_cb), intent(inout) :: self class(swiftest_parameters), intent(inout) :: param end subroutine io_read_cb_in module subroutine io_read_param_in(self, param_file_name) - class(swiftest_parameters),intent(out) :: self !! Current run configuration parameters of parameters - character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) + implicit none + class(swiftest_parameters), intent(out) :: self !! Current run configuration parameters of parameters + character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) end subroutine io_read_param_in module function io_read_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & xh1, xh2, vh1, vh2, encounter_file, out_type) result(ierr) implicit none - integer(I4B) :: ierr - integer(I4B), intent(out) :: name1, name2 - real(DP), intent(out) :: t, mass1, mass2, radius1, radius2 - real(DP), dimension(NDIM), intent(out) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: encounter_file,out_type + integer(I4B) :: ierr + integer(I4B), intent(out) :: name1, name2 + real(DP), intent(out) :: t, mass1, mass2, radius1, radius2 + real(DP), dimension(NDIM), intent(out) :: xh1, xh2, vh1, vh2 + character(*), intent(in) :: encounter_file,out_type end function io_read_encounter module subroutine io_read_frame_body(self, iu, param, form, t, ierr) implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest particle object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_body), intent(inout) :: self !! Swiftest particle object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - real(DP), intent(out) :: t !! Simulation time - integer(I4B), intent(out) :: ierr !! Error code + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + real(DP), intent(out) :: t !! Simulation time + integer(I4B), intent(out) :: ierr !! Error code end subroutine io_read_frame_body module subroutine io_read_frame_cb(self, iu, param, form, t, ierr) implicit none - class(swiftest_cb), intent(inout) :: self !! Swiftest central body object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_cb), intent(inout) :: self !! Swiftest central body object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - real(DP), intent(out) :: t !! Simulation time - integer(I4B), intent(out) :: ierr !! Error code + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + real(DP), intent(out) :: t !! Simulation time + integer(I4B), intent(out) :: ierr !! Error code end subroutine io_read_frame_cb module subroutine io_read_frame_system(self, iu, param, form, t, ierr) implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - real(DP), intent(out) :: t !! Current simulation time - integer(I4B), intent(out) :: ierr !! Error code + class(swiftest_nbody_system),intent(inout) :: self !! Swiftest system object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + real(DP), intent(out) :: t !! Current simulation time + integer(I4B), intent(out) :: ierr !! Error code end subroutine io_read_frame_system module function io_read_hdr(iu, t, npl, ntp, out_form, out_type) result(ierr) @@ -560,31 +540,25 @@ module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, rad character(*), intent(in) :: encounter_file, out_type end subroutine io_write_encounter - module subroutine io_write_frame_body(self, iu, param, t, dt) + module subroutine io_write_frame_body(self, iu, param) implicit none class(swiftest_body), intent(in) :: self !! Swiftest particle object integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Step size end subroutine io_write_frame_body - module subroutine io_write_frame_cb(self, iu, param, t, dt) + module subroutine io_write_frame_cb(self, iu, param) implicit none class(swiftest_cb), intent(in) :: self !! Swiftest central body object integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Step size end subroutine io_write_frame_cb - module subroutine io_write_frame_system(self, iu, param, t, dt) + module subroutine io_write_frame_system(self, iu, param) implicit none class(swiftest_nbody_system), intent(in) :: self !! Swiftest system object integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Step size + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters end subroutine io_write_frame_system module subroutine io_write_hdr(iu, t, npl, ntp, out_form, out_type) @@ -688,10 +662,10 @@ end subroutine setup_tp module subroutine user_getacch_body(self, cb, param, t) implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree + class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: t !! Current time end subroutine user_getacch_body module subroutine util_coord_b2h_pl(self, cb) @@ -720,7 +694,7 @@ end subroutine util_coord_h2b_tp module subroutine util_copy_body(self, src, mask) implicit none - class(swiftest_body), intent(inout) :: self + class(swiftest_body), intent(inout) :: self class(swiftest_base), intent(in) :: src logical, dimension(:), intent(in) :: mask end subroutine util_copy_body @@ -756,8 +730,8 @@ end subroutine util_copy_system module subroutine util_fill_body(self, inserts, lfill_list) implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest generic body object - class(swiftest_body), intent(inout) :: inserts !! Insertted object + class(swiftest_body), intent(inout) :: self !! Swiftest generic body object + class(swiftest_body), intent(inout) :: inserts !! Insertted object logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps end subroutine util_fill_body diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index af6d62106..282ced72a 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -107,130 +107,127 @@ module subroutine whm_setup_set_ir3j(self) class(whm_pl), intent(inout) :: self !! WHM massive body object end subroutine whm_setup_set_ir3j - module subroutine whm_step_pl(self, cb, param, t, dt) - use swiftest_classes, only : swiftest_cb, swiftest_parameters + module subroutine whm_step_pl(self, system, param, t, dt) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time - real(DP), intent(in) :: dt !! Stepsize + class(whm_pl), intent(inout) :: self !! WHM massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters + real(DP), intent(in) :: t !! Simulation time + real(DP), intent(in) :: dt !! Current stepsize end subroutine whm_step_pl !> Get heliocentric accelration of massive bodies - module subroutine whm_getacch_pl(self, cb, param, t) + module subroutine whm_getacch_pl(self, system, param, t) use swiftest_classes, only : swiftest_cb, swiftest_parameters implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + class(whm_nbody_system), intent(inout) :: system !! WHM nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of + real(DP), intent(in) :: t !! Current simulation time end subroutine whm_getacch_pl - module subroutine whm_drift_pl(self, cb, param, dt) + module subroutine whm_drift_pl(self, system, param, dt) use swiftest_classes, only : swiftest_cb, swiftest_parameters implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structur - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: dt !! Stepsize + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + class(whm_nbody_system), intent(inout) :: system !! WHM nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of + real(DP), intent(in) :: dt !! Stepsize end subroutine whm_drift_pl module subroutine whm_getacch_int_pl(self, cb) use swiftest_classes, only : swiftest_cb implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure end subroutine whm_getacch_int_pl module subroutine whm_coord_h2j_pl(self, cb) use swiftest_classes, only : swiftest_cb implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree end subroutine whm_coord_h2j_pl module subroutine whm_coord_j2h_pl(self, cb) use swiftest_classes, only : swiftest_cb implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree end subroutine whm_coord_j2h_pl module subroutine whm_coord_vh2vj_pl(self, cb) use swiftest_classes, only : swiftest_cb implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree end subroutine whm_coord_vh2vj_pl - module subroutine whm_gr_getacch_pl(self, cb, param) + module subroutine whm_gr_getacch_pl(self, param) use swiftest_classes, only : swiftest_cb, swiftest_parameters implicit none class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of end subroutine whm_gr_getacch_pl module pure subroutine whm_gr_p4_pl(self, param, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(whm_pl), intent(inout) :: self !! Swiftest particle object + class(whm_pl), intent(inout) :: self !! Swiftest particle object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters - real(DP), intent(in) :: dt !! Step size + real(DP), intent(in) :: dt !! Step size end subroutine whm_gr_p4_pl module pure subroutine whm_gr_vh2pv_pl(self, param) use swiftest_classes, only : swiftest_parameters implicit none - class(whm_pl), intent(inout) :: self !! Swiftest particle object + class(whm_pl), intent(inout) :: self !! Swiftest particle object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters end subroutine whm_gr_vh2pv_pl module pure subroutine whm_gr_pv2vh_pl(self, param) use swiftest_classes, only : swiftest_parameters implicit none - class(whm_pl), intent(inout) :: self !! Swiftest particle object + class(whm_pl), intent(inout) :: self !! Swiftest particle object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters end subroutine whm_gr_pv2vh_pl !> WHM test particle constructor module subroutine whm_setup_tp(self,n) implicit none - class(whm_tp), intent(inout) :: self !! WHM test particle data structure - integer, intent(in) :: n !! Number of test particles to allocate + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + integer, intent(in) :: n !! Number of test particles to allocate end subroutine whm_setup_tp - module subroutine whm_drift_tp(self, cb, param, dt) + module subroutine whm_drift_tp(self, system, param, dt) use swiftest_classes, only : swiftest_cb, swiftest_parameters implicit none - class(whm_tp), intent(inout) :: self !! WHM test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + class(whm_nbody_system), intent(inout) :: system !! WHM nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: dt !! Stepsize + real(DP), intent(in) :: dt !! Stepsize end subroutine whm_drift_tp !> Get heliocentric accelration of the test particle - module subroutine whm_getacch_tp(self, cb, pl, param, t, xh) + module subroutine whm_getacch_tp(self, system, param, t, xh) use swiftest_classes, only : swiftest_cb, swiftest_parameters implicit none class(whm_tp), intent(inout) :: self !! WHM test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Generic Swiftest central body particle data structuree - class(whm_pl), intent(inout) :: pl !! WHM massive body particle data structure. - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of + class(whm_nbody_system), intent(inout) :: system !! WHM nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of real(DP), intent(in) :: t !! Current time real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets end subroutine whm_getacch_tp - module subroutine whm_step_tp(self, cb, pl, param, t, dt) - use swiftest_classes, only : swiftest_cb, swiftest_parameters + module subroutine whm_step_tp(self, system, param, t, dt) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none - class(whm_tp), intent(inout) :: self !! WHM test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure - class(whm_pl), intent(inout) :: pl !! WHM massive body data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time - real(DP), intent(in) :: dt !! Stepsize + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest central body particle data structure + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of + real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Stepsize end subroutine whm_step_tp module subroutine whm_setup_set_beg_end(self, xbeg, xend, vbeg) @@ -240,26 +237,25 @@ module subroutine whm_setup_set_beg_end(self, xbeg, xend, vbeg) real(DP), dimension(:,:), intent(in), optional :: vbeg ! vbeg is an unused variable to keep this method forward compatible with RMVS end subroutine whm_setup_set_beg_end - module subroutine whm_gr_getacch_tp(self, cb, param) + module subroutine whm_gr_getacch_tp(self, param) use swiftest_classes, only : swiftest_cb, swiftest_parameters implicit none class(whm_tp), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of end subroutine whm_gr_getacch_tp module pure subroutine whm_gr_p4_tp(self, param, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(whm_tp), intent(inout) :: self !! Swiftest particle object + class(whm_tp), intent(inout) :: self !! Swiftest particle object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters - real(DP), intent(in) :: dt !! Step size + real(DP), intent(in) :: dt !! Step size end subroutine whm_gr_p4_tp module pure subroutine whm_gr_vh2pv_tp(self, param) use swiftest_classes, only : swiftest_parameters implicit none - class(whm_tp), intent(inout) :: self !! Swiftest particle object + class(whm_tp), intent(inout) :: self !! Swiftest particle object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters end subroutine whm_gr_vh2pv_tp @@ -294,11 +290,13 @@ module subroutine whm_fill_pl(self, inserts, lfill_list) end subroutine whm_fill_pl !> Steps the Swiftest nbody system forward in time one stepsize - module subroutine whm_step_system(self, param) + module subroutine whm_step_system(self, param, t, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(whm_nbody_system), intent(inout) :: self !! WHM system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters + class(whm_nbody_system), intent(inout) :: self !! WHM system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters + real(DP), intent(in) :: t !! Simulation time + real(DP), intent(in) :: dt !! Current stepsize end subroutine whm_step_system end interface diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 0fad19590..82c6d28c4 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -1,7 +1,7 @@ submodule(rmvs_classes) s_rmvs_step use swiftest contains - module subroutine rmvs_step_system(self, param) + module subroutine rmvs_step_system(self, param, dt) !! author: David A. Minton !! !! Step massive bodies and and active test particles ahead in heliocentric coordinates @@ -10,23 +10,24 @@ module subroutine rmvs_step_system(self, param) !! Adapted from David E. Kaufmann's Swifter routine rmvs_step.f90 implicit none ! Arguments - class(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters + class(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters + integer(I4B), intent(in) :: dt !! Current stepsize ! Internals logical :: lencounter, lfirstpl, lfirsttp real(DP) :: rts real(DP), dimension(:,:), allocatable :: xbeg, xend, vbeg integer(I4B) :: i - + + select type(system => self) + class is (rmvs_nbody_system) select type(cb => self%cb) class is (rmvs_cb) select type(pl => self%pl) class is (rmvs_pl) select type(tp => self%tp) class is (rmvs_tp) - associate(ntp => tp%nbody, npl => pl%nbody, t => param%t, dt => param%dt, & - xhpl => pl%xh, vhpl => pl%vh, xjpl => pl%xj, vjpl => pl%vj, & - xhtp => tp%xh, vhtp => tp%vh) + associate(ntp => tp%nbody, npl => pl%nbody, t => param%t) allocate(xbeg, source=pl%xh) allocate(vbeg, source=pl%vh) call pl%set_rhill(cb) @@ -39,16 +40,16 @@ module subroutine rmvs_step_system(self, param) lfirsttp = tp%lfirst pl%outer(0)%x(:,:) = xbeg(:,:) pl%outer(0)%v(:,:) = vbeg(:,:) - call pl%step(cb, param, t, dt) + call pl%step(system, param, dt) pl%outer(NTENC)%x(:,:) = pl%xh(:,:) pl%outer(NTENC)%v(:,:) = pl%vh(:,:) call tp%set_beg_end(xend = pl%xh) - call rmvs_interp_out(pl,cb, dt, param) - call rmvs_step_out(pl, cb, tp, dt, param) + call rmvs_interp_out(system, param, dt) + call rmvs_step_out(system, param, dt) call tp%reverse_status() call tp%set_beg_end(xbeg = xbeg, xend = xend) tp%lfirst = .true. - call tp%step(cb, pl, param, t, dt) + call tp%step(system, param, dt) where (tp%status(:) == INACTIVE) tp%status(:) = ACTIVE pl%lfirst = lfirstpl tp%lfirst = lfirsttp @@ -59,11 +60,12 @@ module subroutine rmvs_step_system(self, param) end select end select end select + end select return end subroutine rmvs_step_system - subroutine rmvs_step_out(pl, cb, tp, dt, param) + subroutine rmvs_step_out(system, param, dt) !! author: David A. Minton !! !! Step ACTIVE test particles ahead in the outer encounter region, setting up and calling the inner region @@ -73,17 +75,15 @@ subroutine rmvs_step_out(pl, cb, tp, dt, param) !! Adapted from David E. Kaufmann's Swifter routines rmvs_step_out.f90 and rmvs_step_out2.f90 implicit none ! Arguments - class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - real(DP), intent(in) :: dt !! Step size - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameters + class(rmvs_nbody_system), intent(inout) :: system !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters + integer(I4B), intent(in) :: dt !! Current stepsize ! Internals - integer(I4B) :: outer_index, j, k - real(DP) :: dto, outer_time, rts - logical :: lencounter, lfirsttp + integer(I4B) :: outer_index, j, k + real(DP) :: dto, outer_time, rts + logical :: lencounter, lfirsttp - associate(npl => pl%nbody, ntp => tp%nbody, t => param%t) + associate(cb => system%cb, pl => system%pl, npl => system%pl%nbody, tp => system%tp, ntp => system%tp%nbody, t => param%t) dto = dt / NTENC where(tp%plencP(:) == 0) tp%status(:) = INACTIVE @@ -99,9 +99,9 @@ subroutine rmvs_step_out(pl, cb, tp, dt, param) lencounter = tp%encounter_check(cb, pl, dt, rts) if (lencounter) then ! Interpolate planets in inner encounter region - call rmvs_interp_in(pl, cb, dto, outer_index, param) + call rmvs_interp_in(system, param, dto, outer_index) ! Step through the inner region - call rmvs_step_in(pl, cb, tp, param, outer_time, dto) + call rmvs_step_in(system, param, outer_time, dto) lfirsttp = tp%lfirst tp%lfirst = .true. call tp%step(cb, pl, param, outer_time, dto) diff --git a/src/whm/whm_drift.f90 b/src/whm/whm_drift.f90 index 19213bc23..c18648364 100644 --- a/src/whm/whm_drift.f90 +++ b/src/whm/whm_drift.f90 @@ -1,7 +1,7 @@ submodule(whm_classes) whm_drift use swiftest contains - module subroutine whm_drift_pl(self, cb, param, dt) + module subroutine whm_drift_pl(self, system, param, dt) !! author: David A. Minton !! !! Loop through planets and call Danby drift routine @@ -10,15 +10,15 @@ module subroutine whm_drift_pl(self, cb, param, dt) !! Adapted from David E. Kaufmann's Swifter routine whm_drift.f90 implicit none ! Arguments - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structur - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: dt !! Stepsize + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + class(whm_nbody_system), intent(inout) :: system !! WHM nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of + real(DP), intent(in) :: dt !! Stepsize ! Internals - integer(I4B) :: i - real(DP) :: energy, vmag2, rmag !! Variables used in GR calculation - integer(I4B), dimension(:), allocatable :: iflag - real(DP), dimension(:), allocatable :: dtp + integer(I4B) :: i + real(DP) :: energy, vmag2, rmag !! Variables used in GR calculation + integer(I4B), dimension(:), allocatable :: iflag + real(DP), dimension(:), allocatable :: dtp associate(npl => self%nbody, & xj => self%xj, & @@ -63,7 +63,7 @@ module subroutine whm_drift_pl(self, cb, param, dt) end subroutine whm_drift_pl - module subroutine whm_drift_tp(self, cb, param, dt) + module subroutine whm_drift_tp(self, system, param, dt) !! author: David A. Minton !! !! Loop through test particles and call Danby drift routine @@ -73,22 +73,22 @@ module subroutine whm_drift_tp(self, cb, param, dt) !! Adapted from David E. Kaufmann's Swifter routine whm_drift_tp.f90 implicit none ! Arguments - class(whm_tp), intent(inout) :: self !! WHM test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + class(whm_nbody_system), intent(inout) :: system !! WHM nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: dt !! Stepsize + real(DP), intent(in) :: dt !! Stepsize ! Internals - integer(I4B) :: i - real(DP) :: energy, vmag2, rmag !! Variables used in GR calculation - integer(I4B), dimension(:), allocatable :: iflag - real(DP), dimension(:), allocatable :: dtp + integer(I4B) :: i + real(DP) :: energy, vmag2, rmag !! Variables used in GR calculation + integer(I4B), dimension(:), allocatable :: iflag + real(DP), dimension(:), allocatable :: dtp associate(ntp => self%nbody, & xh => self%xh, & vh => self%vh, & status => self%status,& - mu => self%mu) + mu => self%mu, GMcb => system%cb%Gmass) if (ntp == 0) return allocate(iflag(ntp)) @@ -98,7 +98,7 @@ module subroutine whm_drift_tp(self, cb, param, dt) do i = 1,ntp rmag = norm2(xh(:, i)) vmag2 = dot_product(vh(:, i), vh(:, i)) - energy = 0.5_DP * vmag2 - cb%Gmass / rmag + energy = 0.5_DP * vmag2 - GMcb / rmag dtp(i) = dt * (1.0_DP + 3 * param%inv_c2 * energy) end do else diff --git a/src/whm/whm_getacch.f90 b/src/whm/whm_getacch.f90 index 128915360..53cd9b97a 100644 --- a/src/whm/whm_getacch.f90 +++ b/src/whm/whm_getacch.f90 @@ -1,7 +1,7 @@ submodule(whm_classes) s_whm_getacch use swiftest contains - module subroutine whm_getacch_pl(self, cb, param, t) + module subroutine whm_getacch_pl(self, system, param, t) !! author: David A. Minton !! !! Compute heliocentric accelerations of planets @@ -10,16 +10,15 @@ module subroutine whm_getacch_pl(self, cb, param, t) !! Adapted from David E. Kaufmann's Swifter routine whm_getacch.f90 implicit none ! Arguments - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + class(whm_nbody_system), intent(inout) :: system !! Swiftest central body particle data structure + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of + real(DP), intent(in) :: t !! Current time ! Internals integer(I4B) :: i real(DP), dimension(NDIM) :: ah0 - associate(pl => self, npl => self%nbody, j2rp2 => cb%j2rp2, & - ah => self%ah, xh => self%xh, xj => self%xj, vh => self%vh, vj => self%vj) + associate(pl => self, cb => system%cb, npl => self%nbody) if (npl == 0) return call pl%set_ir3() @@ -33,13 +32,13 @@ module subroutine whm_getacch_pl(self, cb, param, t) if (param%loblatecb) call pl%obl_acc(cb) if (param%lextra_force) call pl%user_getacch(cb, param, t) - if (param%lgr) call pl%gr_getacch(cb, param) + if (param%lgr) call pl%gr_getacch(param) end associate return end subroutine whm_getacch_pl - module subroutine whm_getacch_tp(self, cb, pl, param, t, xh) + module subroutine whm_getacch_tp(self, system, param, t, xh) !! author: David A. Minton !! !! Compute heliocentric accelerations of test particles @@ -48,28 +47,27 @@ module subroutine whm_getacch_tp(self, cb, pl, param, t, xh) !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_tp.f90 implicit none ! Arguments - class(whm_tp), intent(inout) :: self !! WHM test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Generic Swiftest central body particle data structuree - class(whm_pl), intent(inout) :: pl !! Generic Swiftest massive body particle data structure. + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + class(whm_nbody_system), intent(inout) :: system !! Swiftest central body particle data structure + class(whm_pl), intent(inout) :: pl !! Generic Swiftest massive body particle data structure. class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets + real(DP), intent(in) :: t !! Current time + real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets ! Internals - integer(I4B) :: i - real(DP), dimension(NDIM) :: ah0 + integer(I4B) :: i + real(DP), dimension(NDIM) :: ah0 - associate(tp => self, ntp => self%nbody, npl => pl%nbody, j2rp2 => cb%j2rp2, aht => self%ah, & - ir3h => pl%ir3h, GMpl => pl%Gmass) + associate(tp => self, ntp => self%nbody, pl => system%pl, npl => system%pl%nbody) if (ntp == 0 .or. npl == 0) return ah0 = whm_getacch_ah0(pl%Gmass(:), xh(:,:), npl) do i = 1, ntp tp%ah(:, i) = ah0(:) end do - call whm_getacch_ah3_tp(cb, pl, tp, xh) + call whm_getacch_ah3_tp(system, xh) if (param%loblatecb) call tp%obl_acc(cb) if (param%lextra_force) call tp%user_getacch(cb, param, t) - if (param%lgr) call tp%gr_getacch(cb, param) + if (param%lgr) call tp%gr_getacch(param) end associate return end subroutine whm_getacch_tp @@ -198,7 +196,7 @@ pure subroutine whm_getacch_ah3(pl) return end subroutine whm_getacch_ah3 - pure subroutine whm_getacch_ah3_tp(cb, pl, tp, xh) + pure subroutine whm_getacch_ah3_tp(system, xh) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations of test particles @@ -207,16 +205,14 @@ pure subroutine whm_getacch_ah3_tp(cb, pl, tp, xh) !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah3.f90 implicit none ! Arguments - class(swiftest_cb), intent(in) :: cb !! Swiftest central body object - class(whm_pl), intent(in) :: pl !! WHM massive body object - class(whm_tp), intent(inout) :: tp !! WHM test particle object + class(whm_nbody_system) :: system !! WHM nbody system object real(DP), dimension(:,:), intent(in) :: xh !! Position vector of massive bodies at required point in step ! Internals integer(I4B) :: i, j real(DP) :: rji2, irij3, fac real(DP), dimension(NDIM) :: dx, acc - associate(ntp => tp%nbody, npl => pl%nbody, msun => cb%Gmass, GMpl => pl%Gmass, xht => tp%xh, aht => tp%ah) + associate(ntp => system%tp%nbody, npl => system%pl%nbody, msun => system%cb%Gmass, GMpl => system%pl%Gmass, xht => system%tp%xh, aht => system%tp%ah) if (ntp == 0) return do i = 1, ntp acc(:) = 0.0_DP diff --git a/src/whm/whm_gr.f90 b/src/whm/whm_gr.f90 index 073ddc2ae..47f4bd449 100644 --- a/src/whm/whm_gr.f90 +++ b/src/whm/whm_gr.f90 @@ -1,8 +1,7 @@ submodule(whm_classes) s_whm_gr use swiftest contains - module subroutine whm_gr_getacch_pl(self, cb, param) - !! author: David A. Minton + module subroutine whm_gr_getacch_pl(self, param) !! author: David A. Minton !! !! Compute relativisitic accelerations of massive bodies !! Based on Saha & Tremaine (1994) Eq. 28 @@ -10,8 +9,7 @@ module subroutine whm_gr_getacch_pl(self, cb, param) !! Adapted from David A. Minton's Swifter routine routine gr_whm_getacch.f90 implicit none ! Arguments - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of ! Internals integer(I4B) :: i @@ -19,11 +17,10 @@ module subroutine whm_gr_getacch_pl(self, cb, param) real(DP), dimension(:, :), allocatable :: aj real(DP) :: beta, rjmag4 - associate(n => self%nbody, msun => cb%Gmass, mu => self%muj, c2 => param%inv_c2, & + associate(n => self%nbody, mu => self%muj, c2 => param%inv_c2, & ah => self%ah, xj => self%xj, GMpl => self%Gmass, eta => self%eta) if (n == 0) return allocate(aj, mold = ah) - !do concurrent(i = 1:n) do i = 1, n rjmag4 = (dot_product(xj(:, i), xj(:, i)))**2 beta = - mu(i)**2 * c2 @@ -39,7 +36,7 @@ module subroutine whm_gr_getacch_pl(self, cb, param) return end subroutine whm_gr_getacch_pl - module subroutine whm_gr_getacch_tp(self, cb, param) + module subroutine whm_gr_getacch_tp(self, param) !! author: David A. Minton !! !! Compute relativisitic accelerations of test particles @@ -48,17 +45,15 @@ module subroutine whm_gr_getacch_tp(self, cb, param) !! Adapted from David A. Minton's Swifter routine routine gr_whm_getacch.f90 implicit none ! Arguments - class(whm_tp), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of + class(whm_tp), intent(inout) :: self !! WHM massive body particle data structure + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i real(DP) :: rjmag4, beta - associate(n => self%nbody, msun => cb%Gmass, mu => self%mu,& + associate(n => self%nbody, mu => self%mu,& c2 => param%inv_c2, ah => self%ah, xh => self%xh, status => self%status) if (n == 0) return - !do concurrent (i = 1:n, status(i) == active) do i = 1, n rjmag4 = (dot_product(xh(:, i), xh(:, i)))**2 beta = - mu(i)**2 * c2 @@ -77,15 +72,14 @@ module pure subroutine whm_gr_p4_pl(self, param, dt) !! Adapted from David A. Minton's Swifter routine routine gr_whm_p4.f90 implicit none ! Arguments - class(whm_pl), intent(inout) :: self !! Swiftest particle object + class(whm_pl), intent(inout) :: self !! Swiftest particle object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters - real(DP), intent(in) :: dt !! Step size + real(DP), intent(in) :: dt !! Step size ! Internals integer(I4B) :: i associate(n => self%nbody, xj => self%xj, vj => self%vj, status => self%status, c2 => param%inv_c2) if (n == 0) return - !do concurrent (i = 1:n, status(i) == ACTIVE) do i = 1,n call p4_func(xj(:, i), vj(:, i), dt, c2) end do diff --git a/src/whm/whm_step.f90 b/src/whm/whm_step.f90 index 05aeebde9..41d3394fa 100644 --- a/src/whm/whm_step.f90 +++ b/src/whm/whm_step.f90 @@ -2,7 +2,7 @@ use swiftest contains - module subroutine whm_step_system(self, param) + module subroutine whm_step_system(self, param, dt) !! author: David A. Minton !! !! Step massive bodies and and active test particles ahead in heliocentric coordinates @@ -11,9 +11,12 @@ module subroutine whm_step_system(self, param) !! Adapted from David E. Kaufmann's Swifter routine whm_step.f90 implicit none ! Arguments - class(whm_nbody_system), intent(inout) :: self !! WHM nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters + class(whm_nbody_system), intent(inout) :: self !! WHM nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of on parameters + integer(I4B), intent(in) :: dt !! Current stepsize + select type (system => self) + class is (whm_nbody_system) select type(cb => self%cb) class is (whm_cb) select type(pl => self%pl) @@ -23,19 +26,20 @@ module subroutine whm_step_system(self, param) associate(ntp => tp%nbody, npl => pl%nbody, t => param%t, dt => param%dt) call pl%set_rhill(cb) call tp%set_beg_end(xbeg = pl%xh) - call pl%step(cb, param, t, dt) + call pl%step(system, param, dt) if (ntp > 0) then call tp%set_beg_end(xend = pl%xh) - call tp%step(cb, pl, param, t, dt) + call tp%step(system, param, dt) end if end associate end select end select end select + end select return end subroutine whm_step_system - module subroutine whm_step_pl(self, cb, param, t, dt) + module subroutine whm_step_pl(self, system, param, dt) !! author: David A. Minton !! !! Step planets ahead using kick-drift-kick algorithm @@ -45,20 +49,18 @@ module subroutine whm_step_pl(self, cb, param, t, dt) !logical, save :: lfirst = .true. implicit none ! Arguments - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time - real(DP), intent(in) :: dt !! Stepsize + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters + integer(I4B), intent(in) :: dt !! Current stepsize ! Internals real(DP) :: dth - associate(pl => self, xh => self%xh, vh => self%vh, ah => self%ah, & - xj => self%xj, vj => self%vj) + associate(cb => system%cb, t => param%t) dth = 0.5_DP * dt if (pl%lfirst) then call pl%h2j(cb) - call pl%getacch(cb, param, t) + call pl%getacch(system, param, t) pl%lfirst = .false. end if @@ -69,13 +71,13 @@ module subroutine whm_step_pl(self, cb, param, t, dt) call pl%drift(cb, param, dt) if (param%lgr) call pl%gr_p4(param, dth) call pl%j2h(cb) - call pl%getacch(cb, param, t + dt) + call pl%getacch(system, param, t + dt) call pl%kickvh(dth) end associate return end subroutine whm_step_pl - module subroutine whm_step_tp(self, cb, pl, param, t, dt) + module subroutine whm_step_tp(self, system, param, dt) !! author: David A. Minton !! !! Step active test particles ahead using kick-drift-kick algorithm @@ -84,20 +86,17 @@ module subroutine whm_step_tp(self, cb, pl, param, t, dt) !! Adapted from David E. Kaufmann's Swifter routine whm_step_tp.f90 implicit none ! Arguments - class(whm_tp), intent(inout) :: self !! WHM test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure - class(whm_pl), intent(inout) :: pl !! WHM massive body data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time - real(DP), intent(in) :: dt !! Stepsize + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of parameters + integer(I4B), intent(in) :: dt !! Current stepsize ! Internals real(DP) :: dth - associate(tp => self, xht => self%xh, vht => self%vh, aht => self%ah, & - xbeg => self%xbeg, xend => self%xend) + associate(cb => system%cb, pl => system%pl, t => param%t, xbeg => self%xbeg, xend => self%xend) dth = 0.5_DP * dt if (tp%lfirst) then - call tp%getacch(cb, pl, param, t, xbeg) + call tp%getacch(system, param, t, xbeg) tp%lfirst = .false. end if call tp%kickvh(dth) @@ -105,7 +104,7 @@ module subroutine whm_step_tp(self, cb, pl, param, t, dt) if (param%lgr) call tp%gr_p4(param, dth) call tp%drift(cb, param, dt) if (param%lgr) call tp%gr_p4(param, dth) - call tp%getacch(cb, pl, param, t + dt, xend) + call tp%getacch(system, param, t + dt, xend) call tp%kickvh(dth) end associate return