IBM: restore the IB force and torque with the dual-pass hypoelastic HLLD

docs/documentation/case.md

@@ -434,7 +434,7 @@ Additional details on this specification can be found in [NACA airfoil](https://
 
   with \f$m\f$ = `rho_solid`\f$\cdot V\f$ (\f$V\f$ from the geometry: circle per unit depth, sphere, or cylinder), \f$c_s\f$ = `cp_solid`, \f$P\f$ = `heat_power`, \f$\varepsilon\f$ = `emissivity` and \f$T_\mathrm{rad}\f$ = `T_rad`. The solid energy is \f$m c_s (T - T_\mathrm{ref})\f$ with \f$T_\mathrm{ref}\f$ = 298.15 K, the reference of the formation enthalpies, so the energy the solid's lost mass carries is consistent with the gas energy. \f$\dot Q_\mathrm{in}\f$ is the total energy flux into the body summed over every face between its cells and fluid cells, taken from the same Riemann, viscous, conductive and diffusive fluxes the flow update applies and weighted by the Runge-Kutta stages, so over a step the body gains exactly the energy those fluxes remove from the fluid. \f$\dot M_\mathrm{out}\f$ is the body's gasification rate from the surface solve (\f$\sum\f$ `area`\f$\cdot\f$`mdot` over the `ib_surface_wrt` surface points; 0 for an inert surface), not the face mass flux, which also carries the mass the ghost-cell reconstruction does not conserve. \f$A\f$ is the surface area from the same surface points. \f$T\f$ is advanced once per time step, explicitly, which is stable while the body's thermal time \f$m c_s/(hA)\f$ greatly exceeds the flow time step. The temperature is clamped to \f$[200, 5000]\f$ K (`T_surface_min`, `T_surface_max`) and is carried across restarts in `restart_data/ib_state`.
 
-  Limitations. One temperature per body is valid only while the Biot number \f$hR/k_s\f$ is small (graphite particles, millimetre rods). The geometry and \f$m\f$ are fixed: mass lost enters the energy balance but the body does not recede. The balance is exact in the fluxes the right-hand side applies, not in the total energy of the discretised system, because the ghost-cell reconstruction is not itself conservative; the energy error this leaves has not been quantified, and the balance has been checked against a first-order estimate, not an analytic solution. It needs `chemistry = T` and a circle (2), sphere (8) or cylinder (10) that is not moving (`moving_ibm = 0`), and is not supported with `igr` or hypoelastic HLLD, which form no single set of face fluxes. `emissivity` \f$> 0\f$ needs `T_rad` \f$> 0\f$.
+  Limitations. One temperature per body is valid only while the Biot number \f$hR/k_s\f$ is small (graphite particles, millimetre rods). The geometry and \f$m\f$ are fixed: mass lost enters the energy balance but the body does not recede. The balance is exact in the fluxes the right-hand side applies, not in the total energy of the discretised system, because the ghost-cell reconstruction is not itself conservative; the energy error this leaves has not been quantified, and the balance has been checked against a first-order estimate, not an analytic solution. It needs `chemistry = T` and a circle (2), sphere (8) or cylinder (10) that is not moving (`moving_ibm = 0`). `emissivity` \f$> 0\f$ needs `T_rad` \f$> 0\f$.
 
 - `Twall` specifies the prescribed surface temperature when `thermal_bc = 1` (the initial temperature when `thermal_bc = 3`) and must be positive in that case. With chemistry it must also be within `T_surface_min` = 200 K and `T_surface_max` = 5000 K (`m_constants.fpp`), the range of the NASA polynomial fits.
 
@@ -894,7 +894,7 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu
 | 15–17   | angles about x/y/z |
 | 18–20   | centroid, x/y/z |
 
-The force and torque are the momentum the fluid exchanges with the body: the inviscid and viscous fluxes of the latest RK stage, summed over the faces between body cells and fluid cells, with the torque taken about the centroid using each face's centre. Collision and body forces (`bf_x` etc. times `mass`) are added to them.
+The force and torque are the momentum the fluid exchanges with the body: the inviscid and viscous fluxes of the latest RK stage, summed over the faces between body cells and fluid cells, with the torque taken about the centroid using each face's centre. These are the fluxes the right-hand side differences; with hypoelastic HLLD, which assembles two flux sets, the hat-L flux at a fluid cell's right face and the hat-R flux at its left. Collision and body forces (`bf_x` etc. times `mass`) are added to them.
 
 The file carries no header line, because every record sits at a computed byte offset and a header would shift them all. Each record is exactly 353 bytes including its newline (`I10` followed by nineteen `1X,ES17.9E3` fields), so the whole file loads with `numpy.loadtxt` and a single body or step can be read without scanning it:
 
@@ -922,7 +922,7 @@ Rows are written in global body-id order, so the file is byte-identical however
 | 5       | heat rate Q [W], positive into the body |
 | 6       | `Twall` [K] |
 
-The gasification rate is the surface-solve mass flux integrated over the body, \f$\sum\f$ `area`\f$\cdot\f$`mdot` of the `ib_surface_wrt` points (0 for an inert body): the body's mass loss. Columns 4 and 5 are the step's mass out of and energy into the body summed over its faces, as for `thermal_bc = 3` above; for a `thermal_bc = 3` body, Q also includes \f$c_s (T - 298.15\,\mathrm{K})\f$ times column 3, so there Q is exactly the fluid's part of \f$m c_s\, dT/dt\f$ (add `heat_power` and radiation for the whole balance). The face mass outflow is not the gasification rate: it also carries the mass the ghost-cell reconstruction does not conserve, which on a reacting cylinder at 20 cells per diameter can be a net inflow several times the gasification rate. As for the force, a row carries the exchange over the previous step and `Twall` at its end, so with `ib_force_stride = 1` \f$\sum\f$ (column 3)\f$\cdot\Delta t\f$ over the rows approximates the mass the body lost (the gasification rate is the last Runge-Kutta stage's, not stage-weighted like columns 4 and 5); with a larger stride a row holds one step of every `ib_force_stride`, so the sum no longer applies. Every body is summed, so a fixed-`Twall` reacting surface gets its mass loss and a heated wall its heat transfer. `ib_heat_wrt` and `thermal_bc = 3` require `chemistry`, Cartesian coordinates and fixed bodies, and are rejected with `igr` and hypoelastic HLLD, which form no face fluxes. A restart reads `Twall` from `restart_data/ib_state` only for `thermal_bc = 3` bodies; state files from before the record grew to 21 fields cannot be restarted.
+The gasification rate is the surface-solve mass flux integrated over the body, \f$\sum\f$ `area`\f$\cdot\f$`mdot` of the `ib_surface_wrt` points (0 for an inert body): the body's mass loss. Columns 4 and 5 are the step's mass out of and energy into the body summed over its faces, as for `thermal_bc = 3` above; for a `thermal_bc = 3` body, Q also includes \f$c_s (T - 298.15\,\mathrm{K})\f$ times column 3, so there Q is exactly the fluid's part of \f$m c_s\, dT/dt\f$ (add `heat_power` and radiation for the whole balance). The face mass outflow is not the gasification rate: it also carries the mass the ghost-cell reconstruction does not conserve, which on a reacting cylinder at 20 cells per diameter can be a net inflow several times the gasification rate. As for the force, a row carries the exchange over the previous step and `Twall` at its end, so with `ib_force_stride = 1` \f$\sum\f$ (column 3)\f$\cdot\Delta t\f$ over the rows approximates the mass the body lost (the gasification rate is the last Runge-Kutta stage's, not stage-weighted like columns 4 and 5); with a larger stride a row holds one step of every `ib_force_stride`, so the sum no longer applies. Every body is summed, so a fixed-`Twall` reacting surface gets its mass loss and a heated wall its heat transfer. `ib_heat_wrt` and `thermal_bc = 3` require `chemistry`, Cartesian coordinates and fixed bodies. A restart reads `Twall` from `restart_data/ib_state` only for `thermal_bc = 3` bodies; state files from before the record grew to 21 fields cannot be restarted.
 
 - `output_partial_domain` activates the output of part of the domain specified by `[x,y,z]_output%%beg` and `[x,y,z]_output%%end`.
 This is useful for large domains where only a portion of the domain is of interest.

