HLLC: opt-in smoothing of the tangential interface velocity in the viscous work (ib_smooth_contact)

docs/documentation/case.md

@@ -361,6 +361,7 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi
 | `ib_flux_wall`        | Logical | Fixed reacting bodies: the mass, species and energy fluxes through the faces between fluid and body cells are the surface solve's rates, so the face fluxes carry the surface rates, and the gas budget closes to within the time and surface-solve error (examples/2D_ibm_wall_mass_budget). |
 | `ib_ip_interp`        | Integer | Image-point interpolation: 1 inverse-distance weights over 2x2(x2) cells (default), 2 triangular-shaped-cloud (TSC) weights over 3x3(x3) cells. |
 | `ib_wall_order`       | Integer | Ghost velocity at no-slip walls: 1 default (wall value, with the `ib_second_order_vel` mirror if set), 2 two-probe quadratic wall-normal fill. |
+| `ib_smooth_contact`   | Logical | HLLC only: blend the tangential interface velocity of the viscous work across the contact wave (the sign of the contact speed is arbitrary where it is round-off). |
 | `ib_visc_compact`     | Logical | With `weno_Re_flux`, compact velocity gradients on the viscous faces within two cells of a body. |
 | `geometry`           | Integer | Geometry configuration of the patch.|
 | `x[y,z]_centroid`    | Real    | Centroid of the applied geometry in the [x,y,z]-direction. |
