Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Removed simd directives and replaced with do concurrent with block fo…
Browse files Browse the repository at this point in the history 
…r private variables
  • Loading branch information
@daminton
daminton committed Aug 10, 2021
1 parent 0d158a5 commit e48c71b
Showing 1 changed file with 41 additions and 47 deletions.
88 changes: 41 additions & 47 deletions src/util/util_get_energy_momentum.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,12 +15,12 @@ module subroutine util_get_energy_momentum_system(self, param)
! Internals
integer(I4B) :: i, j
integer(I8B) :: k
real(DP) :: rmag, v2, rot2, oblpot, hx, hy, hz, hsx, hsy, hsz
real(DP) :: kecb, kespincb, Lcborbitx, Lcborbity, Lcborbitz, Lcbspinx, Lcbspiny, Lcbspinz
real(DP) :: oblpot, kecb, kespincb
real(DP), dimension(self%pl%nbody) :: irh, kepl, kespinpl, pecb
real(DP), dimension(self%pl%nbody) :: Lplorbitx, Lplorbity, Lplorbitz
real(DP), dimension(self%pl%nbody) :: Lplspinx, Lplspiny, Lplspinz
real(DP), dimension(self%pl%nplpl) :: pepl
real(DP), dimension(NDIM) :: Lcborbit, Lcbspin
logical, dimension(self%pl%nplpl) :: lstatpl
logical, dimension(self%pl%nbody) :: lstatus

Expand All @@ -40,60 +40,55 @@ module subroutine util_get_energy_momentum_system(self, param)
lstatus(1:npl) = pl%status(1:npl) /= INACTIVE

kecb = cb%mass * dot_product(cb%vb(:), cb%vb(:))
hx = cb%xb(2) * cb%vb(3) - cb%xb(3) * cb%vb(2)
hy = cb%xb(3) * cb%vb(1) - cb%xb(1) * cb%vb(3)
hz = cb%xb(1) * cb%vb(2) - cb%xb(2) * cb%vb(1)
Lcborbitx = cb%mass * hx
Lcborbity = cb%mass * hy
Lcborbitz = cb%mass * hz
!!$omp simd private(v2, rot2, hx, hy, hz)
do i = 1, npl
v2 = dot_product(pl%vb(:,i), pl%vb(:,i))
hx = pl%xb(2,i) * pl%vb(3,i) - pl%xb(3,i) * pl%vb(2,i)
hy = pl%xb(3,i) * pl%vb(1,i) - pl%xb(1,i) * pl%vb(3,i)
hz = pl%xb(1,i) * pl%vb(2,i) - pl%xb(2,i) * pl%vb(1,i)
Lcborbit(:) = cb%mass * cb%xb(:) .cross. cb%vb(:)

! Angular momentum from orbit
Lplorbitx(i) = pl%mass(i) * hx
Lplorbity(i) = pl%mass(i) * hy
Lplorbitz(i) = pl%mass(i) * hz
do concurrent (i = 1:npl, lstatus(i))
block ! We use a block construct to prevent generating temporary arrays for local variables
real(DP) :: v2, hx, hy, hz
v2 = dot_product(pl%vb(:,i), pl%vb(:,i))
hx = pl%xb(2,i) * pl%vb(3,i) - pl%xb(3,i) * pl%vb(2,i)
hy = pl%xb(3,i) * pl%vb(1,i) - pl%xb(1,i) * pl%vb(3,i)
hz = pl%xb(1,i) * pl%vb(2,i) - pl%xb(2,i) * pl%vb(1,i)

! Angular momentum from orbit
Lplorbitx(i) = pl%mass(i) * hx
Lplorbity(i) = pl%mass(i) * hy
Lplorbitz(i) = pl%mass(i) * hz

! Kinetic energy from orbit and spin
kepl(i) = pl%mass(i) * v2
! Kinetic energy from orbit and spin
kepl(i) = pl%mass(i) * v2
end block
end do

