diff --git a/docs/documentation/equations.md b/docs/documentation/equations.md index b8fb3b956e..da6216aa0e 100644 --- a/docs/documentation/equations.md +++ b/docs/documentation/equations.md @@ -644,7 +644,7 @@ Reaction mechanisms are code-generated via Pyrometheus (\cite Cisneros26), which ## 10. Surface Tension (`surface_tension = .true.`) (\cite Schmidmayer17; \cite Wilfong26 Sec. 4.1.8) -**Source:** `src/simulation/m_surface_tension.fpp`, `src/simulation/include/inline_capillary.fpp` +**Source:** `src/simulation/m_surface_tension.fpp` **Color function advection:** diff --git a/src/simulation/include/inline_capillary.fpp b/src/simulation/include/inline_capillary.fpp deleted file mode 100644 index ba62d99740..0000000000 --- a/src/simulation/include/inline_capillary.fpp +++ /dev/null @@ -1,21 +0,0 @@ -#:def compute_capillary_stress_tensor() - Omega(1, 1) = -sigma*(w2*w2 + w3*w3)/normW - #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 - Omega(2, 1) = sigma*w1*w2/normW - Omega(1, 2) = Omega(2, 1) - - Omega(2, 2) = -sigma*(w1*w1 + w3*w3)/normW - #:endif - - if (p > 0) then - #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 - Omega(3, 1) = sigma*w1*w3/normW - Omega(1, 3) = Omega(3, 1) - - Omega(3, 2) = sigma*w2*w3/normW - Omega(2, 3) = Omega(3, 2) - - Omega(3, 3) = -sigma*(w1*w1 + w2*w2)/normW - #:endif - end if -#:enddef compute_capillary_stress_tensor diff --git a/src/simulation/include/inline_riemann.fpp b/src/simulation/include/inline_riemann.fpp index 7d8e5c74d0..b5ededfaa1 100644 --- a/src/simulation/include/inline_riemann.fpp +++ b/src/simulation/include/inline_riemann.fpp @@ -91,25 +91,3 @@ end if end if #:enddef compute_low_Mach_correction - -#:def compute_hypo_elastic_energy(EL, ER, shear_cond) - if (G_L > verysmall) then - ${EL}$ = ${EL}$ + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall) - if (${shear_cond}$) then - ${EL}$ = ${EL}$ + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall) - end if - end if - if (G_R > verysmall) then - ${ER}$ = ${ER}$ + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall) - if (${shear_cond}$) then - ${ER}$ = ${ER}$ + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall) - end if - end if -#:enddef compute_hypo_elastic_energy - -#:def compute_elastic_wave_speeds_lr() - s_L = min(vel_L(dir_idx(1)) - sqrt(max(verysmall, c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)), & - & vel_R(dir_idx(1)) - sqrt(max(verysmall, c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R))) - s_R = max(vel_R(dir_idx(1)) + sqrt(max(verysmall, c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)), & - & vel_L(dir_idx(1)) + sqrt(max(verysmall, c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L))) -#:enddef compute_elastic_wave_speeds_lr diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index b778adb642..e0db235cde 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -358,7 +358,12 @@ contains s_R = max(vel_R(dir_idx(1)) + c_fast%R, vel_L(dir_idx(1)) + c_fast%L) else if (hypoelasticity) then ! Elastic wave speed, Rodriguez et al. JCP (2019) - @:compute_elastic_wave_speeds_lr() + s_L = min(vel_L(dir_idx(1)) - f_elastic_signal_speed(c_L, G_L, tau_e_L(dir_idx_tau(1)), & + & rho_L), vel_R(dir_idx(1)) - f_elastic_signal_speed(c_R, G_R, & + & tau_e_R(dir_idx_tau(1)), rho_R)) + s_R = max(vel_R(dir_idx(1)) + f_elastic_signal_speed(c_R, G_R, tau_e_R(dir_idx_tau(1)), & + & rho_R), vel_L(dir_idx(1)) + f_elastic_signal_speed(c_L, G_L, & + & tau_e_L(dir_idx_tau(1)), rho_L)) else s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L) diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index fcc9bdf545..8aca58af48 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -1074,7 +1074,14 @@ contains if (wave_speeds == wave_speeds_direct) then #:if HYPO ! Elastic wave speed, Rodriguez et al. JCP (2019) - @:compute_elastic_wave_speeds_lr() + s_L = min(vel_L(dir_idx(1)) - f_elastic_signal_speed(c_L, G_L, & + & tau_e_L(dir_idx_tau(1)), rho_L), & + & vel_R(dir_idx(1)) - f_elastic_signal_speed(c_R, G_R, & + & tau_e_R(dir_idx_tau(1)), rho_R)) + s_R = max(vel_R(dir_idx(1)) + f_elastic_signal_speed(c_R, G_R, & + & tau_e_R(dir_idx_tau(1)), rho_R), & + & vel_L(dir_idx(1)) + f_elastic_signal_speed(c_L, G_L, & + & tau_e_L(dir_idx_tau(1)), rho_L)) s_S = (pres_R - tau_e_R(dir_idx_tau(1)) - pres_L + tau_e_L(dir_idx_tau(1)) & & + rho_L*vel_L(dir_idx(1))*(s_L - vel_L(dir_idx(1))) - rho_R*vel_R(dir_idx(1)) & & *(s_R - vel_R(dir_idx(1))))/(rho_L*(s_L - vel_L(dir_idx(1))) - rho_R*(s_R & diff --git a/src/simulation/m_riemann_solver_hypo_hlld.fpp b/src/simulation/m_riemann_solver_hypo_hlld.fpp index 885d68ece8..d7cc0f0121 100644 --- a/src/simulation/m_riemann_solver_hypo_hlld.fpp +++ b/src/simulation/m_riemann_solver_hypo_hlld.fpp @@ -7,13 +7,6 @@ #:include 'macros.fpp' #:include 'inline_riemann.fpp' -! Single source of truth for the per-component HLL flux on the (U_L, U_R, F_L, F_R) compact -! basis: used by the degenerate-fan fallback and the ADC blend, which must stay consistent. -! Textual inlining; codegen is identical to the materialized expression. -#:def hll_flux_component(LHS, I) - ${LHS}$ = (S_R*F_L(${I}$) - S_L*F_R(${I}$) + S_L*S_R*(U_R(${I}$) - U_L(${I}$)))/(S_R - S_L + verysmall) -#:enddef - module m_riemann_solver_hypo_hlld use m_derived_types @@ -54,6 +47,19 @@ contains end function f_hlld_wave_zone + !> Per-component HLL flux on the compact (U_L, U_R, F_L, F_R) basis. Shared by the degenerate-fan fallback and the ADC blend, + !! which must stay consistent. + function f_hll_flux(S_L, S_R, F_L_i, F_R_i, U_L_i, U_R_i) result(flux) + + $:GPU_ROUTINE(function_name='f_hll_flux', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: S_L, S_R, F_L_i, F_R_i, U_L_i, U_R_i + real(wp) :: flux + + flux = (S_R*F_L_i - S_L*F_R_i + S_L*S_R*(U_R_i - U_L_i))/(S_R - S_L + verysmall) + + end function f_hll_flux + !> HLLD Riemann solver resolves all 5 waves for the hypoelastic equations: 1 entropy wave, 2 shear stress waves, 2 fast waves. subroutine s_hypo_hlld_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, qL_prim_vf, & & qR_prim_rsx_vf, dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, & @@ -157,7 +163,7 @@ contains real(wp) :: phi real(wp), parameter :: ADC_power = 1.0_wp real(wp) :: alpha_L_sum, alpha_R_sum - logical :: degenerate, shear_degenerate, fan_fallback + logical :: degenerate, shear_degenerate, fan_fallback, shear_cond integer :: i, j, k, l, ipass, zone call s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, & @@ -180,7 +186,7 @@ contains ! (_hlld_p1..p4) are ONLY for source readability; fypp concatenates them into one ! clause below. That wrapping is FOLD_DIRECTIVE's job -- its within-clause comma split ! exists for exactly this case -- not the fragments'. - #:set _hlld_p1 = '[i,j,k,l,ipass,degenerate,shear_degenerate,fan_fallback,alpha_rho_L,alpha_rho_R,vel,alpha_L,alpha_R,rho,pres,E,gamma,pi_inf,qv,vel_rms,c,S_L,S_R,s_M,S_Lstar,S_Rstar,pTot_L,pTot_R,rhoL_star,rhoR_star,U_L,U_R,F_L,F_R,F_hlld,us_c,uss_c,zone,F_HLL_c,U_HLL_c,rho_HLL,u_n_HLL_cons,tau_nn_HLL,u_n_HLL_trace,u_t_HLL_trace,p_face_HLL,tau_qq_face_HLL,ncomp,G_eff,G_eff_tol,C_NC,sqrtC_NC,A_L,A_R,denomA,fac_L,fac_R,' + #:set _hlld_p1 = '[i,j,k,l,ipass,degenerate,shear_degenerate,fan_fallback,shear_cond,alpha_rho_L,alpha_rho_R,vel,alpha_L,alpha_R,rho,pres,E,gamma,pi_inf,qv,vel_rms,c,S_L,S_R,s_M,S_Lstar,S_Rstar,pTot_L,pTot_R,rhoL_star,rhoR_star,U_L,U_R,F_L,F_R,F_hlld,us_c,uss_c,zone,F_HLL_c,U_HLL_c,rho_HLL,u_n_HLL_cons,tau_nn_HLL,u_n_HLL_trace,u_t_HLL_trace,p_face_HLL,tau_qq_face_HLL,ncomp,G_eff,G_eff_tol,C_NC,sqrtC_NC,A_L,A_R,denomA,fac_L,fac_R,' #:set _hlld_p2 = 'u_n_L,u_t_L,u_n_R,u_t_R,u_t2_L,u_t2_R,tau_nn_L,tau_nt_L,tau_tt_L,tau_nn_R,tau_nt_R,tau_tt_R,tau_nt2_L,tau_nt2_R,tau_t2t2_L,tau_t2t2_R,tau_t1t2_L,tau_t1t2_R,tau_qq_L,tau_qq_R,G_L,G_R,tau_e_L,tau_e_R,alpha1_L_star,alpha1_R_star,alpha2_L_star,alpha2_R_star,u_t_star,tau_nt_star,u_t2_star,tau_nt2_star,tau_nn_L_star,tau_nn_R_star,tau_tt_L_star,tau_tt_R_star,tau_tt_L_starstar,tau_tt_R_starstar,' #:set _hlld_p3 = 'tau_t2t2_L_star,tau_t2t2_R_star,tau_t2t2_L_starstar,tau_t2t2_R_starstar,tau_t1t2_L_star,tau_t1t2_R_star,tau_t1t2_L_starstar,tau_t1t2_R_starstar,tau_qq_L_star,tau_qq_R_star,pTot_star,E_L_star,E_R_star,E_L_starstar,E_R_starstar,p_face,tau_qq_face,u_n_face,u_t_face,G_hat,rho_hat,tau_nn_hat,tau_nt_hat,tau_tt_hat,tau_qq_hat,tau_nt2_hat,tau_t2t2_hat,tau_t1t2_hat,' #:set _hlld_p4 = 'alpha_hat,alpha_rho_hat,tau_e_hat,pres_hat,blkmod1_hat,blkmod2_hat,K_hat,C_hat_1,C_hat_2,Sigma_L,Sigma_R,dSigma,Sigma_ref,a_L_ref,a_R_ref,a_ref,du_t,dtau_nt,du_t2,dtau_nt2,sensor_ptot,sensor_vt,sensor_tnt,sensor_combined,phi,alpha_L_sum,alpha_R_sum]' @@ -335,9 +341,9 @@ contains $:GPU_LOOP(parallelism='[seq]') do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1 - @:compute_hypo_elastic_energy(E%L, E%R, & - & (n > 0 .and. p == 0 .and. i == 2) .or. (p > 0 .and. (i == 2 & - & .or. i == 4 .or. i == 5))) + shear_cond = (n > 0 .and. p == 0 .and. i == 2) .or. (p > 0 .and. (i == 2 .or. i == 4 .or. i == 5)) + E%L = E%L + f_elastic_energy(tau_e_L(i), G_L, shear_cond) + E%R = E%R + f_elastic_energy(tau_e_R(i), G_R, shear_cond) end do ! Compute Riemann states @@ -612,7 +618,7 @@ contains ! HLL (or one-sided) fallback for an invalid wave structure if (S_L < 0._wp .and. S_R > 0._wp) then do i = 1, ncomp - @:hll_flux_component(F_hlld(i), i) + F_hlld(i) = f_hll_flux(S_L, S_R, F_L(i), F_R(i), U_L(i), U_R(i)) end do else if (S_L >= 0._wp) then F_hlld(1:ncomp) = F_L(1:ncomp) @@ -746,7 +752,7 @@ contains ! phi is anchor-independent: computed once in the shared section above if (S_L < 0._wp .and. S_R > 0._wp) then do i = 1, ncomp - @:hll_flux_component(F_HLL_c, i) + F_HLL_c = f_hll_flux(S_L, S_R, F_L(i), F_R(i), U_L(i), U_R(i)) F_hlld(i) = F_HLL_c + phi*(F_hlld(i) - F_HLL_c) end do else @@ -784,7 +790,7 @@ contains U_HLL_c = (S_R*U_R(11) - S_L*U_L(11) - (F_R(11) - F_L(11)))/(S_R - S_L + verysmall) tau_qq_face_HLL = U_HLL_c/(rho_HLL + verysmall) ! This branch implies S_L < 0 < S_R, so component 3 of F_HLL is the interior HLL flux - @:hll_flux_component(F_HLL_c, 3) + F_HLL_c = f_hll_flux(S_L, S_R, F_L(3), F_R(3), U_L(3), U_R(3)) p_face_HLL = F_HLL_c - rho_HLL*u_n_HLL_cons*u_n_HLL_cons + tau_nn_HLL end if end if diff --git a/src/simulation/m_riemann_state.fpp b/src/simulation/m_riemann_state.fpp index bd978308f2..9edcd1cf32 100644 --- a/src/simulation/m_riemann_state.fpp +++ b/src/simulation/m_riemann_state.fpp @@ -73,6 +73,37 @@ module m_riemann_state contains + !> Elastic signal speed of Rodriguez et al. JCP (2019): the acoustic speed stiffened by the shear modulus and the normal elastic + !! stress. Callers subtract it for the left-going wave and add it for the right-going one. + function f_elastic_signal_speed(c, G, tau, rho) result(a) + + $:GPU_ROUTINE(function_name='f_elastic_signal_speed', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: c, G, tau, rho + real(wp) :: a + + a = sqrt(max(verysmall, c*c + (((4._wp*G)/3._wp) + tau)/rho)) + + end function f_elastic_signal_speed + + !> Elastic strain energy carried by one stress component, doubled for a shear component because the tensor holds it once but the + !! energy counts both off-diagonal entries. Zero where the material has no shear modulus. + function f_elastic_energy(tau, G, is_shear) result(dE) + + $:GPU_ROUTINE(function_name='f_elastic_energy', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: tau, G + logical, intent(in) :: is_shear + real(wp) :: dE + + dE = 0._wp + if (G > verysmall) then + dE = (tau*tau)/max(4._wp*G, verysmall) + if (is_shear) dE = dE + (tau*tau)/max(4._wp*G, verysmall) + end if + + end function f_elastic_energy + !> Dispatch to the subroutines that are utilized to compute the viscous source fluxes for either Cartesian or cylindrical !! geometries. For more information please refer to: 1) s_compute_cartesian_viscous_source_flux 2) !! s_compute_cylindrical_viscous_source_flux diff --git a/src/simulation/m_surface_tension.fpp b/src/simulation/m_surface_tension.fpp index 0b7caf7a47..2fbe62261b 100644 --- a/src/simulation/m_surface_tension.fpp +++ b/src/simulation/m_surface_tension.fpp @@ -4,7 +4,6 @@ #:include 'case.fpp' #:include 'macros.fpp' -#:include 'inline_capillary.fpp' !> @brief Computes capillary source fluxes and color-function gradients for the diffuse-interface surface tension model module m_surface_tension @@ -58,6 +57,41 @@ contains end subroutine s_initialize_surface_tension_module !> Compute the capillary source flux from reconstructed color-gradient fields + !> Capillary stress tensor of Schmidmayer et al. JCP (2017) for one face. Omega is intent(inout): the entries a given + !! dimensionality does not define are left as the caller had them, which is what the macro this replaced did. + subroutine s_compute_capillary_stress_tensor(sigma_c, w1, w2, w3, normW, Omega) + + $:GPU_ROUTINE(function_name='s_compute_capillary_stress_tensor', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: sigma_c, w1, w2, w3, normW + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3, 3), intent(inout) :: Omega + #:else + real(wp), dimension(num_dims, num_dims), intent(inout) :: Omega + #:endif + + Omega(1, 1) = -sigma_c*(w2*w2 + w3*w3)/normW + #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 + Omega(2, 1) = sigma_c*w1*w2/normW + Omega(1, 2) = Omega(2, 1) + + Omega(2, 2) = -sigma_c*(w1*w1 + w3*w3)/normW + #:endif + + if (p > 0) then + #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 + Omega(3, 1) = sigma_c*w1*w3/normW + Omega(1, 3) = Omega(3, 1) + + Omega(3, 2) = sigma_c*w2*w3/normW + Omega(2, 3) = Omega(3, 2) + + Omega(3, 3) = -sigma_c*(w1*w1 + w2*w2)/normW + #:endif + end if + + end subroutine s_compute_capillary_stress_tensor + subroutine s_compute_capillary_source_flux(vSrc_rsx_vf, flux_src_vf, id, isx, isy, isz) real(wp), dimension(-1:,-1:,-1:,1:), intent(in) :: vSrc_rsx_vf @@ -98,7 +132,7 @@ contains normW = (normWL + normWR)/2._wp if (normW > capillary_cutoff) then - @:compute_capillary_stress_tensor() + call s_compute_capillary_stress_tensor(sigma, w1, w2, w3, normW, Omega) do i = 1, num_dims flux_src_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, & @@ -141,7 +175,7 @@ contains normW = (normWL + normWR)/2._wp if (normW > capillary_cutoff) then - @:compute_capillary_stress_tensor() + call s_compute_capillary_stress_tensor(sigma, w1, w2, w3, normW, Omega) do i = 1, num_dims flux_src_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + i - 1)%sf(j, & @@ -184,7 +218,7 @@ contains normW = (normWL + normWR)/2._wp if (normW > capillary_cutoff) then - @:compute_capillary_stress_tensor() + call s_compute_capillary_stress_tensor(sigma, w1, w2, w3, normW, Omega) do i = 1, num_dims flux_src_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + i - 1)%sf(j, & diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp index b34390d02b..981cf25e07 100644 --- a/src/simulation/m_viscous.fpp +++ b/src/simulation/m_viscous.fpp @@ -8,34 +8,6 @@ #! Herschel-Bulkley mixture inv-Re from the local strain rate; Newtonian path is #! the original alpha-weighted harmonic average. Component-3 strain terms are #! gated for compile validity (num_dims <= 2) and runtime safety (p > 0). -#:def compute_axis_inv_re() - if (any_non_newtonian) then - gamma_dot_c = f_compute_shear_rate_from_components(grad_x_vf(1)%sf(j, k, l), grad_y_vf(2)%sf(j, k, l), 0._wp, & - & 0.5_wp*(grad_y_vf(1)%sf(j, k, l) + grad_x_vf(2)%sf(j, k, l)), 0._wp, 0._wp) - #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 - if (p > 0) then - gamma_dot_c = f_compute_shear_rate_from_components(grad_x_vf(1)%sf(j, k, l), grad_y_vf(2)%sf(j, k, l), & - & grad_z_vf(3)%sf(j, k, l), 0.5_wp*(grad_y_vf(1)%sf(j, k, l) + grad_x_vf(2)%sf(j, k, l)), & - & 0.5_wp*(grad_z_vf(1)%sf(j, k, l) + grad_x_vf(3)%sf(j, k, l)), 0.5_wp*(grad_z_vf(2)%sf(j, k, & - & l) + grad_y_vf(3)%sf(j, k, l))) - end if - #:endif - call s_compute_mixture_inv_re(alpha_visc, gamma_dot_c, Res_viscous, Re_visc) - else - $:GPU_LOOP(parallelism='[seq]') - do i = 1, 2 - Re_visc(i) = dflt_real - - if (Re_size(i) > 0) Re_visc(i) = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do q = 1, Re_size(i) - Re_visc(i) = alpha_visc(Re_idx(i, q))/Res_viscous(i, q) + Re_visc(i) - end do - - Re_visc(i) = 1._wp/max(Re_visc(i), sgm_eps) - end do - end if -#:enddef !> @brief Computes viscous stress tensors and diffusive flux contributions for the Navier--Stokes equations module m_viscous @@ -80,6 +52,52 @@ contains end subroutine s_initialize_viscous_module !> Compute viscous stress tensor near cylindrical axis, avoiding 1/r singularity at y_cb(-1)=0 + !> Inverse Reynolds numbers of the mixture at one cell. Non-Newtonian fluids need the local shear rate, which comes from the + !! velocity gradients; Newtonian ones reduce to the volume-fraction-weighted harmonic mean. + subroutine s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc) + + $:GPU_ROUTINE(function_name='s_compute_axis_inv_re', parallelism='[seq]', cray_inline=True) + + type(scalar_field), dimension(num_dims), intent(in) :: grad_x_vf, grad_y_vf, grad_z_vf + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in) :: alpha_visc + #:else + real(wp), dimension(num_fluids), intent(in) :: alpha_visc + #:endif + integer, intent(in) :: j, k, l + real(wp), dimension(2), intent(out) :: Re_visc + real(wp) :: gamma_dot_c + integer :: i, q + + if (any_non_newtonian) then + gamma_dot_c = f_compute_shear_rate_from_components(grad_x_vf(1)%sf(j, k, l), grad_y_vf(2)%sf(j, k, l), 0._wp, & + & 0.5_wp*(grad_y_vf(1)%sf(j, k, l) + grad_x_vf(2)%sf(j, k, l)), 0._wp, 0._wp) + #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 + if (p > 0) then + gamma_dot_c = f_compute_shear_rate_from_components(grad_x_vf(1)%sf(j, k, l), grad_y_vf(2)%sf(j, k, l), & + & grad_z_vf(3)%sf(j, k, l), 0.5_wp*(grad_y_vf(1)%sf(j, k, l) + grad_x_vf(2)%sf(j, k, l)), & + & 0.5_wp*(grad_z_vf(1)%sf(j, k, l) + grad_x_vf(3)%sf(j, k, l)), 0.5_wp*(grad_z_vf(2)%sf(j, k, & + & l) + grad_y_vf(3)%sf(j, k, l))) + end if + #:endif + call s_compute_mixture_inv_re(alpha_visc, gamma_dot_c, Res_viscous, Re_visc) + else + $:GPU_LOOP(parallelism='[seq]') + do i = 1, 2 + Re_visc(i) = dflt_real + + if (Re_size(i) > 0) Re_visc(i) = 0._wp + $:GPU_LOOP(parallelism='[seq]') + do q = 1, Re_size(i) + Re_visc(i) = alpha_visc(Re_idx(i, q))/Res_viscous(i, q) + Re_visc(i) + end do + + Re_visc(i) = 1._wp/max(Re_visc(i), sgm_eps) + end do + end if + + end subroutine s_compute_axis_inv_re + subroutine s_compute_viscous_stress_cylindrical_boundary(q_prim_vf, grad_x_vf, grad_y_vf, grad_z_vf, tau_Re_vf, ix, iy, iz) type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf @@ -88,7 +106,6 @@ contains type(int_bounds_info), intent(in) :: ix, iy, iz real(wp) :: rho_visc, gamma_visc, pi_inf_visc, qv_visc, alpha_visc_sum !< Mixture variables real(wp), dimension(2) :: Re_visc - real(wp) :: gamma_dot_c !< Effective shear rate for non-Newtonian mixture inv-Re. #:if not MFC_CASE_OPTIMIZATION and USING_AMD real(wp), dimension(3) :: alpha_visc, alpha_rho_visc @@ -120,7 +137,7 @@ contains #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 if (shear_stress) then ! Shear stresses $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, q, rho_visc, gamma_visc, pi_inf_visc, qv_visc, & - & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re, gamma_dot_c]') + & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re]') do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 do j = is1_viscous%beg, is1_viscous%end @@ -147,7 +164,7 @@ contains call s_compute_mixture_coefficients(alpha_rho_visc, alpha_visc, rho_visc, gamma_visc, pi_inf_visc, & & qv_visc) if (viscous) then - @:compute_axis_inv_re() + call s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc) end if ! Shear stress near cylindrical axis: includes v/r hoop term @@ -174,7 +191,7 @@ contains #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 if (bulk_stress) then ! Bulk stresses $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, q, rho_visc, gamma_visc, pi_inf_visc, qv_visc, & - & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re, gamma_dot_c]') + & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re]') do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 do j = is1_viscous%beg, is1_viscous%end @@ -201,7 +218,7 @@ contains call s_compute_mixture_coefficients(alpha_rho_visc, alpha_visc, rho_visc, gamma_visc, pi_inf_visc, & & qv_visc) if (viscous) then - @:compute_axis_inv_re() + call s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc) end if tau_Re(2, 2) = (grad_x_vf(1)%sf(j, k, l) + grad_y_vf(2)%sf(j, k, & @@ -222,7 +239,7 @@ contains #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 if (shear_stress) then ! Shear stresses $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, q, rho_visc, gamma_visc, pi_inf_visc, qv_visc, & - & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re, gamma_dot_c]') + & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re]') do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 do j = is1_viscous%beg, is1_viscous%end @@ -249,7 +266,7 @@ contains call s_compute_mixture_coefficients(alpha_rho_visc, alpha_visc, rho_visc, gamma_visc, pi_inf_visc, & & qv_visc) if (viscous) then - @:compute_axis_inv_re() + call s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc) end if tau_Re(2, 2) = -(2._wp/3._wp)*grad_z_vf(3)%sf(j, k, l)/y_cc(k)/Re_visc(1) @@ -273,7 +290,7 @@ contains if (bulk_stress) then ! Bulk stresses $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, q, rho_visc, gamma_visc, pi_inf_visc, qv_visc, & - & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re, gamma_dot_c]') + & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re]') do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 do j = is1_viscous%beg, is1_viscous%end @@ -300,7 +317,7 @@ contains call s_compute_mixture_coefficients(alpha_rho_visc, alpha_visc, rho_visc, gamma_visc, pi_inf_visc, & & qv_visc) if (viscous) then - @:compute_axis_inv_re() + call s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc) end if tau_Re(2, 2) = grad_z_vf(3)%sf(j, k, l)/y_cc(k)/Re_visc(2)