Allow a GRCBC inflow to start up smoothly (#1850)

Author: sbryngelson

docs/documentation/case.md

@@ -1328,6 +1328,9 @@ The entries labeled "Characteristic." are characteristic boundary conditions bas
 | `bc_[x,y,z]%%grcbc_out`        | Logical | Enable grcbc for subsonic outflow (pressure)|
 | `bc_[x,y,z]%%grcbc_vel_out`    | Logical | Enable grcbc for subsonic outflow (pressure + normal velocity) |
 | `bc_[x,y,z]%%vel_in`           | Real Array | Inflow velocities in x, y and z directions |
+| `bc_[x,y,z]%%vel_in_ramp`      | Real | Duration of a smooth start-up of the inflow velocity (0 = none) |
+| `bc_[x,y,z]%%vel_in_t0`        | Real | Time at which that ramp begins |
+| `bc_[x,y,z]%%vel_in_frac0`     | Real | Fraction of the final inflow velocity held before the ramp |
 | `bc_[x,y,z]%%vel_out`          | Real Array | Outflow velocities in x, y and z directions |
 | `bc_[x,y,z]%%pres_in`          | Real    | Inflow pressure |
 | `bc_[x,y,z]%%pres_out`         | Real    | Outflow pressure |
@@ -1336,6 +1339,8 @@ The entries labeled "Characteristic." are characteristic boundary conditions bas
 
 This boundary condition can be used for subsonic inflow (`bc_[x,y,z]%[beg,end]` = -7) and subsonic outflow (`bc_[x,y,z]%[beg,end]` = -8) characteristic boundary conditions. These are based on \cite Pirozzoli13. This enables to provide inflow and outflow conditions outside the computational domain.
 
+`bc_[x,y,z]%%vel_in_ramp` starts the inflow smoothly instead of holding it constant, which is what a jet or a tunnel accelerating from rest requires: the start-up is the event of interest, not a transient to be discarded. The inflow velocity is scaled by \f$f(t) = f_0 + (1 - f_0)\left[1 + \tanh\left(6 (t - t_0)/\tau - 3\right)\right]/2\f$, with \f$\tau\f$ = `vel_in_ramp`, \f$t_0\f$ = `vel_in_t0` and \f$f_0\f$ = `vel_in_frac0`, so it leaves \f$f_0\f$ of the final velocity at \f$t_0\f$ and is within half a percent of it at \f$t_0 + \tau\f$. A boundary with `vel_in_ramp = 0` is held constant, as before.
+
 ### Patch types {#patch-types}
 
 | #    | Name                    | Dim.  | Smooth | Description |

src/common/m_boundary_primitives.fpp

@@ -21,8 +21,27 @@ module m_boundary_primitives
     logical :: dirichlet_from_buffers = .false.
     $:GPU_DECLARE(create='[dirichlet_from_buffers]')
 
+    public :: f_vel_ramp
+
 contains
 
+    !> Velocity scaling for a GRCBC inflow that is ramping up: unity unless `bc_[x,y,z]%vel_in_ramp` is set, so an unramped case is
+    !! untouched. Takes the time as an argument so the caller can evaluate it on the device from `mytime`, which the time stepper
+    !! already places there, rather than computing it on the host and copying the result every Runge-Kutta stage.
+    pure function f_vel_ramp(tau, t0, frac0, t) result(f)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+        real(wp), intent(in) :: tau, t0, frac0, t
+        real(wp)             :: f
+
+        if (tau > 0._wp) then
+            f = frac0 + (1._wp - frac0)*0.5_wp*(1._wp + tanh(6._wp*(t - t0)/tau - 3._wp))
+        else
+            f = 1._wp
+        end if
+
+    end function f_vel_ramp
+
     !> Fill ghost cells by copying the nearest boundary cell value along the specified direction.
     subroutine s_ghost_cell_extrapolation(q_prim_vf, bc_dir, bc_loc, k, l, q_T_sf)
 

src/common/m_derived_types.fpp

@@ -112,8 +112,12 @@ module m_derived_types
         real(wp), dimension(3)              :: vel_in, vel_out
         real(wp), dimension(num_fluids_max) :: alpha_rho_in, alpha_in
         logical                             :: grcbc_in, grcbc_out, grcbc_vel_out
-        logical                             :: isothermal_in, isothermal_out
-        real(wp)                            :: Twall_in, Twall_out
+        !> Smooth start-up of a GRCBC inflow: the inflow velocity is scaled by f(t) = vel_in_frac0 + (1 - vel_in_frac0) (1 + tanh(6
+        !! (t - t0)/tau - 3))/2, so it leaves vel_in_frac0 of its final value at t0 and reaches it after vel_in_ramp. Inactive when
+        !! the ramp duration is zero.
+        real(wp) :: vel_in_ramp, vel_in_t0, vel_in_frac0
+        logical  :: isothermal_in, isothermal_out
+        real(wp) :: Twall_in, Twall_out
     end type int_bounds_info
 
     !> Groups the x, y, z boundary condition begin/end codes for passing as a single argument.

src/simulation/m_cbc.fpp

@@ -12,6 +12,7 @@ module m_cbc
     use m_global_parameters
     use m_variables_conversion
     use m_compute_cbc
+    use m_boundary_primitives, only: f_vel_ramp
     use m_constants, only: riemann_solver_hll, model_eqns_gamma_law, recon_type_weno, recon_type_muscl
     use m_thermochem, only: get_mixture_energy_mass, get_mixture_specific_heat_cv_mass, get_mixture_specific_heat_cp_mass, &
         & gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights, get_species_specific_heats_r, &
@@ -473,6 +474,7 @@ contains
         real(wp)                                               :: dpi_inf_dt
         real(wp)                                               :: dqv_dt
         real(wp)                                               :: dpres_ds
+        real(wp)                                               :: ramp  !< inflow ramp factor; unity unless a ramp is set
 
         #:if USING_AMD
             real(wp), dimension(${AMD_SYS_SIZE_MAX}$) :: L
@@ -596,9 +598,13 @@ contains
                                     & dalpha_rho_ds, dpres_ds, dvel_dt, dadv_dt, dalpha_rho_dt, L, lambda, Ys, dYs_dt, dYs_ds, &
                                     & h_k, Cp_i, Gamma_i, Xs, drho_dt, dpres_dt, dpi_inf_dt, dqv_dt, dgamma_dt, rho, pres, E, &
                                     & gamma, pi_inf, qv, c, Ma, T, sum_Enthalpies, Cv, Cp, e_mix, Mw, R_gas, vel_K_sum, &
-                                    & vel_dv_dt_sum, i, j]', copyin='[dir_idx]')
+                                    & vel_dv_dt_sum, i, j, ramp]', copyin='[dir_idx]')
                 do r = is3%beg, is3%end
                     do k = is2%beg, is2%end