src/simulation/m_ibm.fpp

@@ -59,7 +59,7 @@ module m_ibm
     use m_model
     use m_patch_geometries
     use m_collisions
-    use m_riemann_state, only: flux_rsx_vf, flux_src_rsx_vf
+    use m_riemann_state, only: flux_rsx_vf, flux_src_rsx_vf, flux_hatR_rsx_vf
     use m_thermochem, only: num_species, gas_constant, molecular_weights, get_mixture_molecular_weight, get_mixture_energy_mass, &
         & get_mixture_thermal_conductivity_mixavg, get_species_mass_diffusivities_mixavg
 
@@ -2138,6 +2138,23 @@ contains
 
     end subroutine s_prescribed_kinematics
 
+    !> Face flux of component c. The dual-pass HLLD RHS takes the hat_R set at the left face of a cell (hatR) and hat_L at its
+    !! right.
+    function f_ib_face_flux(hatR, j, k, l, c) result(f)
+
+        $:GPU_ROUTINE(function_name='f_ib_face_flux', parallelism='[seq]', cray_inline=True)
+        logical, intent(in) :: hatR
+        integer, intent(in) :: j, k, l, c
+        real(wp)            :: f
+
+        if (hatR) then
+            f = flux_hatR_rsx_vf(j, k, l, c)
+        else
+            f = flux_rsx_vf(j, k, l, c)
+        end if
+
+    end function f_ib_face_flux
+
     !> Add direction idir's fluid-body face fluxes to each IB. These are the fluxes the RHS differences, Riemann (flux_rsx_vf) plus
     !! viscous, capillary, conductive and diffusive (flux_src_rsx_vf), so each body receives exactly what the fluid loses. A face
     !! counts on the rank owning its fluid cell. With IB forces: force (1:3) and torque about the centroid (4:6) of each local IB
