diff --git a/src/io/io.f90 b/src/io/io.f90 index 3a5a931b6..a4b619847 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -189,7 +189,6 @@ module subroutine io_dump_system(self, param) dump_param%incbfile = trim(adjustl(DUMP_CB_FILE(idx))) dump_param%inplfile = trim(adjustl(DUMP_PL_FILE(idx))) dump_param%intpfile = trim(adjustl(DUMP_TP_FILE(idx))) - dump_param%out_form = XV dump_param%in_form = XV dump_param%out_stat = 'APPEND' dump_param%in_type = REAL8_TYPE @@ -197,6 +196,7 @@ module subroutine io_dump_system(self, param) call dump_param%dump(param_file_name) + dump_param%out_form = XV call self%cb%dump(dump_param) call self%pl%dump(dump_param) call self%tp%dump(dump_param) @@ -454,7 +454,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) call io_toupper(param_value) param%out_stat = param_value case ("ISTEP_DUMP") - read(param_value, *) param%istep_dump + read(param_value, *, err = 667, iomsg = iomsg) param%istep_dump case ("CHK_CLOSE") call io_toupper(param_value) if (param_value == "YES" .or. param_value == 'T') param%lclose = .true. @@ -551,6 +551,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) read(param_value, *, err = 667, iomsg = iomsg) param%Ecollisions case("EUNTRACKED") read(param_value, *, err = 667, iomsg = iomsg) param%Euntracked + case ("MAXID") + read(param_value, *, err = 667, iomsg = iomsg) param%maxid case ("NPLMAX", "NTPMAX", "GMTINY", "MIN_GMFRAG", "PARTICLE_OUT", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters case default write(iomsg,*) "Unknown parameter -> ",param_name @@ -760,6 +762,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) write(param_name, Afmt) "TP_IN"; write(param_value, Afmt) trim(adjustl(param%intpfile)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "IN_TYPE"; write(param_value, Afmt) trim(adjustl(param%in_type)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "IN_FORM"; write(param_value, Afmt) trim(adjustl(param%in_form)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) + if (param%istep_dump > 0) write(param_name, Afmt) "ISTEP_DUMP"; write(param_value, Ifmt) param%istep_dump; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) if (param%istep_out > 0) then write(param_name, Afmt) "ISTEP_OUT"; write(param_value, Ifmt) param%istep_out; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "BIN_OUT"; write(param_value, Afmt) trim(adjustl(param%outfile)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) @@ -768,9 +771,6 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) write(param_name, Afmt) "OUT_STAT"; write(param_value, Afmt) "APPEND"; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) end if write(param_name, Afmt) "ENC_OUT"; write(param_value, Afmt) trim(adjustl(param%enc_out)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) - if (param%istep_dump > 0) then - write(param_name, Afmt) "ISTEP_DUMP"; write(param_value, Ifmt) param%istep_dump; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) - end if write(param_name, Afmt) "CHK_RMIN"; write(param_value, Rfmt) param%rmin; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "CHK_RMAX"; write(param_value, Rfmt) param%rmax; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "CHK_EJECT"; write(param_value, Rfmt) param%rmaxu; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) @@ -810,6 +810,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) write(param_name, Afmt) "EUNTRACKED"; write(param_value, Rfmt) param%Euntracked; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) end if write(param_name, Afmt) "FIRSTKICK"; write(param_value, Lfmt) param%lfirstkick; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "MAXID"; write(param_value, Ifmt) param%maxid; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) iostat = 0 iomsg = "UDIO not implemented" @@ -1053,7 +1054,7 @@ module function io_read_frame_body(self, iu, param) result(ierr) read(iu, err = 667, iomsg = errmsg) pl%Gmass(:) pl%mass(:) = pl%Gmass(:) / param%GU if (param%lrhill_present) read(iu, err = 667, iomsg = errmsg) pl%rhill(:) - if (param%lclose) read(iu, err = 667, iomsg = errmsg) pl%radius(:) + read(iu, err = 667, iomsg = errmsg) pl%radius(:) if (param%lrotation) then read(iu, err = 667, iomsg = errmsg) pl%Ip(1, :) read(iu, err = 667, iomsg = errmsg) pl%Ip(2, :) @@ -1080,7 +1081,7 @@ module function io_read_frame_body(self, iu, param) result(ierr) end if self%Gmass(i) = real(val, kind=DP) self%mass(i) = real(val / param%GU, kind=DP) - if (param%lclose) read(iu, *, err = 667, iomsg = errmsg) self%radius(i) + read(iu, *, err = 667, iomsg = errmsg) self%radius(i) class is (swiftest_tp) read(iu, *, err = 667, iomsg = errmsg) self%id(i) end select @@ -1507,7 +1508,7 @@ module subroutine io_write_frame_body(self, iu, param) class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body write(iu, err = 667, iomsg = errmsg) pl%Gmass(1:n) if (param%lrhill_present) write(iu, err = 667, iomsg = errmsg) pl%rhill(1:n) - if (param%lclose) write(iu, err = 667, iomsg = errmsg) pl%radius(1:n) + write(iu, err = 667, iomsg = errmsg) pl%radius(1:n) if (param%lrotation) then write(iu, err = 667, iomsg = errmsg) pl%Ip(1, 1:n) write(iu, err = 667, iomsg = errmsg) pl%Ip(2, 1:n) diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index c351d0656..19d717b0d 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -12,37 +12,8 @@ module subroutine kick_getacch_int_pl(self) implicit none ! Arguments class(swiftest_pl), intent(inout) :: self - ! Internals - integer(I8B) :: k, nplpl - real(DP) :: rji2 - real(DP) :: dx, dy, dz - integer(I4B) :: i, j - real(DP), dimension(:,:), pointer :: ah, xh - real(DP), dimension(NDIM,self%nbody) :: ahi, ahj - integer(I4B), dimension(:,:), pointer :: k_plpl - logical, dimension(:), pointer :: lmask - real(DP), dimension(:), pointer :: Gmass - - associate(ah => self%ah, xh => self%xh, k_plpl => self%k_plpl, lmask => self%lmask, Gmass => self%Gmass) - nplpl = self%nplpl - ahi(:,:) = 0.0_DP - ahj(:,:) = 0.0_DP - !$omp parallel do default(shared)& - !$omp private(k, i, j, dx, dy, dz, rji2) & - !$omp reduction(+:ahi) & - !$omp reduction(-:ahj) - do k = 1_I8B, nplpl - i = k_plpl(1,k) - j = k_plpl(2,k) - dx = xh(1, j) - xh(1, i) - dy = xh(2, j) - xh(2, i) - dz = xh(3, j) - xh(3, i) - rji2 = dx**2 + dy**2 + dz**2 - if (lmask(i) .and. lmask(j)) call kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) - end do - !$omp end parallel do - ah(:,1:self%nbody) = ah(:,1:self%nbody) + ahi(:,1:self%nbody) + ahj(:,1:self%nbody) - end associate + + call kick_getacch_int_all_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) return end subroutine kick_getacch_int_pl @@ -63,7 +34,6 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) integer(I4B), intent(in) :: npl !! Number of active massive bodies ! Internals integer(I4B) :: i, j - !real(DP), dimension(:,:), allocatable :: aht if ((self%nbody == 0) .or. (npl == 0)) return @@ -73,60 +43,116 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) do j = 1, npl block real(DP) :: rji2 - real(DP) :: dx, dy, dz - dx = tp%xh(1, i) - xhp(1, j) - dy = tp%xh(2, i) - xhp(1, j) - dz = tp%xh(3, i) - xhp(1, j) - rji2 = dx**2 + dy**2 + dz**2 - call kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl(i), tp%ah(1,i), tp%ah(2,i), tp%ah(3,i)) + real(DP) :: xr, yr, zr + xr = tp%xh(1, i) - xhp(1, j) + yr = tp%xh(2, i) - xhp(1, j) + zr = tp%xh(3, i) - xhp(1, j) + rji2 = xr**2 + yr**2 + zr**2 + call kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl(i), tp%ah(1,i), tp%ah(2,i), tp%ah(3,i)) end block end do end if end do - !call move_alloc(aht, tp%ah) end associate return end subroutine kick_getacch_int_tp - module pure elemental subroutine kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) + + module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc) + !! author: David A. Minton + !! + !! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization + !! + !! Adapted from Hal Levison's Swift routine getacch_ah3.f + !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f9 implicit none - real(DP), intent(in) :: rji2 - real(DP), intent(in) :: dx, dy, dz - real(DP), intent(in) :: Gmi - real(DP), intent(in) :: Gmj - real(DP), intent(inout) :: axi, ayi, azi - real(DP), intent(inout) :: axj, ayj, azj + integer(I4B), intent(in) :: npl !! Number of massive bodies + integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute + integer(I4B), dimension(:,:), intent(in) :: k_plpl !! Array of interaction pair indices (flattened upper triangular matrix) + real(DP), dimension(:,:), intent(in) :: x !! Position vector array + real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass + real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array + ! Internals + integer(I8B) :: k + real(DP) :: rji2, rlim2 + real(DP) :: xr, yr, zr + integer(I4B) :: i, j + real(DP), dimension(NDIM,npl) :: ahi, ahj + + ahi(:,:) = 0.0_DP + ahj(:,:) = 0.0_DP + !$omp parallel do default(private)& + !$omp shared(nplpl, k_plpl, x, Gmass, radius) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do k = 1_I8B, nplpl + i = k_plpl(1,k) + j = k_plpl(2,k) + xr = x(1, j) - x(1, i) + yr = x(2, j) - x(2, i) + zr = x(3, j) - x(3, i) + rji2 = xr**2 + yr**2 + zr**2 + rlim2 = (radius(i) + radius(j))**2 + if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) + end do + !$omp end parallel do + acc(:,1:npl) = acc(:,1:npl) + ahi(:,1:npl) + ahj(:,1:npl) + return + end subroutine kick_getacch_int_all_pl + + + module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) + !! author: David A. Minton + !! + !! Compute direct cross (third) term heliocentric accelerations for a single pair of massive bodies + !! + !! Adapted from Hal Levison's Swift routine getacch_ah3.f + !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f9 + implicit none + real(DP), intent(in) :: rji2 !! Square of distance between the two bodies + real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions + real(DP), intent(in) :: Gmi !! G*mass of body i + real(DP), intent(in) :: Gmj !! G*mass of body j + real(DP), intent(inout) :: axi, ayi, azi !! Acceleration vector components of body i + real(DP), intent(inout) :: axj, ayj, azj !! Acceleration vector components of body j ! Internals real(DP) :: faci, facj, irij3 irij3 = 1.0_DP / (rji2 * sqrt(rji2)) faci = Gmi * irij3 facj = Gmj * irij3 - axi = axi + facj * dx - ayi = ayi + facj * dy - azi = azi + facj * dz - axj = axj - faci * dx - ayj = ayj - faci * dy - azj = azj - faci * dz + axi = axi + facj * xr + ayi = ayi + facj * yr + azi = azi + facj * zr + axj = axj - faci * xr + ayj = ayj - faci * yr + azj = azj - faci * zr + return end subroutine kick_getacch_int_one_pl - !module pure elemental subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl, ax, ay, az) - module pure subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl, ax, ay, az) + module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl, ax, ay, az) + !! author: David A. Minton + !! + !! Compute direct cross (third) term heliocentric accelerations of a single test particle massive body pair. + !! + !! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f + !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int_tp.f90 implicit none - real(DP), intent(in) :: rji2 - real(DP), intent(in) :: dx, dy, dz - real(DP), intent(in) :: GMpl - real(DP), intent(inout) :: ax, ay, az + real(DP), intent(in) :: rji2 !! Square of distance between the test particle and massive body + real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions + real(DP), intent(in) :: Gmpl !! G*mass of massive body + real(DP), intent(inout) :: ax, ay, az !! Acceleration vector components of test particle ! Internals real(DP) :: fac fac = GMpl / (rji2 * sqrt(rji2)) - ax = ax - fac * dx - ay = ay - fac * dy - az = az - fac * dz + ax = ax - fac * xr + ay = ay - fac * yr + az = az - fac * zr return end subroutine kick_getacch_int_one_tp diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 3b6c2d591..19c7967bd 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -15,8 +15,7 @@ module swiftest_classes !> Each paramter is initialized to a default values. type :: swiftest_parameters integer(I4B) :: integrator = UNKNOWN_INTEGRATOR !! Symbolic name of the nbody integrator used - integer(I4B) :: nplmax = -1 !! Maximum allowed number of massive bodies - integer(I4B) :: ntpmax = -1 !! Maximum allowed number of test particles + integer(I4B) :: maxid = -1 !! The current maximum particle id number real(DP) :: t0 = -1.0_DP !! Integration start time real(DP) :: t = -1.0_DP !! Integration current time real(DP) :: tstop = -1.0_DP !! Integration stop time @@ -295,7 +294,6 @@ module swiftest_classes logical :: lbeg !! True if this is the beginning of a step. This is used so that test particle steps can be calculated !! separately from massive bodies. Massive body variables are saved at half steps, and passed to !! the test particles - integer(I4B) :: maxid = -1 !! The current maximum particle id number contains !> Each integrator will have its own version of the step procedure(abstract_step_system), deferred :: step @@ -711,29 +709,40 @@ module subroutine io_write_frame_system(self, iu, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine io_write_frame_system - module pure elemental subroutine kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) - implicit none - real(DP), intent(in) :: rji2 - real(DP), intent(in) :: dx, dy, dz - real(DP), intent(in) :: Gmi - real(DP), intent(in) :: Gmj - real(DP), intent(inout) :: axi, ayi, azi - real(DP), intent(inout) :: axj, ayj, azj + module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc) + implicit none + integer(I4B), intent(in) :: npl !! Number of massive bodies + integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute + integer(I4B), dimension(:,:), intent(in) :: k_plpl !! Array of interaction pair indices (flattened upper triangular matrix) + real(DP), dimension(:,:), intent(in) :: x !! Position vector array + real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass + real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array + end subroutine kick_getacch_int_all_pl + + module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) + !$omp declare simd(kick_getacch_int_one_pl) + implicit none + real(DP), intent(in) :: rji2 !! Square of distance between the two bodies + real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions + real(DP), intent(in) :: Gmi !! G*mass of body i + real(DP), intent(in) :: Gmj !! G*mass of body j + real(DP), intent(inout) :: axi, ayi, azi !! Acceleration vector components of body i + real(DP), intent(inout) :: axj, ayj, azj !! Acceleration vector components of body j end subroutine kick_getacch_int_one_pl - !module pure elemental subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, Gmpl, ax, ay, az) - module pure subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, Gmpl, ax, ay, az) + module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, Gmpl, ax, ay, az) !$omp declare simd(kick_getacch_int_one_tp) implicit none - real(DP), intent(in) :: rji2 - real(DP), intent(in) :: dx, dy, dz - real(DP), intent(in) :: Gmpl - real(DP), intent(inout) :: ax, ay, az + real(DP), intent(in) :: rji2 !! Square of distance between the test particle and massive body + real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions + real(DP), intent(in) :: Gmpl !! G*mass of massive body + real(DP), intent(inout) :: ax, ay, az !! Acceleration vector components of test particle end subroutine kick_getacch_int_one_tp module subroutine kick_getacch_int_pl(self) implicit none - class(swiftest_pl), intent(inout) :: self + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object end subroutine kick_getacch_int_pl module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) @@ -1448,7 +1457,7 @@ end subroutine util_spill_tp module subroutine util_valid_id_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine util_valid_id_system module subroutine util_version() diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index b259dfba9..5d95ba44f 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -68,7 +68,7 @@ module function symba_collision_casedisruption(system, param, family, x, v, mass ! Distribute any residual mass if there is any and set the radius m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) - id_frag(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] + id_frag(:) = [(i, i = param%maxid + 1, param%maxid + nfrag)] do i = 1, nfrag Ip_frag(:, i) = Ip_new(:) @@ -169,7 +169,7 @@ module function symba_collision_casehitandrun(system, param, family, x, v, mass, m_frag(2:nfrag) = (mtot - m_frag(1)) / (nfrag - 1) rad_frag(2:nfrag) = (3 * m_frag(2:nfrag) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) - id_frag(1:nfrag) = [(i, i = system%maxid + 1, system%maxid + nfrag)] + id_frag(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] do i = 1, nfrag Ip_frag(:, i) = Ip(:, jproj) @@ -372,7 +372,7 @@ module function symba_collision_casesupercatastrophic(system, param, family, x, ! Distribute any residual mass if there is any and set the radius m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) - id_frag(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] + id_frag(:) = [(i, i = param%maxid + 1, param%maxid + nfrag)] do i = 1, nfrag Ip_frag(:, i) = Ip_new(:) @@ -821,7 +821,7 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, implicit none ! Arguments class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision integer(I4B), dimension(:), intent(in) :: id_frag !! List of fragment ids real(DP), dimension(:), intent(in) :: m_frag, rad_frag !! Distribution of fragment mass and radii @@ -861,7 +861,7 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, ! Copy over identification, information, and physical properties of the new bodies from the fragment list plnew%id(1:nfrag) = id_frag(1:nfrag) - system%maxid = system%maxid + nfrag + param%maxid = param%maxid + nfrag plnew%xb(:, 1:nfrag) = xb_frag(:, 1:nfrag) plnew%vb(:, 1:nfrag) = vb_frag(:, 1:nfrag) do i = 1, nfrag diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index d282f9c53..af440eada 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -12,38 +12,8 @@ module subroutine symba_kick_getacch_int_pl(self) implicit none ! Arguments class(symba_pl), intent(inout) :: self - ! Internals - integer(I8B) :: k, nplplm - real(DP) :: rji2, rlim2 - real(DP) :: dx, dy, dz - integer(I4B) :: i, j - real(DP), dimension(:,:), pointer :: ah, xh - real(DP), dimension(NDIM,self%nbody) :: ahi, ahj - integer(I4B), dimension(:,:), pointer :: k_plpl - logical, dimension(:), pointer :: lmask - real(DP), dimension(:), pointer :: Gmass, radius - associate(ah => self%ah, xh => self%xh, k_plpl => self%k_plpl, lmask => self%lmask, Gmass => self%Gmass, radius => self%radius) - nplplm = self%nplplm - ahi(:,:) = 0.0_DP - ahj(:,:) = 0.0_DP - !$omp parallel do default(shared)& - !$omp private(k, i, j, dx, dy, dz, rji2, rlim2) & - !$omp reduction(+:ahi) & - !$omp reduction(-:ahj) - do k = 1_I8B, nplplm - i = k_plpl(1,k) - j = k_plpl(2,k) - dx = xh(1, j) - xh(1, i) - dy = xh(2, j) - xh(2, i) - dz = xh(3, j) - xh(3, i) - rji2 = dx**2 + dy**2 + dz**2 - rlim2 = (radius(i)+radius(j))**2 - if ((rji2 > rlim2) .and. lmask(i) .and. lmask(j)) call kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) - end do - !$omp end parallel do - ah(:,1:self%nbody) = ah(:,1:self%nbody) + ahi(:,1:self%nbody) + ahj(:,1:self%nbody) - end associate + call kick_getacch_int_all_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) return end subroutine symba_kick_getacch_int_pl @@ -66,7 +36,7 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) ! Internals integer(I4B) :: i, j integer(I8B) :: k, nplplenc - real(DP) :: rji2, rlim2, dx, dy, dz + real(DP) :: rji2, rlim2, xr, yr, zr real(DP), dimension(NDIM,self%nbody) :: ahi, ahj class(symba_plplenc), pointer :: plplenc_list real(DP), dimension(:), pointer :: Gmass, radius @@ -82,18 +52,18 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) ahi(:,:) = 0.0_DP ahj(:,:) = 0.0_DP !$omp parallel do default(shared)& - !$omp private(k, i, j, dx, dy, dz, rji2, rlim2) & + !$omp private(k, i, j, xr, yr, zr, rji2, rlim2) & !$omp reduction(+:ahi) & !$omp reduction(-:ahj) do k = 1_I8B, nplplenc i = plplenc_list%index1(k) j = plplenc_list%index2(k) - dx = xh(1, j) - xh(1, i) - dy = xh(2, j) - xh(2, i) - dz = xh(3, j) - xh(3, i) - rji2 = dx**2 + dy**2 + dz**2 + xr = xh(1, j) - xh(1, i) + yr = xh(2, j) - xh(2, i) + zr = xh(3, j) - xh(3, i) + rji2 = xr**2 + yr**2 + zr**2 rlim2 = (radius(i)+radius(j))**2 - if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) + if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) end do !$omp end parallel do ah(:,1:self%nbody) = ah(:,1:self%nbody) - ahi(:,1:self%nbody) - ahj(:,1:self%nbody) diff --git a/src/util/util_valid.f90 b/src/util/util_valid.f90 index f05c81f35..eb6df7577 100644 --- a/src/util/util_valid.f90 +++ b/src/util/util_valid.f90 @@ -12,7 +12,7 @@ module subroutine util_valid_id_system(self, param) implicit none ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i integer(I4B), dimension(:), allocatable :: idarr @@ -34,7 +34,7 @@ module subroutine util_valid_id_system(self, param) call util_exit(FAILURE) end if end do - self%maxid = maxval(idarr) + param%maxid = max(param%maxid, maxval(idarr)) end associate return