diff --git a/src/simulation/m_qbmm.fpp b/src/simulation/m_qbmm.fpp index abb77a1f35..0a85a0a156 100644 --- a/src/simulation/m_qbmm.fpp +++ b/src/simulation/m_qbmm.fpp @@ -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