From 597cf75f14c50b74f19e35410dffa0290590d765 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 12 Sep 2026 00:53:29 -0400 Subject: [PATCH 1/2] Keep the GRCBC inflow ramp on the device and restore what the branch had dropped The ramp is evaluated in the CBC kernel from mytime, which the time stepper already places on the device, so nothing is copied per Runge-Kutta stage. Three regressions the branch had introduced are undone: the GRCBC device update in s_initialize_cbc_module (which every existing GRCBC GPU case depends on), its matching deallocate, and the missing device residency for the three new bc_* members the CBC kernel reads. The Dirichlet path is dropped from this branch -- it referenced mytime from src/common, which pre_process and post_process also compile, so it could not build. --- docs/documentation/case.md | 5 +++++ src/common/m_boundary_primitives.fpp | 19 +++++++++++++++++++ src/common/m_derived_types.fpp | 8 ++++++-- src/simulation/m_cbc.fpp | 14 ++++++++++---- src/simulation/m_global_parameters.fpp | 6 ++++++ src/simulation/m_mpi_proxy.fpp | 5 ++++- src/simulation/m_start_up.fpp | 4 ++++ toolchain/mfc/case_validator.py | 21 +++++++++++++++++++++ toolchain/mfc/params/definitions.py | 5 ++++- 9 files changed, 79 insertions(+), 8 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 33301a6296..195696e407 100644 --- a/docs/documentation/case.md +++ b/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 | diff --git a/src/common/m_boundary_primitives.fpp b/src/common/m_boundary_primitives.fpp index 90deaf09a7..75b5f6b2c2 100644 --- a/src/common/m_boundary_primitives.fpp +++ b/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) diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index 3e58206cf7..9d26774588 100644 --- a/src/common/m_derived_types.fpp +++ b/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. diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index f730a49604..75745fb635 100644 --- a/src/simulation/m_cbc.fpp +++ b/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(20) :: 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 diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 482d993690..43354cce20 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/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 diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index f8fe85a9b7..f37e28d43d 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/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 diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 6a5434f5b4..11e8dc4440 100644 --- a/src/simulation/m_start_up.fpp +++ b/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]') diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 1da790a90e..7a85eac199 100644 --- a/toolchain/mfc/case_validator.py +++ b/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() diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index ae1ecc028a..2236db081c 100644 --- a/toolchain/mfc/params/definitions.py +++ b/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"}) From a9de108070abc941b1a06f8e3021e711c77fe44e Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 12 Sep 2026 11:36:28 -0500 Subject: [PATCH 2/2] Push mytime to the device before the RHS that reads it --- src/simulation/m_time_steppers.fpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 8e13a687bb..ab6abb6349 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/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, &