2nd order IBM

docs/documentation/case.md

@@ -356,6 +356,8 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi
 | `num_particle_clouds` | Integer | Number of particle bed specifications to generate immersed boundary patches from |
 | `ib_neighborhood_radius`    | Integer | Parameter that controls the neighborhood size for IB detection. |
 | `many_ib_patch_parallelism` | Logical | Parallelize over IB patches instead of grid cells (better for many small patches). |
+| `ib_second_order_vel` | Logical | Extrapolate the ghost-point velocity linearly through the boundary intercept. |
+| `ib_ip_min_dist`      | Real    | Minimum boundary-to-image-point distance, in local cell widths, for `ib_second_order_vel`. |
 | `geometry`           | Integer | Geometry configuration of the patch.|
 | `x[y,z]_centroid`    | Real    | Centroid of the applied geometry in the [x,y,z]-direction. |
 | `length_x[y,z]`      | Real    | Length, if applicable, in the [x,y,z]-direction. |
@@ -450,6 +452,10 @@ Additional details on this specification can be found in [NACA airfoil](https://
 
 - `ib_coefficient_of_friction` is the coefficient of friction used in IB collisions.
 
+- `ib_second_order_vel` replaces the default ghost-point velocity, which imposes the wall velocity at the ghost point, with the linear extrapolation of Mittal et al. (2008) through the boundary intercept: \f$u_{GP} = u_{BI} - \frac{d_{GP}}{d_{IP}}(u_{IP} - u_{BI})\f$. Here \f$d_{GP}\f$ is the ghost-point distance to the wall and \f$d_{IP}\f$ the wall distance to the image point. Slip walls extrapolate only the normal component. Pressure and density keep the zero-gradient mirror.
+
+- `ib_ip_min_dist` sets \f$d_{IP} = \max(d_{GP}, \text{ib\_ip\_min\_dist}\,\Delta)\f$, with \f$\Delta\f$ the smallest local cell width, so image points of near-wall ghost points interpolate from fluid cells. The default 0 is the pure mirror (\f$d_{IP} = d_{GP}\f$). About \f$\sqrt{2}\f$ in 2D (\f$\sqrt{3}\f$ in 3D) keeps the stencil off a locally planar wall.
+
 - `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.
 
 #### Particle Clouds

examples/2D_ibm_second_order_airfoil/case.py

@@ -0,0 +1,119 @@
+#!/usr/bin/env python3
+"""
+NACA 0012 at M = 0.3, alpha = 2 deg, 200 cells per chord: the 2D_ibm_airfoil_surface_pressure example
+switched to a viscous no-slip wall (Re = 1000 per chord), with frequent snapshots for time averaging.
+Env: RE, T_STOP, T_SAVE, SECOND_ORDER=1 (PR branch only: second-order IB velocities).
+"""
+
+import json
+import math
+import os
+
+Ma = 0.3
+alpha_deg = 2.0
+gamma = 1.4
+rho_inf, U_inf, chord = 1.0, 1.0, 1.0
+Re = float(os.environ.get("RE", 1000.0))
+P_inf = rho_inf * U_inf**2 / (gamma * Ma**2)
+
+case = {
+    # --- Output ---
+    "run_time_info": "T",
+    "format": 2,
+    "precision": 2,
+    "parallel_io": "T",
+    "prim_vars_wrt": "T",
+    # --- Domain: pre-stretch extents; the stretching maps them outward ---
+    "x_domain%beg": -3.0,
+    "x_domain%end": 4.0,
+    "y_domain%beg": -3.0,
+    "y_domain%end": 3.0,
+    "m": 1399,
+    "n": 1199,
+    "p": 0,
+    "cyl_coord": "F",
+    "stretch_x": "T",
+    "a_x": 15.0,
+    "x_a": -0.8,
+    "x_b": 1.8,
+    "loops_x": 2,
+    "stretch_y": "T",
+    "a_y": 15.0,
+    "y_a": -0.7,
+    "y_b": 0.7,
+    "loops_y": 2,
+    # --- Time stepping ---
+    "cfl_adap_dt": "T",
+    "cfl_target": 0.5,
+    "n_start": 0,
+    "t_save": 0.2,
+    "t_stop": 6.0,
+    # --- Numerics ---
+    "num_patches": 1,
+    "num_fluids": 1,
+    "model_eqns": 2,
+    "alt_soundspeed": "F",
+    "mpp_lim": "F",
+    "mixture_err": "T",
+    "time_stepper": 3,
+    "weno_order": 5,
+    "weno_eps": 1.0e-10,
+    "weno_Re_flux": "F",
+    "weno_avg": "T",
+    "avg_state": 2,
+    "mapped_weno": "T",
+    "null_weights": "F",
+    "mp_weno": "F",
+    "riemann_solver": 2,
+    "low_Mach": 2,
+    "wave_speeds": 1,
+    "viscous": "T",
+    "fd_order": 4,
+    # --- Uniform freestream ---
+    "patch_icpp(1)%geometry": 3,
+    "patch_icpp(1)%x_centroid": 0.0,
+    "patch_icpp(1)%y_centroid": 0.0,
+    "patch_icpp(1)%length_x": 1.0e3,
+    "patch_icpp(1)%length_y": 1.0e3,
+    "patch_icpp(1)%vel(1)": U_inf,
+    "patch_icpp(1)%vel(2)": 0.0,
+    "patch_icpp(1)%pres": P_inf,
+    "patch_icpp(1)%alpha_rho(1)": rho_inf,
+    "patch_icpp(1)%alpha(1)": 1.0,
+    "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
+    "fluid_pp(1)%eos": "ideal_gas",
+    "fluid_pp(1)%Re(1)": Re,
+    # --- Characteristic far-field boundaries ---
+    "bc_x%beg": -7,
+    "bc_x%grcbc_in": "T",
+    "bc_x%vel_in(1)": U_inf,
+    "bc_x%vel_in(2)": 0.0,
+    "bc_x%pres_in": P_inf,
+    "bc_x%alpha_rho_in(1)": rho_inf,
+    "bc_x%alpha_in(1)": 1.0,
+    "bc_x%end": -8,
+    "bc_x%grcbc_out": "T",
+    "bc_x%pres_out": P_inf,
+    "bc_y%beg": -9,
+    "bc_y%end": -9,
+    # --- IB: NACA 0012 (m must be > 0, so a negligible camber stands in for 0) ---
+    "ib": "T",
+    "num_ibs": 1,
+    "patch_ib(1)%geometry": 4,
+    "patch_ib(1)%x_centroid": 0.0,
+    "patch_ib(1)%y_centroid": 0.0,
+    "patch_ib(1)%airfoil_id": 1,
+    "patch_ib(1)%angles(3)": -math.radians(alpha_deg),
+    "patch_ib(1)%slip": "F",
+    "patch_ib(1)%moving_ibm": 0,
+    "ib_airfoil(1)%c": chord,
+    "ib_airfoil(1)%t": 0.12,
+    "ib_airfoil(1)%p": 0.4,
+    "ib_airfoil(1)%m": 1.0e-9,
+}
+
+# the parameters for second-order IBM 
+case.update({"ib_second_order_vel": "T", "ib_ip_min_dist": 1.5})
+
+
+print(json.dumps(case, indent=4))

src/common/m_derived_types.fpp

@@ -529,6 +529,7 @@ module m_derived_types
     type ghost_point
         integer, dimension(3)        :: loc            !< Physical location of the ghost point
         real(wp), dimension(3)       :: ip_loc         !< Physical location of the image point
+        real(wp)                     :: ip_dist        !< Distance from the boundary intercept to the image point
         integer, dimension(3)        :: ip_grid        !< Top left grid point of IP
         real(wp), dimension(2, 2, 2) :: interp_coeffs  !< Interpolation Coefficients of image point
         logical                      :: interp_valid   !< .false. if every image point stencil cell lies inside an IB

src/simulation/m_global_parameters.fpp

@@ -520,6 +520,8 @@ contains
         ib_force_wrt = .false.
         ib_force_stride = 1
         many_ib_patch_parallelism = .false.
+        ib_second_order_vel = .false.
+        ib_ip_min_dist = 0._wp
 
         ! Bubble modeling (sim-specific)
         bubble_model = 1

src/simulation/m_ibm.fpp

@@ -83,6 +83,7 @@ contains
         @:ACC_SETUP_SFs(corrected_gps)
 
         $:GPU_ENTER_DATA(copyin='[num_gps]')
+        $:GPU_UPDATE(device='[ib_second_order_vel, ib_ip_min_dist]')
 
         if (collision_model > 0) call s_initialize_collisions_module()
 
@@ -195,7 +196,7 @@ contains
         real(wp), intent(in)          :: rho, pres_IP
         real(wp), intent(out)         :: pres_GP
 
-        pres_GP = pres_IP/min(max(1._wp - 2._wp*abs(gp%levelset) &
+        pres_GP = pres_IP/min(max(1._wp - (abs(gp%levelset) + gp%ip_dist) &
                               & *rho/pres_IP*dot_product(patch_ib(gp_patch_id)%force/patch_ib(gp_patch_id)%mass, &
                               & gp%levelset_norm), 5.e-1_wp), 2._wp)
 
@@ -259,6 +260,9 @@ contains
             if (buf > 0._wp) vel_GP = vel_GP + v_blow_eff*norm/buf
         end if
 
+        ! extrapolate linearly through the boundary intercept, Mittal et al. (2008)
+        if (ib_second_order_vel) vel_GP = vel_GP - abs(gp%levelset)/max(gp%ip_dist, sgm_eps)*(vel_IP - vel_GP)
+
     end subroutine s_compute_ghost_point_velocity
 
     !> Update the conservative variables at the ghost points
@@ -708,6 +712,7 @@ contains
     impure subroutine s_compute_image_points()
 
         real(wp)                        :: dist
+        real(wp)                        :: min_cell_width
         real(wp), dimension(3)          :: norm
         real(wp), dimension(3)          :: physical_loc
         real(wp)                        :: temp_loc
@@ -723,8 +728,8 @@ contains
 
         bounds_error = .false.
 
-        $:GPU_PARALLEL_LOOP(private='[q, gp, i, j, k, physical_loc, patch_id, dist, norm, dim, bound, dir, index, temp_loc, &
-                            & s_cc]', copy='[bounds_error]', present='[ghost_points]')
+        $:GPU_PARALLEL_LOOP(private='[q, gp, i, j, k, physical_loc, patch_id, dist, min_cell_width, norm, dim, bound, dir, index, &
+                            & temp_loc, s_cc]', copy='[bounds_error]', present='[ghost_points]')
         do q = 1, num_gps
             gp = ghost_points(q)
             i = gp%loc(1)
@@ -742,7 +747,14 @@ contains
             patch_id = gp%ib_patch_id
             dist = abs(real(gp%levelset, kind=wp))
             norm(:) = gp%levelset_norm
-            ghost_points(q)%ip_loc(:) = physical_loc(:) + 2*dist*norm(:)
+            ghost_points(q)%ip_dist = dist
+            ! keep the image point stencil off the wall
+            if (ib_second_order_vel) then
+                min_cell_width = min(dx(i), dy(j))
+                if (p > 0) min_cell_width = min(min_cell_width, dz(k))
+                ghost_points(q)%ip_dist = max(dist, ib_ip_min_dist*min_cell_width)
+            end if
+            ghost_points(q)%ip_loc(:) = physical_loc(:) + (dist + ghost_points(q)%ip_dist)*norm(:)
 
             ! Find the closest grid point to the image point
             do dim = 1, num_dims