Skip to content

qbmm: hoist s_mom_inv's internal procedures to module scope #1688

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
sbryngelson wants to merge 1 commit into
base: master
Choose a base branch
from
+137 −138
Open

qbmm: hoist s_mom_inv's internal procedures to module scope #1688

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
275 changes: 137 additions & 138 deletions src/simulation/m_qbmm.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -925,152 +925,151 @@ contains
end do
$:END_GPU_PARALLEL_LOOP()

contains
!> Select the polytropic or non-polytropic coefficient routine
subroutine s_coeff_selector(pres, rho, c, coeff, polytropic)

$:GPU_ROUTINE(function_name='s_coeff_selector',parallelism='[seq]', cray_inline=True)
real(wp), intent(in) :: pres, rho, c
#:if USING_AMD
real(wp), dimension(32,0:2,0:2), intent(out) :: coeff
#:else
real(wp), dimension(nterms,0:2,0:2), intent(out) :: coeff
#:endif
logical, intent(in) :: polytropic
if (polytropic) then
call s_coeff(pres, rho, c, coeff)
else
call s_coeff_nonpoly(pres, rho, c, coeff)
end if
end subroutine s_mom_inv

end subroutine s_coeff_selector

!> Perform CHyQMOM inversion for bivariate moments
subroutine s_chyqmom(momin, wght, abscX, abscY)

$:GPU_ROUTINE(function_name='s_chyqmom',parallelism='[seq]', cray_inline=True)

real(wp), dimension(nmom), intent(in) :: momin
real(wp), dimension(nnode), intent(inout) :: wght, abscX, abscY

! Local variables
real(wp), dimension(0:2,0:2) :: moms
real(wp), dimension(3) :: M1, M3
real(wp), dimension(2) :: myrho, myrho3, up, up3, Vf
real(wp) :: bu, bv, d20, d11, d_02, c20, c11, c02
real(wp) :: mu2, vp21, vp22, rho21, rho22

! Assign moments to 2D array for clarity
moms(0, 0) = momin(1)
moms(1, 0) = momin(2)
moms(0, 1) = momin(3)
moms(2, 0) = momin(4)
moms(1, 1) = momin(5)
moms(0, 2) = momin(6)

! Compute means and central moments
bu = moms(1, 0)/moms(0, 0)
bv = moms(0, 1)/moms(0, 0)
d20 = moms(2, 0)/moms(0, 0)
d11 = moms(1, 1)/moms(0, 0)
d_02 = moms(0, 2)/moms(0, 0)

c20 = d20 - bu**2._wp
c11 = d11 - bu*bv
c02 = d_02 - bv**2._wp

! First 1D quadrature (X direction)
M1 = (/1._wp, 0._wp, c20/)
call s_hyqmom(myrho, up, M1)
Vf = c11*up/c20

! Second 1D quadrature (Y direction, conditional on X)
mu2 = max(0._wp, c02 - sum(myrho*(Vf**2._wp)))
M3 = (/1._wp, 0._wp, mu2/)
call s_hyqmom(myrho3, up3, M3)

! Assign roots and weights for 2D quadrature
vp21 = up3(1)
vp22 = up3(2)
rho21 = myrho3(1)
rho22 = myrho3(2)

! Compute weights (vectorized)
wght = moms(0, 0)*[myrho(1)*rho21, myrho(1)*rho22, myrho(2)*rho21, myrho(2)*rho22]

! Compute abscissas (vectorized)
abscX = bu + [up(1), up(1), up(2), up(2)]
abscY = bv + [Vf(1) + vp21, Vf(1) + vp22, Vf(2) + vp21, Vf(2) + vp22]

end subroutine s_chyqmom

!> Perform HyQMOM inversion for univariate moments
subroutine s_hyqmom(frho, fup, fmom)

$:GPU_ROUTINE(function_name='s_hyqmom',parallelism='[seq]', cray_inline=True)

real(wp), dimension(2), intent(inout) :: frho, fup
real(wp), dimension(3), intent(in) :: fmom
real(wp) :: bu, d2, c2

