diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 55d05e823..700ecbe40 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -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 @@ -40,52 +40,48 @@ 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 @@ -93,7 +89,6 @@ module subroutine util_get_energy_momentum_system(self, param) 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) @@ -118,7 +113,6 @@ 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 @@ -126,14 +120,14 @@ module subroutine util_get_energy_momentum_system(self, param) 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