diff --git a/src/simulation/m_qbmm.fpp b/src/simulation/m_qbmm.fpp index 1c52cee2fa..c952a7745d 100644 --- a/src/simulation/m_qbmm.fpp +++ b/src/simulation/m_qbmm.fpp @@ -421,258 +421,143 @@ contains real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: mv, rhs_mv integer :: i, j, k, l, q - real(wp) :: nb_q, nb_dot, R, R2, nR, nR2, nR_dot, nR2_dot, var, AX + logical :: is_axisym - if (idir == 1) then - - !Non-polytropic qbmm needs to account for change in bubble radius due to a change in nb - if (.not. polytropic) then - !$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var, AX) - do i = 1, nb - do q = 1, nnode - do l = 0, p - do k = 0, n - do j = 0, m - nb_q = q_cons_vf(bubxb + (i - 1)*nmom)%sf(j, k, l) - nR = q_cons_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) - nR2 = q_cons_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) - - R = q_prim_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) - R2 = q_prim_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) - - if (R2 - R**2._wp > 0._wp) then - var = R2 - R**2._wp - else - var = verysmall - end if - - if (q <= 2) then - AX = R - sqrt(var) - else - AX = R + sqrt(var) - end if + select case (idir) + case (1) + is_axisym = .false. + case (2) + is_axisym = .false. + case (3) + is_axisym = (grid_geometry == 3) + end select + if (.not. polytropic) then + !$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var, AX) + do i = 1, nb + do q = 1, nnode + do l = 0, p + do k = 0, n + do j = 0, m + nb_q = q_cons_vf(bubxb + (i - 1)*nmom)%sf(j, k, l) + nR = q_cons_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) + nR2 = q_cons_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) + R = q_prim_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) + R2 = q_prim_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) + var = max(R2 - R**2._wp, verysmall) + if (q <= 2) then + AX = R - sqrt(var) + else + AX = R + sqrt(var) + end if + select case (idir) + case (1) nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j - 1, k, l) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l) nR_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j - 1, k, l) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) nR2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j - 1, k, l) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dx(j)*AX*nb_q**2)* & (nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i)) - - if (q <= 2) then - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dx(j)*AX*nb_q**2*sqrt(var)*2._wp)* & - (nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i)) - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dx(j)*AX*nb_q**2*sqrt(var)*2._wp)* & - (-2._wp*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i)) - - else - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dx(j)*AX*nb_q**2*sqrt(var)*2._wp)* & - (nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i)) - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dx(j)*AX*nb_q**2*sqrt(var)*2._wp)* & - (-2._wp*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i)) - end if - - end do - end do - end do - end do - end do - end if - - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do q = 0, n - do i = 0, m - - rhs_vf(alf_idx)%sf(i, q, l) = rhs_vf(alf_idx)%sf(i, q, l) + mom_sp(2)%sf(i, q, l) - j = bubxb - !$acc loop seq - do k = 1, nb - rhs_vf(j)%sf(i, q, l) = & - rhs_vf(j)%sf(i, q, l) + mom_3d(0, 0, k)%sf(i, q, l) - rhs_vf(j + 1)%sf(i, q, l) = & - rhs_vf(j + 1)%sf(i, q, l) + mom_3d(1, 0, k)%sf(i, q, l) - rhs_vf(j + 2)%sf(i, q, l) = & - rhs_vf(j + 2)%sf(i, q, l) + mom_3d(0, 1, k)%sf(i, q, l) - rhs_vf(j + 3)%sf(i, q, l) = & - rhs_vf(j + 3)%sf(i, q, l) + mom_3d(2, 0, k)%sf(i, q, l) - rhs_vf(j + 4)%sf(i, q, l) = & - rhs_vf(j + 4)%sf(i, q, l) + mom_3d(1, 1, k)%sf(i, q, l) - rhs_vf(j + 5)%sf(i, q, l) = & - rhs_vf(j + 5)%sf(i, q, l) + mom_3d(0, 2, k)%sf(i, q, l) - j = j + 6 - end do - - end do - end do - end do - - elseif (idir == 2) then - - !Non-polytropic qbmm needs to account for change in bubble radius due to a change in nb - if (.not. polytropic) then - !$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var, AX) - do i = 1, nb - do q = 1, nnode - do l = 0, p - do k = 0, n - do j = 0, m - nb_q = q_cons_vf(bubxb + (i - 1)*nmom)%sf(j, k, l) - nR = q_cons_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) - nR2 = q_cons_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) - - R = q_prim_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) - R2 = q_prim_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) - - if (R2 - R**2._wp > 0._wp) then - var = R2 - R**2._wp - else - var = verysmall - end if - - if (q <= 2) then - AX = R - sqrt(var) - else - AX = R + sqrt(var) - end if - + case (2) nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k - 1, l) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l) nR_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k - 1, l) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) nR2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k - 1, l) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dy(k)*AX*nb_q**2)* & (nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i)) - - if (q <= 2) then - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dy(k)*AX*nb_q**2*sqrt(var)*2._wp)* & - (nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i)) - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dy(k)*AX*nb_q**2*sqrt(var)*2._wp)* & - (-2._wp*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i)) - - else - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dy(k)*AX*nb_q**2*sqrt(var)*2._wp)* & - (nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i)) - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dy(k)*AX*nb_q**2*sqrt(var)*2._wp)* & - (-2._wp*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i)) - end if - - end do - end do - end do - end do - end do - end if - - elseif (idir == 3) then - - if (.not. polytropic) then - if (grid_geometry == 3) then - !Non-polytropic qbmm needs to account for change in bubble radius due to a change in nb - !$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var, AX) - do i = 1, nb - do q = 1, nnode - do l = 0, p - do k = 0, n - do j = 0, m - nb_q = q_cons_vf(bubxb + (i - 1)*nmom)%sf(j, k, l) - nR = q_cons_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) - nR2 = q_cons_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) - - R = q_prim_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) - R2 = q_prim_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) - - if (R2 - R**2._wp > 0._wp) then - var = R2 - R**2._wp - else - var = verysmall - end if - - if (q <= 2) then - AX = R - sqrt(var) - else - AX = R + sqrt(var) - end if - + case (3) + if (is_axisym) then nb_dot = q_prim_vf(contxe + idir)%sf(j, k, l)*(flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l)) nR_dot = q_prim_vf(contxe + idir)%sf(j, k, l)*(flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)) nR2_dot = q_prim_vf(contxe + idir)%sf(j, k, l)*(flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)) - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dz(l)*y_cc(k)*AX*nb_q**2)* & (nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i)) - - if (q <= 2) then + else + nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l) + nR_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) + nR2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dz(l)*AX*nb_q**2)* & + (nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i)) + end if + end select + if (q <= 2) then + select case (idir) + case (1) + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dx(j)*AX*nb_q**2*sqrt(var)*2._wp)* & + (nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i)) + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dx(j)*AX*nb_q**2*sqrt(var)*2._wp)* & + (-2._wp*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i)) + case (2) + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dy(k)*AX*nb_q**2*sqrt(var)*2._wp)* & + (nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i)) + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dy(k)*AX*nb_q**2*sqrt(var)*2._wp)* & + (-2._wp*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i)) + case (3) + if (is_axisym) then rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dz(l)*y_cc(k)*AX*nb_q**2*sqrt(var)*2._wp)* & (nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i)) rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dz(l)*y_cc(k)*AX*nb_q**2*sqrt(var)*2._wp)* & (-2._wp*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i)) - else - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dz(l)*y_cc(k)*AX*nb_q**2*sqrt(var)*2._wp)* & + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dz(l)*AX*nb_q**2*sqrt(var)*2._wp)* & (nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i)) - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dz(l)*y_cc(k)*AX*nb_q**2*sqrt(var)*2._wp)* & + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dz(l)*AX*nb_q**2*sqrt(var)*2._wp)* & (-2._wp*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i)) end if - end do - end do - end do - end do - end do - else - !Non-polytropic qbmm needs to account for change in bubble radius due to a change in nb - !$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var, AX) - do i = 1, nb - do q = 1, nnode - do l = 0, p - do k = 0, n - do j = 0, m - nb_q = q_cons_vf(bubxb + (i - 1)*nmom)%sf(j, k, l) - nR = q_cons_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) - nR2 = q_cons_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) - - R = q_prim_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) - R2 = q_prim_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) - - if (R2 - R**2._wp > 0._wp) then - var = R2 - R**2._wp - else - var = verysmall - end if - - if (q <= 2) then - AX = R - sqrt(var) - else - AX = R + sqrt(var) - end if - - nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l) - nR_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l) - nR2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l - 1) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l) - - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dz(l)*AX*nb_q**2)* & - (nR_dot*nb_q - nR*nb_dot)*(pb(j, k, l, q, i)) - - if (q <= 2) then - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dz(l)*AX*nb_q**2*sqrt(var)*2._wp)* & + end select + else + select case (idir) + case (1) + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dx(j)*AX*nb_q**2*sqrt(var)*2._wp)* & + (nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i)) + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dx(j)*AX*nb_q**2*sqrt(var)*2._wp)* & + (-2._wp*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i)) + case (2) + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dy(k)*AX*nb_q**2*sqrt(var)*2._wp)* & + (nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i)) + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dy(k)*AX*nb_q**2*sqrt(var)*2._wp)* & + (-2._wp*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i)) + case (3) + if (is_axisym) then + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dz(l)*y_cc(k)*AX*nb_q**2*sqrt(var)*2._wp)* & (nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i)) - rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3._wp*gam/(dz(l)*AX*nb_q**2*sqrt(var)*2._wp)* & + rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dz(l)*y_cc(k)*AX*nb_q**2*sqrt(var)*2._wp)* & (-2._wp*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i)) - else rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dz(l)*AX*nb_q**2*sqrt(var)*2._wp)* & (nR2_dot*nb_q - nR2*nb_dot)*(pb(j, k, l, q, i)) rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3._wp*gam/(dz(l)*AX*nb_q**2*sqrt(var)*2._wp)* & (-2._wp*(nR/nb_q)*(nR_dot*nb_q - nR*nb_dot))*(pb(j, k, l, q, i)) end if - - end do - end do + end select + end if end do end do end do - end if - end if + end do + end do + end if + ! The following block is not repeated and is left as is + if (idir == 1) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do q = 0, n + do i = 0, m + rhs_vf(alf_idx)%sf(i, q, l) = rhs_vf(alf_idx)%sf(i, q, l) + mom_sp(2)%sf(i, q, l) + j = bubxb + !$acc loop seq + do k = 1, nb + rhs_vf(j)%sf(i, q, l) = rhs_vf(j)%sf(i, q, l) + mom_3d(0, 0, k)%sf(i, q, l) + rhs_vf(j + 1)%sf(i, q, l) = rhs_vf(j + 1)%sf(i, q, l) + mom_3d(1, 0, k)%sf(i, q, l) + rhs_vf(j + 2)%sf(i, q, l) = rhs_vf(j + 2)%sf(i, q, l) + mom_3d(0, 1, k)%sf(i, q, l) + rhs_vf(j + 3)%sf(i, q, l) = rhs_vf(j + 3)%sf(i, q, l) + mom_3d(2, 0, k)%sf(i, q, l) + rhs_vf(j + 4)%sf(i, q, l) = rhs_vf(j + 4)%sf(i, q, l) + mom_3d(1, 1, k)%sf(i, q, l) + rhs_vf(j + 5)%sf(i, q, l) = rhs_vf(j + 5)%sf(i, q, l) + mom_3d(0, 2, k)%sf(i, q, l) + j = j + 6 + end do + end do + end do + end do end if end subroutine s_compute_qbmm_rhs @@ -759,7 +644,7 @@ contains !$acc routine seq #endif - real(wp), intent(inout) :: pres, rho, c + real(wp), intent(in) :: pres, rho, c real(wp), dimension(nterms, 0:2, 0:2), intent(out) :: coeffs integer :: i1, i2 @@ -830,13 +715,9 @@ contains real(wp), dimension(nterms, 0:2, 0:2) :: coeff real(wp) :: pres, rho, nbub, c, alf, momsum, drdt, drdt2, chi_vw, x_vw, rho_mw, k_mw, T_bar, grad_T real(wp) :: n_tait, B_tait - - integer :: id1, id2, id3 - integer :: i1, i2 - integer :: j, q, r + integer :: id1, id2, id3, i1, i2, j, q, r is1_qbmm = ix; is2_qbmm = iy; is3_qbmm = iz - !$acc update device(is1_qbmm, is2_qbmm, is3_qbmm) !$acc parallel loop collapse(3) gang vector default(present) private(moms, msum, wght, abscX, abscY, wght_pb, wght_mv, wght_ht, coeff, ht, r, q, n_tait, B_tait, pres, rho, nbub, c, alf, momsum, drdt, drdt2, chi_vw, x_vw, rho_mw, k_mw, T_bar, grad_T) @@ -847,73 +728,51 @@ contains alf = q_prim_vf(alf_idx)%sf(id1, id2, id3) pres = q_prim_vf(E_idx)%sf(id1, id2, id3) rho = q_prim_vf(contxb)%sf(id1, id2, id3) + if (bubble_model == 2) then - n_tait = gammas(1) - n_tait = 1._wp/n_tait + 1._wp !make this the usual little 'gamma' + n_tait = 1._wp/gammas(1) + 1._wp B_tait = pi_infs(1)*(n_tait - 1)/n_tait c = n_tait*(pres + B_tait)*(1._wp - alf)/(rho) - - if (c > 0._wp) then - c = sqrt(c) - else - c = sgm_eps - end if + c = merge(sqrt(c), sgm_eps, c > 0._wp) end if - if (polytropic) then - call s_coeff(pres, rho, c, coeff) - else - call s_coeff_nonpoly(pres, rho, c, coeff) - end if + call s_coeff_selector(pres, rho, c, coeff, polytropic) - ! SHB: Manually adjusted pressure here for no-coupling case - ! pres = 1._wp/0.3_wp if (alf > small_alf) then - nbub = q_cons_vf(bubxb)%sf(id1, id2, id3) - !$acc loop seq do q = 1, nb - !Initialize moment set for each R0 bin - !$acc loop seq + ! Gather moments for this bubble bin do r = 2, nmom moms(r) = q_prim_vf(bubmoms(q, r))%sf(id1, id2, id3) end do - moms(1) = 1._wp - call s_chyqmom(moms, wght(:, q), abscX(:, q), abscY(:, q)) if (polytropic) then - !Account for bubble pressure pb0 at each R0 bin !$acc loop seq do j = 1, nnode wght_pb(j, q) = wght(j, q)*(pb0(q) - pv) end do else - !Account for bubble pressure, mass transfer rate and heat transfer rate in wght_pb, wght_mv and wght_ht using Preston model !$acc loop seq do j = 1, nnode chi_vw = 1._wp/(1._wp + R_v/R_n*(pb(id1, id2, id3, j, q)/pv - 1._wp)) x_vw = M_n*chi_vw/(M_v + (M_n - M_v)*chi_vw) - k_mw = x_vw*k_v(q)/(x_vw + (1._wp - x_vw)*phi_vn) & - + (1._wp - x_vw)*k_n(q)/(x_vw*phi_nv + 1._wp - x_vw) + k_mw = x_vw*k_v(q)/(x_vw + (1._wp - x_vw)*phi_vn) + (1._wp - x_vw)*k_n(q)/(x_vw*phi_nv + 1._wp - x_vw) rho_mw = pv/(chi_vw*R_v*Tw) rhs_mv(id1, id2, id3, j, q) = -Re_trans_c(q)*((mv(id1, id2, id3, j, q)/(mv(id1, id2, id3, j, q) + mass_n0(q))) - chi_vw) rhs_mv(id1, id2, id3, j, q) = rho_mw*rhs_mv(id1, id2, id3, j, q)/Pe_c/(1._wp - chi_vw)/abscX(j, q) - - T_bar = Tw*(pb(id1, id2, id3, j, q)/pb0(q))*(abscX(j, q)/R0(q))**3 & - *(mass_n0(q) + mass_v0(q))/(mass_n0(q) + mv(id1, id2, id3, j, q)) + T_bar = Tw*(pb(id1, id2, id3, j, q)/pb0(q))*(abscX(j, q)/R0(q))**3*(mass_n0(q) + mass_v0(q))/(mass_n0(q) + mv(id1, id2, id3, j, q)) grad_T = -Re_trans_T(q)*(T_bar - Tw) ht(j, q) = pb0(q)*k_mw*grad_T/Pe_T(q)/abscX(j, q) - wght_pb(j, q) = wght(j, q)*(pb(id1, id2, id3, j, q)) wght_mv(j, q) = wght(j, q)*(rhs_mv(id1, id2, id3, j, q)) wght_ht(j, q) = wght(j, q)*ht(j, q) end do end if - !Compute change in moments due to bubble dynamics + ! Compute change in moments due to bubble dynamics r = 1 !$acc loop seq do i2 = 0, 2 @@ -923,67 +782,45 @@ contains momsum = 0._wp !$acc loop seq do j = 1, nterms - ! Account for term with pb in Rayleigh Plesset equation - if (bubble_model == 3 .and. j == 3) then - momsum = momsum + coeff(j, i1, i2)*(R0(q)**momrhs(3, i1, i2, j, q)) & - *f_quad2D(abscX(:, q), abscY(:, q), wght_pb(:, q), momrhs(:, i1, i2, j, q)) - ! Account for terms with pb in Keller-Miksis equation - else if (bubble_model == 2 .and. ((j >= 7 .and. j <= 9) .or. (j >= 22 .and. j <= 23) .or. (j >= 10 .and. j <= 11) .or. (j == 26))) then - momsum = momsum + coeff(j, i1, i2)*(R0(q)**momrhs(3, i1, i2, j, q)) & - *f_quad2D(abscX(:, q), abscY(:, q), wght_pb(:, q), momrhs(:, i1, i2, j, q)) - ! Account for terms with mass transfer rate in Keller-Miksis equation - else if (bubble_model == 2 .and. (j >= 27 .and. j <= 29) .and. (.not. polytropic)) then - momsum = momsum + coeff(j, i1, i2)*(R0(q)**momrhs(3, i1, i2, j, q)) & - *f_quad2D(abscX(:, q), abscY(:, q), wght_mv(:, q), momrhs(:, i1, i2, j, q)) - ! Account for terms with heat transfer rate in Keller-Miksis equation - else if (bubble_model == 2 .and. (j >= 30 .and. j <= 32) .and. (.not. polytropic)) then - momsum = momsum + coeff(j, i1, i2)*(R0(q)**momrhs(3, i1, i2, j, q)) & - *f_quad2D(abscX(:, q), abscY(:, q), wght_ht(:, q), momrhs(:, i1, i2, j, q)) - else - momsum = momsum + coeff(j, i1, i2)*(R0(q)**momrhs(3, i1, i2, j, q)) & - *f_quad2D(abscX(:, q), abscY(:, q), wght(:, q), momrhs(:, i1, i2, j, q)) - end if - + select case (bubble_model) + case (3) + if (j == 3) then + momsum = momsum + coeff(j, i1, i2)*(R0(q)**momrhs(3, i1, i2, j, q))*f_quad2D(abscX(:, q), abscY(:, q), wght_pb(:, q), momrhs(:, i1, i2, j, q)) + else + momsum = momsum + coeff(j, i1, i2)*(R0(q)**momrhs(3, i1, i2, j, q))*f_quad2D(abscX(:, q), abscY(:, q), wght(:, q), momrhs(:, i1, i2, j, q)) + end if + case (2) + if ((j >= 7 .and. j <= 9) .or. (j >= 22 .and. j <= 23) .or. (j >= 10 .and. j <= 11) .or. (j == 26)) then + momsum = momsum + coeff(j, i1, i2)*(R0(q)**momrhs(3, i1, i2, j, q))*f_quad2D(abscX(:, q), abscY(:, q), wght_pb(:, q), momrhs(:, i1, i2, j, q)) + else if ((j >= 27 .and. j <= 29) .and. (.not. polytropic)) then + momsum = momsum + coeff(j, i1, i2)*(R0(q)**momrhs(3, i1, i2, j, q))*f_quad2D(abscX(:, q), abscY(:, q), wght_mv(:, q), momrhs(:, i1, i2, j, q)) + else if ((j >= 30 .and. j <= 32) .and. (.not. polytropic)) then + momsum = momsum + coeff(j, i1, i2)*(R0(q)**momrhs(3, i1, i2, j, q))*f_quad2D(abscX(:, q), abscY(:, q), wght_ht(:, q), momrhs(:, i1, i2, j, q)) + else + momsum = momsum + coeff(j, i1, i2)*(R0(q)**momrhs(3, i1, i2, j, q))*f_quad2D(abscX(:, q), abscY(:, q), wght(:, q), momrhs(:, i1, i2, j, q)) + end if + end select end do - moms3d(i1, i2, q)%sf(id1, id2, id3) = nbub*momsum msum(r) = momsum r = r + 1 - end if end do end do - ! Compute change in pb and mv for non-polytroic model + ! Compute change in pb and mv for non-polytropic model if (.not. polytropic) then !$acc loop seq do j = 1, nnode - ! Compute Rdot (drdt) at quadrature node in the ODE for pb (note this is not the same as bubble variable Rdot) drdt = msum(2) - if (moms(4) - moms(2)**2._wp > 0._wp) then - if (j == 1 .or. j == 2) then - drdt2 = -1._wp/(2._wp*sqrt(moms(4) - moms(2)**2._wp)) - else - drdt2 = 1._wp/(2._wp*sqrt(moms(4) - moms(2)**2._wp)) - end if - else - ! Edge case where variance < 0 - if (j == 1 .or. j == 2) then - drdt2 = -1._wp/(2._wp*sqrt(verysmall)) - else - drdt2 = 1._wp/(2._wp*sqrt(verysmall)) - end if - end if - + drdt2 = merge(-1._wp, 1._wp, j == 1 .or. j == 2)/(2._wp*sqrt(merge(moms(4) - moms(2)**2._wp, verysmall, moms(4) - moms(2)**2._wp > 0._wp))) drdt2 = drdt2*(msum(3) - 2._wp*moms(2)*msum(2)) drdt = drdt + drdt2 - rhs_pb(id1, id2, id3, j, q) = (-3._wp*gam*drdt/abscX(j, q))*(pb(id1, id2, id3, j, q)) rhs_pb(id1, id2, id3, j, q) = rhs_pb(id1, id2, id3, j, q) + (3._wp*gam/abscX(j, q))*rhs_mv(id1, id2, id3, j, q)*R_v*Tw rhs_pb(id1, id2, id3, j, q) = rhs_pb(id1, id2, id3, j, q) + (3._wp*gam/abscX(j, q))*ht(j, q) rhs_mv(id1, id2, id3, j, q) = rhs_mv(id1, id2, id3, j, q)*(4._wp*pi*abscX(j, q)**2._wp) end do - end if end do @@ -992,19 +829,14 @@ contains momsp(2)%sf(id1, id2, id3) = 4._wp*pi*nbub*f_quad(abscX, abscY, wght, 2._wp, 1._wp, 0._wp) momsp(3)%sf(id1, id2, id3) = f_quad(abscX, abscY, wght, 3._wp, 2._wp, 0._wp) if (abs(gam - 1._wp) <= 1.e-4_wp) then - ! Gam \approx 1, don't risk imaginary quadrature momsp(4)%sf(id1, id2, id3) = 1._wp else - !Special moment with bubble pressure pb if (polytropic) then - momsp(4)%sf(id1, id2, id3) = f_quad(abscX, abscY, wght_pb, 3._wp*(1._wp - gam), 0._wp, 3._wp*gam) + pv*f_quad(abscX, abscY, wght, 3._wp, 0._wp, 0._wp) & - - 4._wp*Re_inv*f_quad(abscX, abscY, wght, 2._wp, 1._wp, 0._wp) - (2._wp/Web)*f_quad(abscX, abscY, wght, 2._wp, 0._wp, 0._wp) + momsp(4)%sf(id1, id2, id3) = f_quad(abscX, abscY, wght_pb, 3._wp*(1._wp - gam), 0._wp, 3._wp*gam) + pv*f_quad(abscX, abscY, wght, 3._wp, 0._wp, 0._wp) - 4._wp*Re_inv*f_quad(abscX, abscY, wght, 2._wp, 1._wp, 0._wp) - (2._wp/Web)*f_quad(abscX, abscY, wght, 2._wp, 0._wp, 0._wp) else - momsp(4)%sf(id1, id2, id3) = f_quad(abscX, abscY, wght_pb, 3._wp, 0._wp, 0._wp) & - - 4._wp*Re_inv*f_quad(abscX, abscY, wght, 2._wp, 1._wp, 0._wp) - (2._wp/Web)*f_quad(abscX, abscY, wght, 2._wp, 0._wp, 0._wp) + momsp(4)%sf(id1, id2, id3) = f_quad(abscX, abscY, wght_pb, 3._wp, 0._wp, 0._wp) - 4._wp*Re_inv*f_quad(abscX, abscY, wght, 2._wp, 1._wp, 0._wp) - (2._wp/Web)*f_quad(abscX, abscY, wght, 2._wp, 0._wp, 0._wp) end if end if - else !$acc loop seq do q = 1, nb @@ -1016,133 +848,141 @@ contains end do end do end do - momsp(1)%sf(id1, id2, id3) = 0._wp momsp(2)%sf(id1, id2, id3) = 0._wp momsp(3)%sf(id1, id2, id3) = 0._wp momsp(4)%sf(id1, id2, id3) = 0._wp - end if - end do end do end do - end subroutine s_mom_inv + contains + ! Helper to select the correct coefficient routine + subroutine s_coeff_selector(pres, rho, c, coeff, polytropic) +#ifdef _CRAYFTN + !DIR$ INLINEALWAYS s_chyqmom +#else + !$acc routine seq +#endif + real(wp), intent(in) :: pres, rho, c + real(wp), dimension(nterms, 0:2, 0:2), intent(out) :: coeff + 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 - pure subroutine s_chyqmom(momin, wght, abscX, abscY) + pure subroutine s_chyqmom(momin, wght, abscX, abscY) #ifdef _CRAYFTN - !DIR$ INLINEALWAYS s_chyqmom + !DIR$ INLINEALWAYS s_chyqmom #else - !$acc routine seq + !$acc routine seq #endif - real(wp), dimension(nmom), intent(in) :: momin - real(wp), dimension(nnode), intent(inout) :: wght, abscX, abscY - - 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) :: mu2avg, mu2, vp21, vp22, rho21, rho22 - - 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) - - 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; - M1 = (/1._wp, 0._wp, c20/) - call s_hyqmom(myrho, up, M1) - Vf = c11*up/c20 - - mu2avg = c02 - sum(myrho(:)*(Vf(:)**2._wp)) - - mu2avg = maxval((/mu2avg, 0._wp/)) - mu2 = mu2avg - M3 = (/1._wp, 0._wp, mu2/) - call s_hyqmom(myrho3, up3, M3) - - vp21 = up3(1) - vp22 = up3(2) - rho21 = myrho3(1) - rho22 = myrho3(2) - - wght(1) = myrho(1)*rho21 - wght(2) = myrho(1)*rho22 - wght(3) = myrho(2)*rho21 - wght(4) = myrho(2)*rho22 - wght = moms(0, 0)*wght - - abscX(1) = up(1) - abscX(2) = up(1) - abscX(3) = up(2) - abscX(4) = up(2) - abscX = bu + abscX - - abscY(1) = Vf(1) + vp21 - abscY(2) = Vf(1) + vp22 - abscY(3) = Vf(2) + vp21 - abscY(4) = Vf(2) + vp22 - abscY = bv + abscY - - end subroutine s_chyqmom - - pure subroutine s_hyqmom(frho, fup, fmom) + 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 + + pure subroutine s_hyqmom(frho, fup, fmom) #ifdef _CRAYFTN - !DIR$ INLINEALWAYS s_hyqmom + !DIR$ INLINEALWAYS s_hyqmom #else - !$acc routine seq + !$acc routine seq #endif - real(wp), dimension(2), intent(inout) :: frho, fup - real(wp), dimension(3), intent(in) :: fmom + real(wp), dimension(2), intent(inout) :: frho, fup + real(wp), dimension(3), intent(in) :: fmom - real(wp) :: bu, d2, c2 + 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, verysmall/)) - fup(1) = bu - sqrt(c2) - fup(2) = bu + sqrt(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, verysmall/)) + fup(1) = bu - sqrt(c2) + fup(2) = bu + sqrt(c2) - end subroutine s_hyqmom + end subroutine s_hyqmom - pure function f_quad(abscX, abscY, wght_in, q, r, s) - !$acc routine seq - real(wp), dimension(nnode, nb), intent(in) :: abscX, abscY, wght_in - real(wp), intent(in) :: q, r, s + pure function f_quad(abscX, abscY, wght_in, q, r, s) + !$acc routine seq + real(wp), dimension(nnode, nb), intent(in) :: abscX, abscY, wght_in + real(wp), intent(in) :: q, r, s - real(wp) :: f_quad_RV, f_quad - integer :: i + real(wp) :: f_quad_RV, f_quad + integer :: i - f_quad = 0._wp - do i = 1, nb - f_quad_RV = sum(wght_in(:, i)*(abscX(:, i)**q)*(abscY(:, i)**r)) - f_quad = f_quad + weight(i)*(R0(i)**s)*f_quad_RV - end do + f_quad = 0._wp + do i = 1, nb + f_quad_RV = sum(wght_in(:, i)*(abscX(:, i)**q)*(abscY(:, i)**r)) + f_quad = f_quad + weight(i)*(R0(i)**s)*f_quad_RV + end do - end function f_quad + end function f_quad - pure function f_quad2D(abscX, abscY, wght_in, pow) - !$acc routine seq - real(wp), dimension(nnode), intent(in) :: abscX, abscY, wght_in - real(wp), dimension(3), intent(in) :: pow + pure function f_quad2D(abscX, abscY, wght_in, pow) + !$acc routine seq + real(wp), dimension(nnode), intent(in) :: abscX, abscY, wght_in + real(wp), dimension(3), intent(in) :: pow + + real(wp) :: f_quad2D - real(wp) :: f_quad2D + f_quad2D = sum(wght_in(:)*(abscX(:)**pow(1))*(abscY(:)**pow(2))) + end function f_quad2D - f_quad2D = sum(wght_in(:)*(abscX(:)**pow(1))*(abscY(:)**pow(2))) - end function f_quad2D + end subroutine s_mom_inv end module m_qbmm