diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 33301a629..195696e40 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 90deaf09a..75b5f6b2c 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 3e58206cf..9d2677458 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 f730a4960..75745fb63 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 482d99369..43354cce2 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 f8fe85a9b..f37e28d43 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 6a5434f5b..11e8dc444 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/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 8e13a687b..ab6abb634 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, & diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 1da790a90..7a85eac19 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 ae1ecc028..2236db081 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"})