Replace three enclosing-scope fypp macros with explicit procedures

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:**
 

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

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)

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 &

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

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

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, &