diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index b060d13ef..3c2151413 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -409,7 +409,8 @@ def make_swiftest_labels(param): tlab.append('capm') plab = tlab.copy() plab.append('Gmass') - plab.append('radius') + if param['CHK_CLOSE'] == 'YES': + plab.append('radius') if param['RHILL_PRESENT'] == 'YES': plab.append('rhill') clab = ['Gmass', 'radius', 'J_2', 'J_4'] @@ -817,7 +818,10 @@ def swiftest_xr2infile(ds, param, framenum=-1): tp = frame.where(np.isnan(frame['Gmass']), drop=True).drop_vars(['Gmass', 'radius', 'J_2', 'J_4']) GMSun = np.double(cb['Gmass']) - RSun = np.double(cb['radius']) + if param['CHK_CLOSE'] == 'YES': + RSun = np.double(cb['radius']) + else: + RSun = param['CHK_RMIN'] J2 = np.double(cb['J_2']) J4 = np.double(cb['J_4']) cbname = cb['name'].values[0] @@ -851,7 +855,8 @@ def swiftest_xr2infile(ds, param, framenum=-1): print(pli['name'].values, pli['Gmass'].values, pli['rhill'].values, file=plfile) else: print(pli['name'].values, pli['Gmass'].values, file=plfile) - print(pli['radius'].values, file=plfile) + if param['CHK_CLOSE'] == 'YES': + print(pli['radius'].values, file=plfile) if param['IN_FORM'] == 'XV': print(pli['xhx'].values, pli['xhy'].values, pli['xhz'].values, file=plfile) print(pli['vhx'].values, pli['vhy'].values, pli['vhz'].values, file=plfile) @@ -918,7 +923,8 @@ def swiftest_xr2infile(ds, param, framenum=-1): else: print(f"{param['IN_FORM']} is not a valid input format type.") Gmass = pl['Gmass'].values - radius = pl['radius'].values + if param['CHK_CLOSE'] == 'YES': + radius = pl['radius'].values plfile.write_record(npl) plfile.write_record(plid) @@ -931,7 +937,8 @@ def swiftest_xr2infile(ds, param, framenum=-1): plfile.write_record(Gmass) if param['RHILL_PRESENT'] == 'YES': plfile.write_record(pl['rhill'].values) - plfile.write_record(radius) + if param['CHK_CLOSE'] == 'YES': + plfile.write_record(radius) if param['ROTATION'] == 'YES': plfile.write_record(pl['Ip1'].values) plfile.write_record(pl['Ip2'].values) @@ -995,7 +1002,10 @@ def swifter_xr2infile(ds, param, framenum=-1): tp = frame.where(np.isnan(frame['Gmass']), drop=True).drop_vars(['Gmass', 'radius', 'J_2', 'J_4']) GMSun = np.double(cb['Gmass']) - RSun = np.double(cb['radius']) + if param['CHK_CLOSE'] == 'YES': + RSun = np.double(cb['radius']) + else: + RSun = param['CHK_RMIN'] param['J2'] = np.double(cb['J_2']) param['J4'] = np.double(cb['J_4']) diff --git a/src/io/io.f90 b/src/io/io.f90 index f225c6523..99e33a7db 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -1525,7 +1525,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(:) - read(iu, err = 667, iomsg = errmsg) pl%radius(:) + if (param%lclose) 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, :) @@ -1553,7 +1553,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) - read(iu, *, err = 667, iomsg = errmsg) self%radius(i) + if (param%lclose) read(iu, *, err = 667, iomsg = errmsg) self%radius(i) class is (swiftest_tp) read(iu, *, err = 667, iomsg = errmsg) name(i) end select @@ -2027,7 +2027,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) - write(iu, err = 667, iomsg = errmsg) pl%radius(1:n) + if (param%lclose) 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 59b837408..d739427c4 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -33,9 +33,17 @@ module subroutine kick_getacch_int_pl(self, param) end if if (param%lflatten_interactions) then - call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + if (param%lclose) then + call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + else + call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, acc=self%ah) + end if else - call kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%xh, self%Gmass, self%radius, self%ah) + if (param%lclose) then + call kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%xh, self%Gmass, self%radius, self%ah) + else + call kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%xh, self%Gmass, acc=self%ah) + end if end if if (param%ladaptive_interactions .and. self%nplpl > 0) then @@ -78,13 +86,13 @@ module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, rad !! 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 - 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 + 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), optional :: radius !! Array of massive body radii + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array ! Internals integer(I8B) :: k real(DP), dimension(NDIM,npl) :: ahi, ahj @@ -95,23 +103,42 @@ module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, rad ahi(:,:) = 0.0_DP ahj(:,:) = 0.0_DP - !$omp parallel do default(private) schedule(static)& - !$omp shared(nplpl, k_plpl, x, Gmass, radius) & - !$omp lastprivate(rji2, rlim2, xr, yr, zr) & - !$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 + if (present(radius)) then + !$omp parallel do default(private) schedule(static)& + !$omp shared(nplpl, k_plpl, x, Gmass, radius) & + !$omp lastprivate(rji2, rlim2, xr, yr, zr) & + !$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 + else + !$omp parallel do default(private) schedule(static)& + !$omp shared(nplpl, k_plpl, x, Gmass, radius) & + !$omp lastprivate(rji2, xr, yr, zr) & + !$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 + 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 + end if do concurrent(i = 1:npl) acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i) @@ -130,12 +157,12 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius !! 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 - integer(I4B), intent(in) :: npl !! Total number of massive bodies - integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies - 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 + integer(I4B), intent(in) :: npl !! Total number of massive bodies + integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies + 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), optional :: radius !! Array of massive body radii + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array ! Internals real(DP), dimension(NDIM,npl) :: ahi, ahj integer(I4B) :: i, j @@ -145,23 +172,42 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius ahi(:,:) = 0.0_DP ahj(:,:) = 0.0_DP - !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, nplm, x, Gmass, radius) & - !$omp lastprivate(rji2, rlim2, xr, yr, zr) & - !$omp reduction(+:ahi) & - !$omp reduction(-:ahj) - do i = 1, nplm - do concurrent(j = i+1:npl) - 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)) + if (present(radius)) then + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl, nplm, x, Gmass, radius) & + !$omp lastprivate(rji2, rlim2, xr, yr, zr) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do i = 1, nplm + do concurrent(j = i+1:npl) + 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 end do - end do - !$omp end parallel do + !$omp end parallel do + else + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl, nplm, x, Gmass, radius) & + !$omp lastprivate(rji2, xr, yr, zr) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do i = 1, nplm + do concurrent(j = i+1:npl) + 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 + 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 + end do + !$omp end parallel do + end if do concurrent(i = 1:npl) acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 643d045f9..46e2f3824 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -906,23 +906,23 @@ end subroutine kick_getacch_int_tp module subroutine kick_getacch_int_all_flat_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 + 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), optional :: radius !! Array of massive body radii + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array end subroutine kick_getacch_int_all_flat_pl module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc) implicit none - integer(I4B), intent(in) :: npl !! Total number of massive bodies - integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies - 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 + integer(I4B), intent(in) :: npl !! Total number of massive bodies + integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies + 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), optional :: radius !! Array of massive body radii + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array end subroutine kick_getacch_int_all_triangular_pl module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index f8363d976..5165b2c8c 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -567,6 +567,9 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr) call check( nf90_get_var(iu%ncid, iu%radius_varid, rtemp, start=[1, tslot]) ) cb%radius = rtemp(1) pl%radius(:) = pack(rtemp, plmask) + else + cb%radius = param%rmin + pl%radius(:) = 0.0_DP end if if (param%lrotation) then @@ -896,12 +899,8 @@ module subroutine netcdf_write_frame_base(self, iu, param) select type(self) class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body call check( nf90_put_var(iu%ncid, iu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]) ) - if (param%lrhill_present) then - call check( nf90_put_var(iu%ncid, iu%rhill_varid, self%rhill(j), start=[idslot, tslot]) ) - end if - if (param%lclose) then - call check( nf90_put_var(iu%ncid, iu%radius_varid, self%radius(j), start=[idslot, tslot]) ) - end if + if (param%lrhill_present) call check( nf90_put_var(iu%ncid, iu%rhill_varid, self%rhill(j), start=[idslot, tslot]) ) + if (param%lclose) call check( nf90_put_var(iu%ncid, iu%radius_varid, self%radius(j), start=[idslot, tslot]) ) if (param%lrotation) then call check( nf90_put_var(iu%ncid, iu%Ip1_varid, self%Ip(1, j), start=[idslot, tslot]) ) call check( nf90_put_var(iu%ncid, iu%Ip2_varid, self%Ip(2, j), start=[idslot, tslot]) ) @@ -923,7 +922,7 @@ module subroutine netcdf_write_frame_base(self, iu, param) call check( nf90_put_var(iu%ncid, iu%id_varid, self%id, start=[idslot]) ) call check( nf90_put_var(iu%ncid, iu%Gmass_varid, self%Gmass, start=[idslot, tslot]) ) - call check( nf90_put_var(iu%ncid, iu%radius_varid, self%radius, start=[idslot, tslot]) ) + if (param%lclose) call check( nf90_put_var(iu%ncid, iu%radius_varid, self%radius, start=[idslot, tslot]) ) if (param%lrotation) then call check( nf90_put_var(iu%ncid, iu%Ip1_varid, self%Ip(1), start=[idslot, tslot]) ) call check( nf90_put_var(iu%ncid, iu%Ip2_varid, self%Ip(2), start=[idslot, tslot]) ) diff --git a/src/util/util_dealloc.f90 b/src/util/util_dealloc.f90 index bd648212e..09f1f29bf 100644 --- a/src/util/util_dealloc.f90 +++ b/src/util/util_dealloc.f90 @@ -22,8 +22,8 @@ module subroutine util_dealloc_body(self) if (allocated(self%vb)) deallocate(self%vb) if (allocated(self%ah)) deallocate(self%ah) if (allocated(self%aobl)) deallocate(self%aobl) - if (allocated(self%agr)) deallocate(self%lmask) - if (allocated(self%atide)) deallocate(self%lmask) + if (allocated(self%agr)) deallocate(self%agr) + if (allocated(self%atide)) deallocate(self%atide) if (allocated(self%ir3h)) deallocate(self%ir3h) if (allocated(self%a)) deallocate(self%a) if (allocated(self%e)) deallocate(self%e)