From 149189f6465dfa571cdbfc88c8b032eddb15ff81 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 16 Dec 2022 08:38:59 -0500 Subject: [PATCH] Added new .unit. operator for computing unit vectors --- src/CMakeLists.txt | 1 + src/fraggle/fraggle_generate.f90 | 6 -- src/modules/swiftest_operators.f90 | 71 +++++++++++++--- src/operators/operator_cross.f90 | 1 + src/operators/operator_mag.f90 | 1 + src/operators/operator_unit.f90 | 129 +++++++++++++++++++++++++++++ 6 files changed, 191 insertions(+), 18 deletions(-) create mode 100644 src/operators/operator_unit.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 594850a50..c74ea07a0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -50,6 +50,7 @@ SET(FAST_MATH_FILES ${SRC}/obl/obl.f90 ${SRC}/operators/operator_cross.f90 ${SRC}/operators/operator_mag.f90 + ${SRC}/operators/operator_unit.f90 ${SRC}/orbel/orbel.f90 ${SRC}/rmvs/rmvs_discard.f90 ${SRC}/rmvs/rmvs_encounter_check.f90 diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index bd1df3227..ac5f166ba 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -57,12 +57,6 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa return end if - ! Get the unit vectors for the relative position and velocity vectors. These are used to shift the fragment cloud depending on the - runit(:) = colliders%rb(:,2) - colliders%rb(:,1) - runit(:) = runit(:) / (.mag. runit(:)) - - vunit(:) = colliders%vb(:,2) - colliders%vb(:,1) - vunit(:) = vunit(:) / (.mag. vunit(:)) ! This is a factor that will "distort" the shape of the frgment cloud in the direction of the impact velocity f_spin= .mag. (runit(:) .cross. vunit(:)) diff --git a/src/modules/swiftest_operators.f90 b/src/modules/swiftest_operators.f90 index 30e5b26a6..165c7b283 100644 --- a/src/modules/swiftest_operators.f90 +++ b/src/modules/swiftest_operators.f90 @@ -1,7 +1,7 @@ !! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh !! This file is part of Swiftest. !! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! as published by the Free Software Foundation, either version NDIM of the License, or (at your option) any later version. !! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty !! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. !! You should have received a copy of the GNU General Public License along with Swiftest. @@ -11,15 +11,15 @@ module swiftest_operators !! author: David A. Minton !! !! Custom operators, including - !! A .cross. B = Cross product of A(1:3) and B(1:3) + !! A .cross. B = Cross product of A(1:NDIM) and B(1:NDIM) !! - !! Each operator can also do element-wise computation on arrays of the form .mag. A(1:3, 1:n) + !! Each operator can also do element-wise computation on arrays of the form .mag. A(1:NDIM, 1:n) use swiftest_globals implicit none public !******************************************************************************************************************************** - ! Interfaces for .cross. operator + ! Interfaces for .cross. operator: Computes the cross product of two (NDIM) vectors or (NDIM,:) arrays !******************************************************************************************************************************** interface operator(.cross.) @@ -27,49 +27,49 @@ pure module function operator_cross_sp(A, B) result(C) !$omp declare simd(operator_cross_sp) implicit none real(SP), dimension(:), intent(in) :: A, B - real(SP), dimension(3) :: C + real(SP), dimension(NDIM) :: C end function operator_cross_sp pure module function operator_cross_dp(A, B) result(C) !$omp declare simd(operator_cross_dp) implicit none real(DP), dimension(:), intent(in) :: A, B - real(DP), dimension(3) :: C + real(DP), dimension(NDIM) :: C end function operator_cross_dp pure module function operator_cross_qp(A, B) result(C) !$omp declare simd(operator_cross_qp) implicit none real(QP), dimension(:), intent(in) :: A, B - real(QP), dimension(3) :: C + real(QP), dimension(NDIM) :: C end function operator_cross_qp pure module function operator_cross_i1b(A, B) result(C) !$omp declare simd(operator_cross_i1b) implicit none integer(I1B), dimension(:), intent(in) :: A, B - integer(I1B), dimension(3) :: C + integer(I1B), dimension(NDIM) :: C end function operator_cross_i1b pure module function operator_cross_i2b(A, B) result(C) !$omp declare simd(operator_cross_i2b) implicit none integer(I2B), dimension(:), intent(in) :: A, B - integer(I2B), dimension(3) :: C + integer(I2B), dimension(NDIM) :: C end function operator_cross_i2b pure module function operator_cross_i4b(A, B) result(C) !$omp declare simd(operator_cross_i4b) implicit none integer(I4B), dimension(:), intent(in) :: A, B - integer(I4B), dimension(3) :: C + integer(I4B), dimension(NDIM) :: C end function operator_cross_i4b pure module function operator_cross_i8b(A, B) result(C) !$omp declare simd(operator_cross_i8b) implicit none integer(I8B), dimension(:), intent(in) :: A, B - integer(I8B), dimension(3) :: C + integer(I8B), dimension(NDIM) :: C end function operator_cross_i8b pure module function operator_cross_el_sp(A, B) result(C) @@ -116,7 +116,7 @@ end function operator_cross_el_i8b end interface !******************************************************************************************************************************** - ! Interfaces for .mag. operator + ! Interfaces for .mag. operator: Computes the magnitude of a vector or array of vectors using norm2 !******************************************************************************************************************************** interface operator(.mag.) @@ -160,4 +160,51 @@ pure module function operator_mag_el_qp(A) result(B) end function operator_mag_el_qp end interface + + !******************************************************************************************************************************** + ! Interfaces for .unit. operator: Returns a unit vector or array of unit vectors from an input vector or array of vectors + !******************************************************************************************************************************** + + interface operator(.unit.) + pure module function operator_unit_sp(A) result(B) + !$omp declare simd(operator_unit_sp) + implicit none + real(SP), dimension(:), intent(in) :: A + real(SP), dimension(NDIM) :: B + end function operator_unit_sp + + pure module function operator_unit_dp(A) result(B) + !$omp declare simd(operator_unit_dp) + implicit none + real(DP), dimension(:), intent(in) :: A + real(DP), dimension(NDIM) :: B + end function operator_unit_dp + + pure module function operator_unit_qp(A) result(B) + !$omp declare simd(operator_unit_qp) + implicit none + real(QP), dimension(:), intent(in) :: A + real(QP), dimension(NDIM) :: B + end function operator_unit_qp + + pure module function operator_unit_el_sp(A) result(B) + implicit none + real(SP), dimension(:,:), intent(in) :: A + real(SP), dimension(:,:), allocatable :: B + end function operator_unit_el_sp + + pure module function operator_unit_el_dp(A) result(B) + implicit none + real(DP), dimension(:,:), intent(in) :: A + real(DP), dimension(:,:), allocatable :: B + end function operator_unit_el_dp + + pure module function operator_unit_el_qp(A) result(B) + implicit none + real(QP), dimension(:,:), intent(in) :: A + real(QP), dimension(:,:), allocatable :: B + end function operator_unit_el_qp + end interface + + end module swiftest_operators diff --git a/src/operators/operator_cross.f90 b/src/operators/operator_cross.f90 index ba9582828..2a9af1ecf 100644 --- a/src/operators/operator_cross.f90 +++ b/src/operators/operator_cross.f90 @@ -12,6 +12,7 @@ !! author: David A. Minton !! !! Contains implementations for the .cross. operator for all defined integer and real types + !! Computes the cross product of two (3) vectors or (3,:) arrays !! Single vector implementations: C(1:3) = A(1:3) .cross. B(1:3) !! Vector list implementations: C(1:3, :) = A(1:3, :) .cross. B(1:3, :) contains diff --git a/src/operators/operator_mag.f90 b/src/operators/operator_mag.f90 index bea89d55b..2cf9e643c 100644 --- a/src/operators/operator_mag.f90 +++ b/src/operators/operator_mag.f90 @@ -11,6 +11,7 @@ !! author: David A. Minton !! !! Contains implementations for the .mag. operator for all defined real types + !! Computes the magnitude of a vector or array of vectors using norm2 !! Single vector implementations: B = .mag. A(1:3) !! Vector list implementations: B(:) = .mag. A(1:3, :) contains diff --git a/src/operators/operator_unit.f90 b/src/operators/operator_unit.f90 new file mode 100644 index 000000000..e9c68f28f --- /dev/null +++ b/src/operators/operator_unit.f90 @@ -0,0 +1,129 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version NDIM of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +submodule(swiftest_operators) s_operator_unit + !! author: David A. Minton + !! + !! Contains implementations for the .unit. operator for all defined real types + !! Returns a unit vector or array of unit vectors from an input vector or array of vectors + !! Single vector implementations: B = .unit. A(1:NDIM) + !! Vector list implementations: B(:) = .unit. A(1:NDIM, :) + contains + + pure module function operator_unit_sp(A) result(B) + implicit none + ! Arguments + real(SP), dimension(:), intent(in) :: A + real(SP), dimension(NDIM) :: B + ! Internals + real(SP) :: Amag + + Amag = norm2(A(:)) + B(:) = A(:) / Amag + + return + end function operator_unit_sp + + + pure module function operator_unit_dp(A) result(B) + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: A + real(DP), dimension(NDIM) :: B + ! Internals + real(DP) :: Amag + + Amag = norm2(A(:)) + B(:) = A(:) / Amag + + return + end function operator_unit_dp + + + pure module function operator_unit_qp(A) result(B) + implicit none + ! Arguments + real(QP), dimension(:), intent(in) :: A + real(QP), dimension(NDIM) :: B + ! Internals + real(QP) :: Amag + + Amag = norm2(A(:)) + B(:) = A(:) / Amag + + return + end function operator_unit_qp + + + pure module function operator_unit_el_sp(A) result(B) + implicit none + ! Arguments + real(SP), dimension(:,:), intent(in) :: A + real(SP), dimension(:,:), allocatable :: B + ! Internals + real(SP) :: Amag + integer(I4B) :: i,n + + n = size(A, 2) + if (allocated(B)) deallocate(B) + allocate(B(NDIM,n)) + + do concurrent (i=1:n) + Amag = norm2(A(:, i)) + B(:,i) = A(:,i) / Amag + end do + + return + end function operator_unit_el_sp + + + pure module function operator_unit_el_dp(A) result(B) + implicit none + ! Arguments + real(DP), dimension(:,:), intent(in) :: A + real(DP), dimension(:,:), allocatable :: B + ! Internals + real(DP) :: Amag + integer(I4B) :: i,n + + n = size(A, 2) + if (allocated(B)) deallocate(B) + allocate(B(NDIM,n)) + + do concurrent (i=1:n) + Amag = norm2(A(:, i)) + B(:,i) = A(:,i) / Amag + end do + + return + end function operator_unit_el_dp + + pure module function operator_unit_el_qp(A) result(B) + implicit none + ! Arguments + real(QP), dimension(:,:), intent(in) :: A + real(QP), dimension(:,:), allocatable :: B + ! Internals + real(QP) :: Amag + integer(I4B) :: i,n + + n = size(A, 2) + if (allocated(B)) deallocate(B) + allocate(B(NDIM,n)) + + do concurrent (i=1:n) + Amag = norm2(A(:, i)) + B(:,i) = A(:,i) / Amag + end do + + return + end function operator_unit_el_qp + +end submodule s_operator_unit +