Skip to content

Commit 0bcfc2d

Browse files
sbryngelsonclaude
andauthored
Ramp Dirichlet inflows as well as GRCBC ones (#1931)
Co-authored-by: Claude <noreply@anthropic.com>
1 parent c3c45f6 commit 0bcfc2d

4 files changed

Lines changed: 41 additions & 13 deletions

File tree

‎docs/documentation/case.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1390,7 +1390,7 @@ The entries labeled "Characteristic." are characteristic boundary conditions bas
13901390

13911391
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.
13921392

1393-
`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.
1393+
`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. The ramp acts on a GRCBC inflow (`grcbc_in`) or on the ghost-cell velocity of a Dirichlet (-17) boundary or boundary patch normal to that direction, such as a nozzle set into a wall.
13941394

13951395
### Patch types {#patch-types}
13961396

‎src/common/m_boundary_primitives.fpp‎

Lines changed: 22 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,11 @@ module m_boundary_primitives
2121
logical :: dirichlet_from_buffers = .false.
2222
$:GPU_DECLARE(create='[dirichlet_from_buffers]')
2323

24-
public :: f_vel_ramp
24+
!> f_vel_ramp per direction, applied to Dirichlet ghost velocities; set by the time stepper when a ramp is active
25+
real(wp), dimension(3) :: bc_vel_ramp = 1._wp
26+
$:GPU_DECLARE(create='[bc_vel_ramp]')
27+
28+
public :: f_vel_ramp, bc_vel_ramp
2529

2630
contains
2731

@@ -921,7 +925,7 @@ contains
921925
if (bc_loc == -1) then ! bc_x%beg
922926
do i = 1, sys_size
923927
do j = 1, buff_size
924-
q_prim_vf(i)%sf(-j, k, l) = bc_buffers(1, 1)%sf(i, k, l)
928+
q_prim_vf(i)%sf(-j, k, l) = bc_buffers(1, 1)%sf(i, k, l)*f_dir_vel(i, 1)
925929
end do
926930
end do
927931
if ((chemistry .or. heat_conduction) .and. present(q_T_sf)) then
@@ -932,7 +936,7 @@ contains
932936
else !< bc_x%end
933937
do i = 1, sys_size
934938
do j = 1, buff_size
935-
q_prim_vf(i)%sf(m + j, k, l) = bc_buffers(1, 2)%sf(i, k, l)
939+
q_prim_vf(i)%sf(m + j, k, l) = bc_buffers(1, 2)%sf(i, k, l)*f_dir_vel(i, 1)
936940
end do
937941
end do
938942
if ((chemistry .or. heat_conduction) .and. present(q_T_sf)) then
@@ -946,7 +950,7 @@ contains
946950
if (bc_loc == -1) then !< bc_y%beg
947951
do i = 1, sys_size
948952
do j = 1, buff_size
949-
q_prim_vf(i)%sf(k, -j, l) = bc_buffers(2, 1)%sf(k, i, l)
953+
q_prim_vf(i)%sf(k, -j, l) = bc_buffers(2, 1)%sf(k, i, l)*f_dir_vel(i, 2)
950954
end do
951955
end do
952956
if ((chemistry .or. heat_conduction) .and. present(q_T_sf)) then
@@ -957,7 +961,7 @@ contains
957961
else !< bc_y%end
958962
do i = 1, sys_size
959963
do j = 1, buff_size
960-
q_prim_vf(i)%sf(k, n + j, l) = bc_buffers(2, 2)%sf(k, i, l)
964+
q_prim_vf(i)%sf(k, n + j, l) = bc_buffers(2, 2)%sf(k, i, l)*f_dir_vel(i, 2)
961965
end do
962966
end do
963967
if ((chemistry .or. heat_conduction) .and. present(q_T_sf)) then
@@ -972,7 +976,7 @@ contains
972976
if (bc_loc == -1) then !< bc_z%beg
973977
do i = 1, sys_size
974978
do j = 1, buff_size
975-
q_prim_vf(i)%sf(k, l, -j) = bc_buffers(3, 1)%sf(k, l, i)
979+
q_prim_vf(i)%sf(k, l, -j) = bc_buffers(3, 1)%sf(k, l, i)*f_dir_vel(i, 3)
976980
end do
977981
end do
978982
if ((chemistry .or. heat_conduction) .and. present(q_T_sf)) then
@@ -983,7 +987,7 @@ contains
983987
else !< bc_z%end
984988
do i = 1, sys_size
985989
do j = 1, buff_size
986-
q_prim_vf(i)%sf(k, l, p + j) = bc_buffers(3, 2)%sf(k, l, i)
990+
q_prim_vf(i)%sf(k, l, p + j) = bc_buffers(3, 2)%sf(k, l, i)*f_dir_vel(i, 3)
987991
end do
988992
end do
989993
if ((chemistry .or. heat_conduction) .and. present(q_T_sf)) then
@@ -997,6 +1001,17 @@ contains
9971001

9981002
end subroutine s_dirichlet
9991003

1004+
!> Ramp factor for variable i of a Dirichlet ghost cell on a dir-normal face: bc_vel_ramp(dir) for velocities, else 1
1005+
pure function f_dir_vel(i, dir) result(f)
1006+
1007+
$:GPU_ROUTINE(parallelism='[seq]')
1008+
integer, intent(in) :: i, dir
1009+
real(wp) :: f
1010+
1011+
f = merge(bc_vel_ramp(dir), 1._wp, i >= eqn_idx%mom%beg .and. i <= eqn_idx%mom%end)
1012+
1013+
end function f_dir_vel
1014+
10001015
!> Extrapolate QBMM bubble pressure and mass-vapor variables into ghost cells by copying boundary values.
10011016
subroutine s_qbmm_extrapolation(bc_dir, bc_loc, k, l, pb_in, mv_in)
10021017

‎src/simulation/m_time_steppers.fpp‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@ module m_time_steppers
2222
use m_collisions, only: collisions_active
2323
use m_mpi_proxy
2424
use m_boundary_common
25+
use m_boundary_primitives, only: f_vel_ramp, bc_vel_ramp
2526
use m_helper
2627
use m_sim_helpers
2728
use m_fftw
@@ -466,6 +467,13 @@ contains
466467
! reads it, not after. Its GPU_DECLARE only creates device storage and never copies the host value, so
467468
! without this the first RHS of a run reads uninitialised memory and later stages read a stale time.
468469
$:GPU_UPDATE(device='[mytime]')
470+
! Dirichlet inflows ramp too; their ghost fill lives in common code, which cannot read mytime
471+
if (any([bc_x%vel_in_ramp, bc_y%vel_in_ramp, bc_z%vel_in_ramp] > 0._wp)) then
472+
bc_vel_ramp = [f_vel_ramp(bc_x%vel_in_ramp, bc_x%vel_in_t0, bc_x%vel_in_frac0, mytime), &
473+
& f_vel_ramp(bc_y%vel_in_ramp, bc_y%vel_in_t0, bc_y%vel_in_frac0, mytime), &
474+
& f_vel_ramp(bc_z%vel_in_ramp, bc_z%vel_in_t0, bc_z%vel_in_frac0, mytime)]
475+
$:GPU_UPDATE(device='[bc_vel_ramp]')
476+
end if
469477
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, &
470478
& t_step, s)
471479