bu = fmom(2)/fmom(1)
d2 = fmom(3)/fmom(1)
c2 = d2 - bu**2._wp
frho(1) = fmom(1)/2._wp
frho(2) = fmom(1)/2._wp
c2 = maxval((/c2, sgm_eps/))
fup(1) = bu - sqrt(c2)
fup(2) = bu + sqrt(c2)

end subroutine s_hyqmom

!> Evaluate a weighted quadrature sum over all bubble size bins and nodes
function f_quad(abscX, abscY, wght_in, q, r, s)

$:GPU_ROUTINE(parallelism='[seq]')
#:if not MFC_CASE_OPTIMIZATION and USING_AMD
real(wp), dimension(4, 3), intent(in) :: abscX, abscY, wght_in
#:else
real(wp), dimension(nnode, nb), intent(in) :: abscX, abscY, wght_in
#:endif
real(wp), intent(in) :: q, r, s
real(wp) :: f_quad_RV, f_quad
integer :: i, i1

f_quad = 0._wp
$:GPU_LOOP(parallelism='[seq]')
do i = 1, nb
f_quad_RV = 0._wp
$:GPU_LOOP(parallelism='[seq]')
do i1 = 1, nnode
f_quad_RV = f_quad_RV + wght_in(i1, i)*(abscX(i1, i)**q)*(abscY(i1, i)**r)
end do
f_quad = f_quad + weight(i)*(R0(i)**s)*f_quad_RV
end do
!> Select the polytropic or non-polytropic coefficient routine
subroutine s_coeff_selector(pres, rho, c, coeff, polytropic)

$:GPU_ROUTINE(function_name='s_coeff_selector',parallelism='[seq]', cray_inline=True)
real(wp), intent(in) :: pres, rho, c
#:if USING_AMD
real(wp), dimension(32,0:2,0:2), intent(out) :: coeff
#:else
real(wp), dimension(nterms,0:2,0:2), intent(out) :: coeff
#:endif
logical, intent(in) :: polytropic
if (polytropic) then
call s_coeff(pres, rho, c, coeff)
else
call s_coeff_nonpoly(pres, rho, c, coeff)
end if

end subroutine s_coeff_selector

!> Perform CHyQMOM inversion for bivariate moments
subroutine s_chyqmom(momin, wght, abscX, abscY)

$:GPU_ROUTINE(function_name='s_chyqmom',parallelism='[seq]', cray_inline=True)

real(wp), dimension(nmom), intent(in) :: momin
real(wp), dimension(nnode), intent(inout) :: wght, abscX, abscY

! Local variables
real(wp), dimension(0:2,0:2) :: moms
real(wp), dimension(3) :: M1, M3
real(wp), dimension(2) :: myrho, myrho3, up, up3, Vf
real(wp) :: bu, bv, d20, d11, d_02, c20, c11, c02
real(wp) :: mu2, vp21, vp22, rho21, rho22

! Assign moments to 2D array for clarity
moms(0, 0) = momin(1)
moms(1, 0) = momin(2)
moms(0, 1) = momin(3)
moms(2, 0) = momin(4)
moms(1, 1) = momin(5)
moms(0, 2) = momin(6)

! Compute means and central moments
bu = moms(1, 0)/moms(0, 0)
bv = moms(0, 1)/moms(0, 0)
d20 = moms(2, 0)/moms(0, 0)
d11 = moms(1, 1)/moms(0, 0)
d_02 = moms(0, 2)/moms(0, 0)

c20 = d20 - bu**2._wp
c11 = d11 - bu*bv
c02 = d_02 - bv**2._wp

end function f_quad
! First 1D quadrature (X direction)
M1 = (/1._wp, 0._wp, c20/)
call s_hyqmom(myrho, up, M1)
Vf = c11*up/c20

!> Evaluate a weighted 2D quadrature sum over quadrature nodes for a single size bin
function f_quad2D(abscX, abscY, wght_in, pow)
! Second 1D quadrature (Y direction, conditional on X)
mu2 = max(0._wp, c02 - sum(myrho*(Vf**2._wp)))
M3 = (/1._wp, 0._wp, mu2/)
call s_hyqmom(myrho3, up3, M3)