@@ -479,6 +480,7 @@ Additional details on this specification can be found in [NACA airfoil](https://
 - `ib_ip_interp` selects how the state at a ghost point's image point is interpolated from the grid. `1` (default) is the inverse-distance weighting over the 2x2(x2) cells around the point, which drops solid cells from the stencil. `2` uses triangular-shaped-cloud (quadratic B-spline) weights over the 3x3(x3) cells around the nearest cell (Ezra, Peles & Kozak 2025, Int. J. Heat Mass Transf. 252:127314). The kernel is a partition of unity and smooth across cell boundaries, so a wall quantity interpolated from it converges one order faster than from inverse-distance weights, which are only continuous. The TSC stencil reads solid cells, which carry the ghost values of the last correction; a point buried in a body gets no fluid weight and takes the neighbor average. It does not combine with `ib_second_order_vel`.
 - `ib_wall_order = 2` fills the ghost velocity of each no-slip wall point from a quadratic wall-normal profile through the wall value and two probes at one and two cell widths from the wall along the surface normal, interpolated with the same TSC weights as `ib_ip_interp = 2` (the option works with either `ib_ip_interp` setting; it does not require 2). The ghost value is what enforces no-slip in the solve. Where the outer probe has too little fluid around it (a narrow gap, a rim, a neighbouring body) a ghost point falls back to the linear mirror, or to the default fill when the first probe also lacks fluid; the share of each tier is printed at the end of the run (`[WALLFILL]`, for every case with `ib_wall_order = 2`). Ghost cells keep the values of the last correction while the option is on, which is also what `ib_ip_interp = 2` does. Slip walls, the pressure and the energy (temperature) of the ghost are unchanged, and the option does not combine with `ib_second_order_vel`. The default `1` leaves every result unchanged.
   Measured gain (non-reacting ideal gas, tip or peak Mach 0.05, 4 to 80 cells across the gap, NVHPC on an A100): between flat walls (plane Poiseuille) the velocity and flow-rate errors are 10 to 20 times smaller than with the default and converge at about order 1.5 to 1.9 until they reach the level of a domain no-slip boundary, whatever the wall position inside a cell; the default converges at order 1 with the no-slip acting about 0.75 cell outside the wall. In a Taylor-Couette annulus (curved walls, one rotating) the gain is 25 to 40 % of the velocity error and a factor 2 or more in the torque, and neither option converges monotonically: an error of about 0.5 % of the torque and 0.3 % of the wall speed remains at 80 cells per inner radius, and the rotating wall keeps a velocity deficit that this option halves but does not remove. With 4 to 6 cells across the gap 15 to 32 % of the ghost points use the linear mirror and the quadratic fill is not accurate (26 % velocity error at 4 cells across a flat channel). The wall shear read back from a fit to the near-wall cells does not converge for either option; use the flow rate or an integrated force instead.
+- `ib_smooth_contact = T` (HLLC, `viscous`) blends the left and right states of the tangential interface velocity in the viscous work over `|s_S| < 0.01 c` instead of switching on the sign of the contact speed `s_S`. A face between a fluid cell and a ghost cell whose `s_S` is round-off would otherwise pick one side arbitrarily, although the two tangential velocities differ by order one across the ghost mirror, so builds and devices could disagree there. Faces away from `s_S = 0` are unchanged. The default `F` leaves every result unchanged.
 - `ib_visc_compact = T` (requires `weno_Re_flux`) replaces the WENO velocity gradient on the viscous faces within two cells of a body cell by the compact one of `weno_Re_flux = F`: the normal difference across the face and the mean of the two cells' central tangential differences. The gradient of `weno_Re_flux` is built from reconstructed face values, a divergence-theorem cell gradient and a second reconstruction, a stencil about five cells wide, while the ghost values extend the flow only a layer or two into the body (the mirror leaves a kink at the wall), so the wide stencil carries the kink into the fluid faces next to the wall. Faces farther than two cells from a body are unchanged. The default `F` leaves every result unchanged.
 
 - `ib_neighborhood_radius` controls the size of the neighborhood size. A value of $r$ indicates that any given rank is aware of IBs up to $r$ ranks away. This value defaults to 0, which leaves the radius unset so that it is selected automatically. This parameter is required to strong-scale a case when IBs eventually grow to be larger than one full processor domain wide.

src/common/m_constants.fpp

@@ -24,6 +24,10 @@ module m_constants
     real(wp), parameter :: T_surface_min = 200._wp
     real(wp), parameter :: T_surface_max = 5000._wp
 
+    !> Half-width of the contact-speed ramp, as a fraction of the sound speed, over which the HLLC tangential interface velocity
+    !! blends between the left and right states (ib_smooth_contact)
+    real(wp), parameter :: contact_ramp = 1.e-2_wp
+
     !> Reals per IB in restart_data/ib_state: time, force(3), torque(3), vel(3), angular_vel(3), angles(3), centroid(3), radius,
     !! Twall. Twall carries a thermal_bc = 3 body's evolving temperature across a restart.
     integer, parameter :: ib_state_nfields = 21

src/simulation/m_global_parameters.fpp

@@ -673,6 +673,7 @@ contains
         ib_flux_wall = .false.
         ib_ip_interp = 1
         ib_wall_order = 1
+        ib_smooth_contact = .false.
         ib_visc_compact = .false.
 
         ! Bubble modeling (sim-specific)

src/simulation/m_ibm.fpp

@@ -212,7 +212,7 @@ contains
         ib_ip_width = merge(3, 2, ib_ip_interp == 2)
         ib_wall_quad = ib_wall_order == 2
         ib_keep_ghost = ib_ip_interp == 2 .or. ib_wall_quad
-        $:GPU_UPDATE(device='[ib_second_order_vel, ib_ip_min_dist, ib_ip_interp, ib_ip_width, ib_keep_ghost, ib_wall_quad]')
+        $:GPU_UPDATE(device='[ib_second_order_vel, ib_ip_min_dist, ib_ip_interp, ib_ip_width, ib_keep_ghost, ib_wall_quad, ib_smooth_contact]')
 
         if (collision_model > 0) call s_initialize_collisions_module()
 

src/simulation/m_riemann_solver_hllc.fpp

@@ -85,7 +85,7 @@ contains
             real(wp)                             :: s_L, s_R, s_M, s_P, s_S
             real(wp)                             :: xi_L, xi_R        !< Left and right wave speeds functions
             real(wp)                             :: xi_L_m1, xi_R_m1  !< xi_L/R - 1, computed without cancellation
-            real(wp)                             :: xi_M, xi_P
+            real(wp)                             :: xi_M, xi_P, xi_tan
             real(wp)                             :: xi_MP, xi_PP
             #:if not MFC_CASE_OPTIMIZATION and USING_AMD
                 real(wp), dimension(3) :: R0_L, R0_R
@@ -876,7 +876,7 @@ contains
                             #:set _hllc_s2 = ['pi_inf_R', 'qv_R', 'alpha_L_sum', 'alpha_R_sum', 'E_L', 'E_R', 'MW_L', 'MW_R', &
                                 & 'R_gas_L', 'R_gas_R', 'Cp_L', 'Cp_R', 'Cv_L', 'Cv_R', 'c_sum_Yi_Phi']
                             #:set _hllc_s3 = ['Gamm_L', 'Gamm_R', 'Y_L', 'Y_R', 'H_L', 'H_R', 'qv_avg', 'rho_avg', 'gamma_avg', &