‎toolchain/mfc/case_validator.py‎

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -34,11 +34,12 @@
3434
# See the contributing guide for how to add entries.
3535
PHYSICS_DOCS = {
3636
"check_inflow_ramp": {
37-
"title": "GRCBC Inflow Ramp",
37+
"title": "Inflow Ramp",
3838
"category": "Boundary Conditions",
3939
"math": r"f(t) = f_0 + (1 - f_0)\left[1 + \tanh\left(6 (t - t_0)/\tau - 3\right)\right]/2",
4040
"explanation": "A ramped inflow scales the inflow velocity from a fraction f_0 of its final value to "
41-
"that value over a duration tau. It requires grcbc_in to act on, a non-negative duration, and f_0 in [0, 1].",
41+
"that value over a duration tau. It requires an inflow to act on (grcbc_in, or a Dirichlet boundary or patch), "
42+
"a non-negative duration, and f_0 in [0, 1].",
4243
},
4344
# Thermodynamic Constraints
4445
"check_stiffened_eos": {
@@ -760,10 +761,14 @@ def check_inflow_ramp(self):
760761
ramp = self.get(f"bc_{d}%vel_in_ramp", 0) or 0
761762
frac0 = self.get(f"bc_{d}%vel_in_frac0", 0) or 0
762763
self.prohibit(ramp < 0, f"bc_{d}%vel_in_ramp must be >= 0")
763-
# a ramp needs an inflow to act on
764+
# a ramp needs an inflow to act on: a GRCBC inflow, or a Dirichlet face or boundary patch
765+
grcbc = self.get(f"bc_{d}%grcbc_in", "F") == "T"
766+
dirichlet = any(self.get(f"bc_{d}%{e}", 0) == -17 for e in ("beg", "end")) or any(
767+
self.get(f"patch_bc({i})%type", 0) == -17 and self.get(f"patch_bc({i})%dir", 0) == "xyz".index(d) + 1 for i in range(1, (self.get("num_bc_patches", 0) or 0) + 1)
768+
)
764769
self.prohibit(
765-
ramp > 0 and self.get(f"bc_{d}%grcbc_in", "F") != "T",
766-
f"bc_{d}%vel_in_ramp requires bc_{d}%grcbc_in",
770+
ramp > 0 and not (grcbc or dirichlet),
771+
f"bc_{d}%vel_in_ramp requires bc_{d}%grcbc_in or a Dirichlet (-17) boundary or patch normal to {d}",
767772
)
768773
self.prohibit(not 0 <= frac0 <= 1, f"bc_{d}%vel_in_frac0 must lie in [0, 1]")
769774

0 commit comments

Comments
 (0)