$:GPU_ROUTINE(parallelism='[seq]')
#:if not MFC_CASE_OPTIMIZATION and USING_AMD
real(wp), dimension(4), intent(in) :: abscX, abscY, wght_in
#:else
real(wp), dimension(nnode), intent(in) :: abscX, abscY, wght_in
#:endif
real(wp), dimension(3), intent(in) :: pow
real(wp) :: f_quad2D
integer :: i
! Assign roots and weights for 2D quadrature
vp21 = up3(1)
vp22 = up3(2)
rho21 = myrho3(1)
rho22 = myrho3(2)

f_quad2D = 0._wp
! Compute weights (vectorized)
wght = moms(0, 0)*[myrho(1)*rho21, myrho(1)*rho22, myrho(2)*rho21, myrho(2)*rho22]

! Compute abscissas (vectorized)
abscX = bu + [up(1), up(1), up(2), up(2)]
abscY = bv + [Vf(1) + vp21, Vf(1) + vp22, Vf(2) + vp21, Vf(2) + vp22]

end subroutine s_chyqmom

!> Perform HyQMOM inversion for univariate moments
subroutine s_hyqmom(frho, fup, fmom)

$:GPU_ROUTINE(function_name='s_hyqmom',parallelism='[seq]', cray_inline=True)

real(wp), dimension(2), intent(inout) :: frho, fup
real(wp), dimension(3), intent(in) :: fmom
real(wp) :: bu, d2, c2

bu = fmom(2)/fmom(1)
d2 = fmom(3)/fmom(1)
c2 = d2 - bu**2._wp
frho(1) = fmom(1)/2._wp
frho(2) = fmom(1)/2._wp
c2 = maxval((/c2, sgm_eps/))
fup(1) = bu - sqrt(c2)
fup(2) = bu + sqrt(c2)

end subroutine s_hyqmom

!> Evaluate a weighted quadrature sum over all bubble size bins and nodes
function f_quad(abscX, abscY, wght_in, q, r, s)

$:GPU_ROUTINE(parallelism='[seq]')
#:if not MFC_CASE_OPTIMIZATION and USING_AMD
real(wp), dimension(4, 3), intent(in) :: abscX, abscY, wght_in
#:else
real(wp), dimension(nnode, nb), intent(in) :: abscX, abscY, wght_in
#:endif
real(wp), intent(in) :: q, r, s
real(wp) :: f_quad_RV, f_quad
integer :: i, i1

f_quad = 0._wp
$:GPU_LOOP(parallelism='[seq]')
do i = 1, nb
f_quad_RV = 0._wp
$:GPU_LOOP(parallelism='[seq]')
do i = 1, nnode
f_quad2D = f_quad2D + wght_in(i)*(abscX(i)**pow(1))*(abscY(i)**pow(2))
do i1 = 1, nnode
f_quad_RV = f_quad_RV + wght_in(i1, i)*(abscX(i1, i)**q)*(abscY(i1, i)**r)
end do
f_quad = f_quad + weight(i)*(R0(i)**s)*f_quad_RV
end do

end function f_quad2D
end function f_quad

end subroutine s_mom_inv
!> Evaluate a weighted 2D quadrature sum over quadrature nodes for a single size bin
function f_quad2D(abscX, abscY, wght_in, pow)

$:GPU_ROUTINE(parallelism='[seq]')
#:if not MFC_CASE_OPTIMIZATION and USING_AMD
real(wp), dimension(4), intent(in) :: abscX, abscY, wght_in
#:else
real(wp), dimension(nnode), intent(in) :: abscX, abscY, wght_in
#:endif
real(wp), dimension(3), intent(in) :: pow
real(wp) :: f_quad2D
integer :: i

f_quad2D = 0._wp
$:GPU_LOOP(parallelism='[seq]')
do i = 1, nnode
f_quad2D = f_quad2D + wght_in(i)*(abscX(i)**pow(1))*(abscY(i)**pow(2))
end do

end function f_quad2D

end module m_qbmm
Loading