@@ -2152,10 +2169,10 @@ contains
         real(wp)                          :: area, w, f_E, f_m
         integer                           :: i, j, k, l, jn, kn, ln, ci, cj, ck, side, ib_id, ib_tmp, xp, yp, zp, g
         integer                           :: j0, k0, l0
-        logical                           :: has_src, has_src_E, do_ft, all_ibs
+        logical                           :: has_src, has_src_E, do_ft, all_ibs, dual_pass, hatR
 
-        if (num_ibs == 0) return
         do_ft = moving_immersed_boundary_flag .or. ib_state_wrt
+        if (num_ibs == 0 .or. .not. (do_ft .or. ib_exchange)) return
         all_ibs = ib_heat_wrt
 
         if (do_ft .and. idir == 1) then
@@ -2171,6 +2188,7 @@ contains
         ! The conditions under which s_compute_additional_physics_rhs differences flux_src for the momentum and the energy
         has_src = viscous .or. surface_tension .or. heat_conduction
         has_src_E = has_src .or. chem_params%diffusion
+        dual_pass = hypo_nc_mode == hypo_nc_mode_dual_pass
         w = 1._wp
         if (time_stepper == time_stepper_rk2) w = 0.5_wp
         if (time_stepper == time_stepper_rk3) w = w_rk3(stage)
@@ -2181,7 +2199,7 @@ contains
 
         ! Face (j, k, l) lies between that cell and the next one along idir
         $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, jn, kn, ln, ci, cj, ck, side, ib_id, ib_tmp, xp, yp, zp, g, area, &
-                            & f, r, tq, f_E, f_m]')
+                            & f, r, tq, f_E, f_m, hatR]')
         do l = l0, p
             do k = k0, n
                 do j = j0, m
@@ -2197,6 +2215,7 @@ contains
                         if (jn <= m .and. kn <= n .and. ln <= p) side = -1
                     end if
                     if (side /= 0) then
+                        hatR = dual_pass .and. side == -1
                         ci = jn; cj = kn; ck = ln
                         if (side == -1) then
                             ci = j; cj = k; ck = l
@@ -2220,7 +2239,7 @@ contains
                                     & wp)*(glb_bounds(3)%end - glb_bounds(3)%beg))
                                 f = 0._wp
                                 do i = 1, num_dims