+                        ! Ramp factor for a smoothly starting inflow, evaluated here from mytime rather than
+                        ! computed on the host and copied every Runge-Kutta stage. Unity unless a ramp is set.
+                        ramp = f_vel_ramp(bc_${XYZ}$%vel_in_ramp, bc_${XYZ}$%vel_in_t0, bc_${XYZ}$%vel_in_frac0, mytime)
+
                         ! Transferring the Primitive Variables
                         $:GPU_LOOP(parallelism='[seq]')
                         do i = 1, eqn_idx%cont%end
@@ -735,10 +741,10 @@ contains
                                       & ${CBC_DIR}$))/Del_in(${CBC_DIR}$) - c*Ma*(pres - pres_in(${CBC_DIR}$))/Del_in(${CBC_DIR}$)
                                 end do
                                 if (n > 0) then
-                                    L(eqn_idx%mom%beg + 1) = c*Ma*(vel(dir_idx(2)) - vel_in(${CBC_DIR}$, &
+                                    L(eqn_idx%mom%beg + 1) = c*Ma*(vel(dir_idx(2)) - ramp*vel_in(${CBC_DIR}$, &
                                       & dir_idx(2)))/Del_in(${CBC_DIR}$)
                                     if (p > 0) then
-                                        L(eqn_idx%mom%beg + 2) = c*Ma*(vel(dir_idx(3)) - vel_in(${CBC_DIR}$, &
+                                        L(eqn_idx%mom%beg + 2) = c*Ma*(vel(dir_idx(3)) - ramp*vel_in(${CBC_DIR}$, &
                                           & dir_idx(3)))/Del_in(${CBC_DIR}$)
                                     end if
                                 end if
@@ -747,7 +753,7 @@ contains
                                     L(i) = c*Ma*(adv_local(i + 1 - eqn_idx%E) - alpha_in(i + 1 - eqn_idx%E, &
                                       & ${CBC_DIR}$))/Del_in(${CBC_DIR}$)
                                 end do
-                                L(eqn_idx%adv%end) = rho*c**2._wp*(1._wp + Ma)*(vel(dir_idx(1)) + vel_in(${CBC_DIR}$, &
+                                L(eqn_idx%adv%end) = rho*c**2._wp*(1._wp + Ma)*(vel(dir_idx(1)) + ramp*vel_in(${CBC_DIR}$, &
                                   & dir_idx(1))*sign(1, &
                                   & cbc_loc))/Del_in(${CBC_DIR}$) + c*(1._wp + Ma)*(pres - pres_in(${CBC_DIR}$))/Del_in(${CBC_DIR}$)
                             end if

src/simulation/m_global_parameters.fpp

@@ -121,6 +121,9 @@ module m_global_parameters
     $:GPU_DECLARE(create='[bc_x%vb1, bc_x%vb2, bc_x%vb3, bc_x%ve1, bc_x%ve2, bc_x%ve3]')
     $:GPU_DECLARE(create='[bc_y%vb1, bc_y%vb2, bc_y%vb3, bc_y%ve1, bc_y%ve2, bc_y%ve3]')
     $:GPU_DECLARE(create='[bc_z%vb1, bc_z%vb2, bc_z%vb3, bc_z%ve1, bc_z%ve2, bc_z%ve3]')
+    $:GPU_DECLARE(create='[bc_x%vel_in_ramp, bc_x%vel_in_t0, bc_x%vel_in_frac0]')
+    $:GPU_DECLARE(create='[bc_y%vel_in_ramp, bc_y%vel_in_t0, bc_y%vel_in_frac0]')
+    $:GPU_DECLARE(create='[bc_z%vel_in_ramp, bc_z%vel_in_t0, bc_z%vel_in_frac0]')
     $:GPU_DECLARE(create='[ib_bc_x%beg, ib_bc_x%end, ib_bc_y%beg, ib_bc_y%end, ib_bc_z%beg, ib_bc_z%end]')
 #elif defined(MFC_OpenMP)
     $:GPU_DECLARE(create='[bc_x, bc_y, bc_z]')
@@ -606,6 +609,9 @@ contains
             bc_${dir}$%grcbc_in = .false.
             bc_${dir}$%grcbc_out = .false.
             bc_${dir}$%grcbc_vel_out = .false.
+            bc_${dir}$%vel_in_ramp = 0._wp
+            bc_${dir}$%vel_in_t0 = 0._wp
+            bc_${dir}$%vel_in_frac0 = 0._wp
         #:endfor
 
         ! Lagrangian subgrid bubble model

src/simulation/m_mpi_proxy.fpp

@@ -135,7 +135,10 @@ contains
             & 'bc_x%pres_in','bc_x%pres_out','bc_y%pres_in','bc_y%pres_out',               &
             & 'bc_z%pres_in','bc_z%pres_out',                                               &
             & 'bc_x%Twall_in', 'bc_x%Twall_out', 'bc_y%Twall_in', 'bc_y%Twall_out',       &
-            & 'bc_z%Twall_in', 'bc_z%Twall_out']
+            & 'bc_z%Twall_in', 'bc_z%Twall_out',                                           &
+            & 'bc_x%vel_in_ramp', 'bc_x%vel_in_t0', 'bc_x%vel_in_frac0',                   &
+            & 'bc_y%vel_in_ramp', 'bc_y%vel_in_t0', 'bc_y%vel_in_frac0',                   &
+            & 'bc_z%vel_in_ramp', 'bc_z%vel_in_t0', 'bc_z%vel_in_frac0']
             call MPI_BCAST(${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 

src/simulation/m_start_up.fpp

@@ -1097,6 +1097,10 @@ contains
         $:GPU_UPDATE(device='[bc_y%grcbc_in, bc_y%grcbc_out, bc_y%grcbc_vel_out]')
         $:GPU_UPDATE(device='[bc_z%grcbc_in, bc_z%grcbc_out, bc_z%grcbc_vel_out]')
 
+        $:GPU_UPDATE(device='[bc_x%vel_in_ramp, bc_x%vel_in_t0, bc_x%vel_in_frac0]')
+        $:GPU_UPDATE(device='[bc_y%vel_in_ramp, bc_y%vel_in_t0, bc_y%vel_in_frac0]')
+        $:GPU_UPDATE(device='[bc_z%vel_in_ramp, bc_z%vel_in_t0, bc_z%vel_in_frac0]')
+
         $:GPU_UPDATE(device='[bc_x%isothermal_in, bc_x%isothermal_out]')
         $:GPU_UPDATE(device='[bc_y%isothermal_in, bc_y%isothermal_out]')
         $:GPU_UPDATE(device='[bc_z%isothermal_in, bc_z%isothermal_out]')

src/simulation/m_time_steppers.fpp

@@ -460,6 +460,10 @@ contains
 
         do s = 1, nstage
             call system_clock(stage_t0)
+            ! mytime is read on the device by the GRCBC inflow ramp, so it has to be current before the RHS that
+            ! reads it, not after. Its GPU_DECLARE only creates device storage and never copies the host value, so
+            ! without this the first RHS of a run reads uninitialised memory and later stages read a stale time.
+            $:GPU_UPDATE(device='[mytime]')
             call s_compute_rhs(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, bc_type, rhs_vf, pb_ts(1)%sf, rhs_pb, mv_ts(1)%sf, rhs_mv, &
                                & t_step, s)
 
@@ -532,7 +536,6 @@ contains
                 $:END_GPU_PARALLEL_LOOP()
             end if
 
-            $:GPU_UPDATE(device='[mytime]')
             if (bodyForces) call s_apply_bodyforces(q_cons_ts(1)%vf, q_prim_vf, rhs_vf, rk_coef(s, 3)*dt/rk_coef(s, 4))
 
             if (synthetic_turbulence) call s_apply_synthetic_turbulence_force(q_cons_ts(1)%vf, q_prim_vf, rhs_vf, rk_coef(s, &

toolchain/mfc/case_validator.py

@@ -32,6 +32,13 @@
 # to auto-generate docs/documentation/physics_constraints.md.
 # See the contributing guide for how to add entries.
 PHYSICS_DOCS = {
+    "check_inflow_ramp": {
+        "title": "GRCBC Inflow Ramp",
+        "category": "Boundary Conditions",
+        "math": r"f(t) = f_0 + (1 - f_0)\left[1 + \tanh\left(6 (t - t_0)/\tau - 3\right)\right]/2",
+        "explanation": "A ramped inflow scales the inflow velocity from a fraction f_0 of its final value to "
+        "that value over a duration tau. It requires grcbc_in to act on, a non-negative duration, and f_0 in [0, 1].",
+    },
     # Thermodynamic Constraints
     "check_stiffened_eos": {
         "title": "Stiffened EOS Positivity",
@@ -734,6 +741,19 @@ def check_phase_change(self):
         self.prohibit(ptgalpha_eps is not None and ptgalpha_eps <= 0, "ptgalpha_eps must be positive")
         self.prohibit(ptgalpha_eps is not None and ptgalpha_eps >= 1, "ptgalpha_eps must be less than 1")
 
+    def check_inflow_ramp(self):
+        """Checks constraints on the smooth start-up of a GRCBC inflow"""
+        for d in ("x", "y", "z"):
+            ramp = self.get(f"bc_{d}%vel_in_ramp", 0) or 0
+            frac0 = self.get(f"bc_{d}%vel_in_frac0", 0) or 0
+            self.prohibit(ramp < 0, f"bc_{d}%vel_in_ramp must be >= 0")
+            # a ramp needs an inflow to act on
+            self.prohibit(
+                ramp > 0 and self.get(f"bc_{d}%grcbc_in", "F") != "T",
+                f"bc_{d}%vel_in_ramp requires bc_{d}%grcbc_in",
+            )
+            self.prohibit(not 0 <= frac0 <= 1, f"bc_{d}%vel_in_frac0 must lie in [0, 1]")
+
     def check_ibm(self):
         """Checks constraints on Immersed Boundaries parameters"""
         ib = self.get("ib", "F") == "T"
@@ -2812,6 +2832,7 @@ def validate_common(self):
         self.check_hypoelasticity()
         self.check_phase_change()
         self.check_ibm()
+        self.check_inflow_ramp()
         self.check_eos_selector()
         self.check_stiffened_eos()
         self.check_eos_parameter_sanity()

toolchain/mfc/params/definitions.py

@@ -69,6 +69,9 @@ def _fc(name: str, default: int) -> int:
 
 HINTS = {
     "bc": {
+        "vel_in_ramp": "Duration of the smooth start-up of the inflow velocity (0 = no ramp)",
+        "vel_in_t0": "Time at which the inflow velocity ramp begins",
+        "vel_in_frac0": "Fraction of the final inflow velocity held before the ramp",
         "grcbc_in": "Enables GRCBC subsonic inflow (bc type -7)",
         "grcbc_out": "Enables GRCBC subsonic outflow (bc type -8)",
         "grcbc_vel_out": "GRCBC velocity outlet (requires `grcbc_out`)",
@@ -1110,7 +1113,7 @@ def _load():
     # Extended BC
     for d in ["x", "y", "z"]:
         px = f"bc_{d}%"
-        for a in ["vb1", "vb2", "vb3", "ve1", "ve2", "ve3", "pres_in", "pres_out"]:
+        for a in ["vb1", "vb2", "vb3", "ve1", "ve2", "ve3", "pres_in", "pres_out", "vel_in_ramp", "vel_in_t0", "vel_in_frac0"]:
             _r(f"{px}{a}", REAL, {"bc"})
         for a in ["grcbc_in", "grcbc_out", "grcbc_vel_out"]:
             _r(f"{px}{a}", LOG, {"bc"})