From f8c3decaac979afa9da0b4e38e74b9816e24f22d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 6 Jul 2021 09:17:17 -0400 Subject: [PATCH] Finished major restructuring of data types. It now compiles. --- Makefile | 18 ++++-- src/{driftkick => drift}/drift.f90 | 0 src/helio/helio_getacch.f90 | 44 +++++++-------- src/helio/helio_step.f90 | 87 +++++++++++++++-------------- src/io/io.f90 | 8 +-- src/{driftkick => kick}/kick.f90 | 0 src/modules/swiftest_classes.f90 | 88 ++++++++++++++++++------------ src/obl/obl.f90 | 35 ++++++------ src/whm/whm_getacch.f90 | 39 ++++++------- src/whm/whm_step.f90 | 6 +- 10 files changed, 171 insertions(+), 154 deletions(-) rename src/{driftkick => drift}/drift.f90 (100%) rename src/{driftkick => kick}/kick.f90 (100%) diff --git a/Makefile b/Makefile index 9bda6503a..dfc6fd5c0 100644 --- a/Makefile +++ b/Makefile @@ -86,7 +86,7 @@ lib: ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ make libdir - cd $(SWIFTEST_HOME)/src/driftkick; \ + cd $(SWIFTEST_HOME)/src/drift; \ rm -f Makefile.Defines Makefile; \ ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ @@ -106,6 +106,11 @@ lib: ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ make libdir + cd $(SWIFTEST_HOME)/src/kick; \ + rm -f Makefile.Defines Makefile; \ + ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ + ln -s $(SWIFTEST_HOME)/Makefile .; \ + make libdir cd $(SWIFTEST_HOME)/src/obl; \ rm -f Makefile.Defines Makefile; \ ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ @@ -179,20 +184,21 @@ bin: *.f90 clean: cd $(SWIFTEST_HOME)/src/modules; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/discard; rm -f Makefile.Defines Makefile *.gc* - cd $(SWIFTEST_HOME)/src/driftkick; rm -f Makefile.Defines Makefile *.gc* + cd $(SWIFTEST_HOME)/src/drift; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/eucl; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/gr; rm -f Makefile.Defines Makefile *.gc* + cd $(SWIFTEST_HOME)/src/helio; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/io; rm -f Makefile.Defines Makefile *.gc* + cd $(SWIFTEST_HOME)/src/kick; rm -f Makefile.Defines Makefile *.gc* + cd $(SWIFTEST_HOME)/src/main; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/obl; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/operators; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/orbel; rm -f Makefile.Defines Makefile *.gc* + cd $(SWIFTEST_HOME)/src/rmvs; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/setup; rm -f Makefile.Defines Makefile *.gc* - cd $(SWIFTEST_HOME)/src/util; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/user; rm -f Makefile.Defines Makefile *.gc* - cd $(SWIFTEST_HOME)/src/main; rm -f Makefile.Defines Makefile *.gc* + cd $(SWIFTEST_HOME)/src/util; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/whm; rm -f Makefile.Defines Makefile *.gc* - cd $(SWIFTEST_HOME)/src/rmvs; rm -f Makefile.Defines Makefile *.gc* - cd $(SWIFTEST_HOME)/src/helio; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/bin; rm -f swiftest_* cd $(SWIFTEST_HOME)/bin; rm -f tool_* cd $(SWIFTEST_HOME)/lib; rm -f lib* diff --git a/src/driftkick/drift.f90 b/src/drift/drift.f90 similarity index 100% rename from src/driftkick/drift.f90 rename to src/drift/drift.f90 diff --git a/src/helio/helio_getacch.f90 b/src/helio/helio_getacch.f90 index 806047cc3..af6ab9e4d 100644 --- a/src/helio/helio_getacch.f90 +++ b/src/helio/helio_getacch.f90 @@ -21,16 +21,13 @@ module subroutine helio_getacch_pl(self, system, param, t) real(DP), dimension(:), allocatable, save :: irh real(DP), dimension(:, :), allocatable, save :: xh_loc, aobl - associate(pl => self, npl => self%nbody) - !if (lflag) then - pl%ahi(:,2:npl) = 0.0_DP - call helio_getacch_int_pl(pl, t) - !end if - !if (param%loblatecb) call self%obl_acc(cb) TODO: Fix this - !else + associate(cb => system%cb, pl => self, npl => self%nbody) + pl%ahi(:,2:npl) = 0.0_DP + call helio_getacch_int_pl(pl, t) pl%ah(:,:) = pl%ahi(:,:) - !end if + if (param%loblatecb) call pl%obl_acc(cb) if (param%lextra_force) call pl%user_getacch(system, param, t) + if (param%lgr) call pl%gr_getacch(param) end associate return @@ -56,18 +53,15 @@ module subroutine helio_getacch_tp(self, system, param, t, xhp) real(DP) :: r2, mu real(DP), dimension(:), allocatable, save :: irh, irht - ! executable code - associate(tp => self, ntp => self%nbody, npl => system%pl%nbody) + associate(tp => self, ntp => self%nbody, cb => system%cb, pl => system%pl, npl => system%pl%nbody) select type(pl => system%pl) - class is (rmvs_pl) - !if (lflag) then - self%ahi(:,:) = 0.0_DP - call helio_getacch_int_tp(tp, pl, t, xhp) - !end if - !if (param%loblatecb) call self%obl_acc(cb) TODO: Fix this + class is (helio_pl) + self%ahi(:,:) = 0.0_DP + call helio_getacch_int_tp(tp, pl, t, xhp) tp%ah(:,:) = tp%ahi(:,:) + if (param%loblatecb) call tp%obl_acc(cb) if (param%lextra_force) call tp%user_getacch(system, param, t) - if (param%lgr) call tp%gr_getacch(system, param) + if (param%lgr) call tp%gr_getacch(param) end select end associate return @@ -95,10 +89,10 @@ subroutine helio_getacch_int_pl(pl, t) dx(:) = pl%xh(:,j) - pl%xh(:,i) rji2 = dot_product(dx(:), dx(:)) irij3 = 1.0_DP / (rji2 * sqrt(rji2)) - faci = self%Gmass(i) * irij3 - facj = self%Gmass(j) * irij3 - self%ahi(:,i) = self%ahi(:,i) + facj * dx(:) - self%ahi(:,i) = self%ahi(:,j) - faci * dx(:) + faci = pl%Gmass(i) * irij3 + facj = pl%Gmass(j) * irij3 + pl%ahi(:,i) = pl%ahi(:,i) + facj * dx(:) + pl%ahi(:,i) = pl%ahi(:,j) - faci * dx(:) end do end do end associate @@ -115,10 +109,10 @@ subroutine helio_getacch_int_tp(tp, pl, t, xhp) !! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f implicit none ! Arguments - class(helio_tp), intent(inout) :: tp !! Helio test particle data structure - class(helio_pl), intent(inout) :: pl !! Helio massive body particle data structure - real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xhp !! Heliocentric positions of planets + class(helio_tp), intent(inout) :: tp !! Helio test particle data structure + class(swiftest_pl), intent(inout) :: pl !! Helio massive body particle data structure + real(DP), intent(in) :: t !! Current time + real(DP), dimension(:,:), intent(in) :: xhp !! Heliocentric positions of planets ! Internals integer(I4B) :: i, j real(DP) :: r2, fac diff --git a/src/helio/helio_step.f90 b/src/helio/helio_step.f90 index 480de821d..53052fc9e 100644 --- a/src/helio/helio_step.f90 +++ b/src/helio/helio_step.f90 @@ -17,22 +17,22 @@ module subroutine helio_step_system(self, param, t, dt) select type(system => self) class is (helio_nbody_system) - select type(cb => self%cb) - class is (helio_cb) - select type(pl => self%pl) - class is (helio_pl) - select type(tp => self%tp) - class is (helio_tp) - call pl%set_rhill(cb) - call tp%set_beg_end(xbeg = pl%xh) - call pl%step(system, param, t, dt) - if (ntp > 0) then - call tp%set_beg_end(xend = pl%xh) - call tp%step(system, param, t, dt) - end if - end select - end select - end select + select type(cb => self%cb) + class is (helio_cb) + select type(pl => self%pl) + class is (helio_pl) + select type(tp => self%tp) + class is (helio_tp) + call pl%set_rhill(cb) + call system%set_beg_end(xbeg = pl%xh) + call pl%step(system, param, t, dt) + if (tp%nbody > 0) then + call system%set_beg_end(xend = pl%xh) + call tp%step(system, param, t, dt) + end if + end select + end select + end select end select return end subroutine helio_step_system @@ -57,21 +57,21 @@ module subroutine helio_step_pl(self, system, param, t, dt) real(DP), dimension(NDIM) :: ptbeg, ptend !! TODO: Incorporate these into the tp structure logical, save :: lfirst = .true. - associate(cb => system%cb) + associate(pl => self, cb => system%cb) dth = 0.5_DP * dt if (lfirst) then - call self%vh2vb(cb) + call pl%vh2vb(cb) lfirst = .false. end if - call self%lindrift(cb, dth, ptbeg) - call self%getacch(system, param, t) - call self%kickvb(dth) + call pl%lindrift(cb, dth, ptbeg) + call pl%getacch(system, param, t) + call pl%kickvb(dth) - call self%drift(system, param, dt) - call self%getacch(system, param, t + dt) - call self%kickvb(dth) - call self%lindrift(cb, dth, ptend) - call self%vb2vh(cb) + call pl%drift(system, param, dt) + call pl%getacch(system, param, t + dt) + call pl%kickvb(dth) + call pl%lindrift(cb, dth, ptend) + call pl%vb2vh(cb) end associate return @@ -79,6 +79,7 @@ module subroutine helio_step_pl(self, system, param, t, dt) end subroutine helio_step_pl module subroutine helio_step_tp(self, system, param, t, dt) + !! author: David A. Minton !! !! Step active test particles ahead using Democratic Heliocentric method @@ -96,22 +97,24 @@ module subroutine helio_step_tp(self, system, param, t, dt) logical, save :: lfirst = .true. !! Flag to indicate that this is the first call real(DP) :: dth !! Half step size - ! executable code - associate(cb => system%cb, pl => system%pl) - dth = 0.5_DP * dt - if (lfirst) then - call self%vh2vb(vbcb = -self%ptbeg) - lfirst = .false. - end if - call self%lindrift(dth, self%ptbeg) - call self%getacch(system, param, t, self%xbeg) - call self%kickvb(dth) - call self%drift(system, param, dt) - call self%getacch(system, param, t + dt, self%xend) - call self%kickvb(dth) - call self%lindrift(dth, self%ptend) - call self%vb2vh(vbcb = -self%ptend) - end associate + select type(system) + class is (helio_nbody_system) + associate(tp => self, cb => system%cb, pl => system%pl, xbeg => system%xbeg, xend => system%xend) + dth = 0.5_DP * dt + if (lfirst) then + call tp%vh2vb(vbcb = -tp%ptbeg) + lfirst = .false. + end if + call tp%lindrift(dth, tp%ptbeg) + call tp%getacch(system, param, t, xbeg) + call tp%kickvb(dth) + call tp%drift(system, param, dt) + call tp%getacch(system, param, t + dt, xend) + call tp%kickvb(dth) + call tp%lindrift(dth, tp%ptend) + call tp%vb2vh(vbcb = -tp%ptend) + end associate + end select return diff --git a/src/io/io.f90 b/src/io/io.f90 index 629a4952f..c1774f7ca 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -517,13 +517,13 @@ end subroutine io_dump_system module function io_get_args(integrator, param_file_name) result(ierr) !! author: David A. Minton !! - !! Reads in the name of the parameters file. + !! Reads in the name of the parameter file from command line arguments. implicit none ! Arguments integer(I4B) :: integrator !! Symbolic code of the requested integrator character(len=:), allocatable :: param_file_name !! Name of the input parameters file ! Result - integer(I4B) :: ierr !! I/O error cod + integer(I4B) :: ierr !! I/O error code ! Internals character(len=STRMAX) :: arg1, arg2 integer :: narg,ierr_arg1, ierr_arg2 @@ -574,7 +574,7 @@ module function io_get_args(integrator, param_file_name) result(ierr) if (ierr /= 0) call util_exit(USAGE) end function io_get_args - module function io_get_token(buffer, ifirst, ilast, ierr) result(token) + function io_get_token(buffer, ifirst, ilast, ierr) result(token) !! author: David A. Minton !! !! Retrieves a character token from an input string. Here a token is defined as any set of contiguous non-blank characters not @@ -1126,7 +1126,7 @@ module subroutine io_write_discard(self, param, discards) end subroutine io_write_discard - subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & + module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & xh1, xh2, vh1, vh2, encounter_file, out_type) !! author: David A. Minton !! diff --git a/src/driftkick/kick.f90 b/src/kick/kick.f90 similarity index 100% rename from src/driftkick/kick.f90 rename to src/kick/kick.f90 diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 150408537..f7ed5beb2 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -5,7 +5,22 @@ module swiftest_classes !! Adapted from David E. Kaufmann's Swifter routine: module_swifter.f90 use swiftest_globals implicit none - public + private + public :: discard_pl_tp, discard_sun_tp, discard_system + public :: drift_one + public :: eucl_dist_index_plpl, eucl_dist_index_pltp, eucl_irij3_plpl + public :: kick_vb_body, kick_vh_body + public :: io_dump_param, io_dump_swiftest, io_dump_system, io_get_args, io_param_reader, io_param_writer, io_read_body_in, & + io_read_cb_in, io_read_param_in, io_read_frame_body, io_read_frame_cb, io_read_frame_system, io_read_initialize_system, & + io_write_discard, io_write_encounter, io_write_frame_body, io_write_frame_cb, io_write_frame_system + public :: obl_acc_body + public :: orbel_el2xv_vec, orbel_xv2el_vec, orbel_scget, orbel_xv2aeq, orbel_xv2aqt + public :: setup_body, setup_construct_system, setup_pl, setup_set_ir3h, setup_set_msys, setup_set_mu_pl, setup_set_mu_tp, & + setup_set_rhill, setup_tp + public :: user_getacch_body + public :: util_coord_b2h_pl, util_coord_b2h_tp, util_coord_h2b_pl, util_coord_h2b_tp, util_copy_body, util_copy_cb, util_copy_pl, & + util_copy_tp, util_copy_system, util_fill_body, util_fill_pl, util_fill_tp, util_reverse_status, util_spill_body, & + util_spill_pl, util_spill_tp !******************************************************************************************************************************** ! swiftest_parameters class definitions @@ -62,10 +77,11 @@ module swiftest_classes logical :: lyarkovsky = .false. !! Turn on Yarkovsky effect logical :: lyorp = .false. !! Turn on YORP effect contains - procedure :: reader => io_param_reader - procedure :: writer => io_param_writer - procedure :: dump => io_dump_param - procedure :: read_from_file => io_read_param_in + private + procedure, public :: reader => io_param_reader + procedure, public :: writer => io_param_writer + procedure, public :: dump => io_dump_param + procedure, public :: read_from_file => io_read_param_in !TODO: Figure out if user-defined derived-type io can be made to work properly !generic :: read(FORMATTED) => param_reader !generic :: write(FORMATTED) => param_writer @@ -379,28 +395,6 @@ module subroutine kick_vh_body(self, dt) real(DP), intent(in) :: dt !! Stepsize 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(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) - 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) implicit none class(swiftest_parameters),intent(in) :: self !! Output collection of parameters @@ -425,16 +419,30 @@ module function io_get_args(integrator, param_file_name) result(ierr) implicit none integer(I4B) :: integrator !! Symbolic code of the requested integrator character(len=:), allocatable :: param_file_name !! Name of the input parameters file - integer(I4B) :: ierr !! I/O error code + integer(I4B) :: ierr !! I/O error code end function io_get_args - module function io_get_token(buffer, ifirst, ilast, ierr) result(token) - character(len=*), intent(in) :: buffer !! Input string buffer - integer(I4B), intent(inout) :: ifirst !! Index of the buffer at which to start the search for a token - integer(I4B), intent(out) :: ilast !! Index of the buffer at the end of the returned token - integer(I4B), intent(out) :: ierr !! Error code - character(len=:),allocatable :: token !! Returned token string - end function io_get_token + module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) + implicit none + class(swiftest_parameters), intent(inout) :: 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(:) !! 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) + 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_read_body_in(self, param) implicit none @@ -494,6 +502,15 @@ module subroutine io_write_discard(self, param, discards) class(swiftest_body), intent(inout) :: discards !! Swiftest discard object end subroutine io_write_discard + module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & + xh1, xh2, vh1, vh2, encounter_file, out_type) + implicit none + integer(I4B), intent(in) :: name1, name2 + real(DP), intent(in) :: t, mass1, mass2, radius1, radius2 + real(DP), dimension(:), intent(in) :: xh1, xh2, vh1, vh2 + character(*), intent(in) :: encounter_file, out_type + end subroutine io_write_encounter + module subroutine io_write_frame_body(self, iu, param) implicit none class(swiftest_body), intent(in) :: self !! Swiftest particle object @@ -651,7 +668,6 @@ module subroutine util_copy_cb(self, src, mask) logical, dimension(:), intent(in) :: mask end subroutine util_copy_cb - module subroutine util_copy_pl(self, src, mask) implicit none class(swiftest_pl), intent(inout) :: self diff --git a/src/obl/obl.f90 b/src/obl/obl.f90 index 1a6b77d09..6d1487872 100644 --- a/src/obl/obl.f90 +++ b/src/obl/obl.f90 @@ -1,7 +1,7 @@ submodule (swiftest_classes) s_obl use swiftest contains - module procedure obl_acc_body + module subroutine obl_acc_body(self, cb) !! author: David A. Minton !! !! Compute the barycentric accelerations of bodies due to the oblateness of the central body @@ -10,39 +10,40 @@ !! Adapted from David E. Kaufmann's Swifter routine: obl_acc.f90 and obl_acc_tp.f90 !! Adapted from Hal Levison's Swift routine obl_acc.f and obl_acc_tp.f implicit none + ! Arguments + class(swiftest_body), intent(inout) :: self !! Swiftest generic body object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + ! Internals integer(I4B) :: i real(DP) :: r2, irh, rinv2, t0, t1, t2, t3, fac1, fac2 - associate(n => self%nbody, aobl => self%aobl, xh => self%xh, j2rp2 => cb%j2rp2, j4rp4 => cb%j4rp4, & - msun => cb%Gmass, aoblcb => cb%aobl, ah => self%ah) + associate(n => self%nbody) do i = 1, n - r2 = dot_product(xh(:, i), xh(:, i)) + r2 = dot_product(self%xh(:, i), self%xh(:, i)) irh = 1.0_DP / sqrt(r2) rinv2 = irh**2 - t0 = -msun * rinv2*rinv2*irh - t1 = 1.5_DP * j2rp2 - t2 = xh(3, i) * xh(3, i) * rinv2 - t3 = 1.875_DP * j4rp4 * rinv2 + t0 = -cb%Gmass * rinv2 * rinv2 * irh + t1 = 1.5_DP * cb%j2rp2 + t2 = self%xh(3, i) * self%xh(3, i) * rinv2 + t3 = 1.875_DP * cb%j4rp4 * rinv2 fac1 = t0 * (t1 - t3 - (5.0_DP * t1 - (14.0_DP - 21.0_DP * t2) * t3) * t2) fac2 = 2.0_DP * t0 * (t1 - (2.0_DP - (14.0_DP * t2 / 3.0_DP)) * t3) - aobl(:, i) = fac1 * xh(:, i) - aobl(3, i) = fac2 * xh(3, i) + aobl(3, i) + self%aobl(:, i) = fac1 * self%xh(:, i) + self%aobl(3, i) = fac2 * self%xh(3, i) + self%aobl(3, i) end do select type(self) class is (swiftest_pl) - associate(Mpl => self%Gmass) - do i = 1, NDIM - aoblcb(i) = -sum(Mpl(1:n) * aobl(i, 1:n)) / msun - end do - end associate + do i = 1, NDIM + cb%aobl(i) = -sum(self%Gmass(1:n) * self%aobl(i, 1:n)) / cb%Gmass + end do end select do i = 1, NDIM - ah(i, 1:n) = ah(i, 1:n) + aobl(i, 1:n) - aoblcb(i) + self%ah(i, 1:n) = self%ah(i, 1:n) + self%aobl(i, 1:n) - cb%aobl(i) end do end associate return - end procedure obl_acc_body + end subroutine obl_acc_body end submodule s_obl diff --git a/src/whm/whm_getacch.f90 b/src/whm/whm_getacch.f90 index d2e1bdc3e..7985fac5f 100644 --- a/src/whm/whm_getacch.f90 +++ b/src/whm/whm_getacch.f90 @@ -18,7 +18,7 @@ module subroutine whm_getacch_pl(self, system, param, t) integer(I4B) :: i real(DP), dimension(NDIM) :: ah0 - associate(pl => self, npl => self%nbody) + associate(cb => system%cb, pl => self, npl => self%nbody) if (npl == 0) return call pl%set_ir3() @@ -26,17 +26,14 @@ module subroutine whm_getacch_pl(self, system, param, t) do i = 1, npl pl%ah(:, i) = ah0(:) end do - select type(cb => system%cb) - class is (whm_cb) - call whm_getacch_ah1(cb, pl) - call whm_getacch_ah2(cb, pl) - call whm_getacch_ah3(pl) + call whm_getacch_ah1(cb, pl) + call whm_getacch_ah2(cb, pl) + call whm_getacch_ah3(pl) - if (param%loblatecb) call pl%obl_acc(cb) - if (param%lextra_force) call pl%user_getacch(system, param, t) - if (param%lgr) call pl%gr_getacch(param) - end select + if (param%loblatecb) call pl%obl_acc(cb) + if (param%lextra_force) call pl%user_getacch(system, param, t) + if (param%lgr) call pl%gr_getacch(param) end associate return @@ -63,7 +60,7 @@ module subroutine whm_getacch_tp(self, system, param, t, xhp) associate(tp => self, ntp => self%nbody, pl => system%pl, cb => system%cb, npl => system%pl%nbody) if (ntp == 0 .or. npl == 0) return - ah0 = whm_getacch_ah0(pl%Gmass(:), xhp(:,:), npl) + ah0(:) = whm_getacch_ah0(pl%Gmass(:), xhp(:,:), npl) do i = 1, ntp tp%ah(:, i) = ah0(:) end do @@ -110,8 +107,8 @@ pure subroutine whm_getacch_ah1(cb, pl) !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah1.f90 implicit none ! Arguments - class(whm_cb), intent(in) :: cb !! WHM central body object - class(whm_pl), intent(inout) :: pl !! WHM massive body object + class(swiftest_cb), intent(in) :: cb !! WHM central body object + class(whm_pl), intent(inout) :: pl !! WHM massive body object ! Internals integer(I4B) :: i real(DP), dimension(NDIM) :: ah1h, ah1j @@ -137,8 +134,8 @@ pure subroutine whm_getacch_ah2(cb, pl) !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah2.f90 implicit none ! Arguments - class(whm_cb), intent(in) :: cb !! Swiftest central body object - class(whm_pl), intent(inout) :: pl !! WHM massive body object + class(swiftest_cb), intent(in) :: cb !! Swiftest central body object + class(whm_pl), intent(inout) :: pl !! WHM massive body object ! Internals integer(I4B) :: i real(DP) :: etaj, fac @@ -169,15 +166,15 @@ pure subroutine whm_getacch_ah3(pl) !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah3.f90 implicit none - class(whm_pl), intent(inout) :: pl - integer(I4B) :: i, j - real(DP) :: rji2, irij3, faci, facj - real(DP), dimension(NDIM) :: dx - real(DP), dimension(:,:), allocatable :: ah3 + class(whm_pl), intent(inout) :: pl + integer(I4B) :: i, j + real(DP) :: rji2, irij3, faci, facj + real(DP), dimension(NDIM) :: dx + real(DP), dimension(:,:), allocatable :: ah3 associate(npl => pl%nbody) allocate(ah3, mold=pl%ah) - ah3(:, 1:npl) = 0.0_DP + ah3(:, :) = 0.0_DP do i = 1, npl - 1 do j = i + 1, npl diff --git a/src/whm/whm_step.f90 b/src/whm/whm_step.f90 index b37d06bbe..8e87796ea 100644 --- a/src/whm/whm_step.f90 +++ b/src/whm/whm_step.f90 @@ -85,10 +85,10 @@ module subroutine whm_step_tp(self, system, param, t, dt) select type(system) class is (whm_nbody_system) - associate(tp => self, cb => system%cb, pl => system%pl) + associate(tp => self, cb => system%cb, pl => system%pl, xbeg => system%xbeg, xend => system%xend) dth = 0.5_DP * dt if (tp%lfirst) then - call tp%getacch(system, param, t, system%xbeg) + call tp%getacch(system, param, t, xbeg) tp%lfirst = .false. end if call tp%kickvh(dth) @@ -96,7 +96,7 @@ module subroutine whm_step_tp(self, system, param, t, dt) if (param%lgr) call tp%gr_p4(param, dth) call tp%drift(system, param, dt) if (param%lgr) call tp%gr_p4(param, dth) - call tp%getacch(system, param, t + dt, system%xend) + call tp%getacch(system, param, t + dt, xend) call tp%kickvh(dth) end associate end select