if (param%lrotation) then
kespincb = cb%mass * cb%Ip(3) * cb%radius**2 * dot_product(cb%rot(:), cb%rot(:))
! For simplicity, we always assume that the rotation pole is the 3rd principal axis
hsx = cb%Ip(3) * cb%radius**2 * cb%rot(1)
hsy = cb%Ip(3) * cb%radius**2 * cb%rot(2)
hsz = cb%Ip(3) * cb%radius**2 * cb%rot(3)

! Angular momentum from spin
Lcbspinx = cb%mass * hsx
Lcbspiny = cb%mass * hsy
Lcbspinz = cb%mass * hsz
! For simplicity, we always assume that the rotation pole is the 3rd principal axis
Lcbspin(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:)

do i = 1, npl
rot2 = dot_product(pl%rot(:,i), pl%rot(:,i))
! For simplicity, we always assume that the rotation pole is the 3rd principal axis
hsx = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(1,i)
hsy = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(2,i)
hsz = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(3,i)
do concurrent (i = 1:npl, lstatus(i))
block
real(DP) :: rot2, hsx, hsy, hsz

! Angular momentum from spin
Lplspinx(i) = pl%mass(i) * hsx
Lplspiny(i) = pl%mass(i) * hsy
Lplspinz(i) = pl%mass(i) * hsz
kespinpl(i) = pl%mass(i) * pl%Ip(3, i) * pl%radius(i)**2 * rot2
rot2 = dot_product(pl%rot(:,i), pl%rot(:,i))
! For simplicity, we always assume that the rotation pole is the 3rd principal axis
hsx = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(1,i)
hsy = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(2,i)
hsz = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(3,i)

! Angular momentum from spin
Lplspinx(i) = pl%mass(i) * hsx
Lplspiny(i) = pl%mass(i) * hsy
Lplspinz(i) = pl%mass(i) * hsz
kespinpl(i) = pl%mass(i) * pl%Ip(3, i) * pl%radius(i)**2 * rot2
end block
end do
else
kespincb = 0.0_DP
kespinpl(:) = 0.0_DP
end if

! Do the central body potential energy component first
!$omp simd
associate(px => pl%xb(1,:), py => pl%xb(2,:), pz => pl%xb(3,:))
do concurrent(i = 1:npl, lstatus(i))
pecb(i) = -cb%Gmass * pl%mass(i) / sqrt(px(i)**2 + py(i)**2 + pz(i)**2)
Expand All @@ -118,22 +113,21 @@ module subroutine util_get_energy_momentum_system(self, param)

! Potential energy from the oblateness term
if (param%loblatecb) then
!$omp simd
do concurrent(i = 1:npl, lstatus(i))
irh(i) = 1.0_DP / norm2(pl%xh(:,i))
end do
call obl_pot(npl, cb%Gmass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot)
system%pe = system%pe + oblpot
end if

system%Lorbit(1) = Lcborbitx + sum(Lplorbitx(1:npl), lstatus(1:npl))
system%Lorbit(2) = Lcborbity + sum(Lplorbity(1:npl), lstatus(1:npl))
system%Lorbit(3) = Lcborbitz + sum(Lplorbitz(1:npl), lstatus(1:npl))
system%Lorbit(1) = Lcborbit(1) + sum(Lplorbitx(1:npl), lstatus(1:npl))
system%Lorbit(2) = Lcborbit(2) + sum(Lplorbity(1:npl), lstatus(1:npl))
system%Lorbit(3) = Lcborbit(3) + sum(Lplorbitz(1:npl), lstatus(1:npl))

if (param%lrotation) then
system%Lspin(1) = Lcbspinx + sum(Lplspinx(1:npl), lstatus(1:npl))
system%Lspin(2) = Lcbspiny + sum(Lplspiny(1:npl), lstatus(1:npl))
system%Lspin(3) = Lcbspinz + sum(Lplspinz(1:npl), lstatus(1:npl))
system%Lspin(1) = Lcbspin(1) + sum(Lplspinx(1:npl), lstatus(1:npl))
system%Lspin(2) = Lcbspin(2) + sum(Lplspiny(1:npl), lstatus(1:npl))
system%Lspin(3) = Lcbspin(3) + sum(Lplspinz(1:npl), lstatus(1:npl))
end if
end associate

Expand Down

0 comments on commit e48c71b

Please sign in to comment.