-                                    f(i) = flux_rsx_vf(j, k, l, eqn_idx%mom%beg + i - 1)
+                                    f(i) = f_ib_face_flux(hatR, j, k, l, eqn_idx%mom%beg + i - 1)
                                     if (has_src) f(i) = f(i) + flux_src_rsx_vf(j, k, l, eqn_idx%mom%beg + i - 1)
                                     f(i) = real(side, wp)*area*f(i)
                                 end do
@@ -2233,11 +2252,11 @@ contains
                                 end do
                             end if
                             if (all_ibs .or. patch_ib(ib_id)%thermal_bc == 3) then
-                                f_E = flux_rsx_vf(j, k, l, eqn_idx%E)
+                                f_E = f_ib_face_flux(hatR, j, k, l, eqn_idx%E)
                                 if (has_src_E) f_E = f_E + flux_src_rsx_vf(j, k, l, eqn_idx%E)
                                 f_m = 0._wp
                                 do i = eqn_idx%cont%beg, eqn_idx%cont%end
-                                    f_m = f_m + flux_rsx_vf(j, k, l, i)
+                                    f_m = f_m + f_ib_face_flux(hatR, j, k, l, i)
                                 end do
                                 g = patch_ib(ib_id)%gbl_patch_id
                                 $:GPU_ATOMIC(atomic='update')

src/simulation/m_rhs.fpp

@@ -804,10 +804,7 @@ contains
 
                     ! Momentum (IB forces) and energy/mass (thermal_bc = 3, ib_heat_wrt) the fluid exchanges with each IB through
                     ! this direction's faces. AMR: coarse faces only; a fine pass would restart the sum on its block.
-                    if (ib .and. (moving_immersed_boundary_flag .or. ib_state_wrt .or. ib_exchange) &
-                        & .and. .not. amr_in_fine_advance) then
-                        call s_accumulate_ib_face_fluxes(id, stage)
-                    end if
+                    if (.not. amr_in_fine_advance) call s_accumulate_ib_face_fluxes(id, stage)
 
                     ! Bubble dynamics source terms
                     if (bubbles_euler) then
@@ -860,6 +857,8 @@ contains
                 call nvtxStartRange("RHS-ADVECTION-SRC")
                 call s_compute_advection_source_term(id, rhs_hatR_vf, q_cons_qp, q_prim_qp, .false.)
                 call nvtxEndRange
+                ! flux_rsx_vf holds the hat_L set and flux_hatR_rsx_vf the hat_R set, each the RHS's at its own faces
+                if (.not. amr_in_fine_advance) call s_accumulate_ib_face_fluxes(id, stage)
             end do
 
             $:GPU_PARALLEL_LOOP(private='[i, j, k, l]', collapse=4)

src/simulation/m_riemann_solvers.fpp

@@ -152,7 +152,10 @@ contains
             @:ALLOCATE(nc_iface_vel_rsx_vf(-1:m_alloc, -1:n_alloc, -1:p_alloc, 1:num_dims))
         end if
 
-        if (hypo_nc_mode == hypo_nc_mode_dual_pass) then
+        if (hypo_nc_mode /= hypo_nc_mode_dual_pass) then
+            ! degenerate placeholder: m_ibm names it in an untaken branch
+            @:ALLOCATE(flux_hatR_rsx_vf(-1:-1, -1:-1, -1:-1, 1:1))
+        else
             @:ALLOCATE(flux_hatR_rsx_vf(-1:m_alloc, -1:n_alloc, -1:p_alloc, 1:sys_size))
             if (use_nc_iface_vel) then
                 @:ALLOCATE(nc_iface_vel_hatR_rsx_vf(-1:m_alloc, -1:n_alloc, -1:p_alloc, 1:num_dims))
@@ -190,8 +193,8 @@ contains
             @:DEALLOCATE(nc_iface_vel_rsx_vf)
         end if
         @:DEALLOCATE(mom_sp_rsx_vf)
+        @:DEALLOCATE(flux_hatR_rsx_vf)
         if (hypo_nc_mode == hypo_nc_mode_dual_pass) then
-            @:DEALLOCATE(flux_hatR_rsx_vf)
             if (use_nc_iface_vel) then
                 @:DEALLOCATE(nc_iface_vel_hatR_rsx_vf)
             end if