-                                & 'H_avg', 'c_L', 'c_R', 'c_avg', 's_P', 's_M', 'xi_P', 'xi_M', 'xi_L']
+                                & 'H_avg', 'c_L', 'c_R', 'c_avg', 's_P', 's_M', 'xi_P', 'xi_M', 'xi_tan', 'xi_L']
                             #:set _hllc_s4 = ['xi_R', 'xi_L_m1', 'xi_R_m1', 'Ms_L', 'Ms_R', 'pres_SL', 'pres_SR', 's_L', 's_R', &
                                 & 's_S', 'vel_avg_rms', 'pcorr']
                             #:set _hllc_s5 = ['ib_fw', 'fw_o', 'fw_a', 'fw_m']
@@ -1401,6 +1401,15 @@ contains
                                                            & *s_M*xi_L_m1) + xi_P*(vel_R(dir_idx(i)) + dir_flg(dir_idx(i)) &
                                                            & *s_P*xi_R_m1)
                                         end do
+                                        if (ib_smooth_contact) then
+                                            xi_tan = f_contact_weight(s_S, max(c_L, c_R))
+                                            $:GPU_LOOP(parallelism='[seq]')
+                                            do i = 2, num_dims
+                                                vel_src_rsx_vf(${SF('')}$, &
+                                                               & dir_idx(i)) = xi_tan*vel_L(dir_idx(i)) + (1._wp - xi_tan) &
+                                                               & *vel_R(dir_idx(i))
+                                            end do
+                                        end if
 
                                         ! COLOR FUNCTION FLUX
                                         if (surface_tension) then

src/simulation/m_riemann_state.fpp

@@ -11,7 +11,7 @@ module m_riemann_state
 
     use m_derived_types
     use m_global_parameters
-    use m_constants, only: verysmall
+    use m_constants, only: verysmall, contact_ramp, sgm_eps
     use m_hb_function
 
     implicit none
@@ -119,6 +119,18 @@ contains
 
     end function f_low_Mach_pcorr_hll
 
+    !> Weight on the left state of the contact-upwinded tangential interface velocity: the step sign(s_S) smoothed to a linear ramp
+    !! of half-width contact_ramp*c, so a face with s_S at round-off averages the two sides rather than picking one.
+    pure function f_contact_weight(s_S, c) result(w)
+
+        $:GPU_ROUTINE(function_name='f_contact_weight', parallelism='[seq]', cray_inline=True)
+        real(wp), intent(in) :: s_S, c
+        real(wp)             :: w
+
+        w = min(max(5.e-1_wp + 5.e-1_wp*s_S/max(contact_ramp*c, sgm_eps), 0._wp), 1._wp)
+
+    end function f_contact_weight
+
     !> The same correction for the HLLC flux, where the star state supplies the pressure jump directly and the correction scales
     !! with the mass flux through the acoustic waves instead. Zero unless low_Mach == 1.
     function f_low_Mach_pcorr_hllc(vel_L_rms, vel_R_rms, c_L, c_R, rho_L, rho_R, s_L, s_R, vel_L_norm, vel_R_norm) result(pcorr)