!> @file boundary_conditions.f90 !> @brief Boundary condition enforcement for the 1D Euler solver. !! !! Provides apply_bcs(state), which must be called at the start of every residual !! evaluation. It reads the current solution array and the BC type strings !! from @p state, then writes the appropriate ghost values directly into the !! halo cells of state%ub. !! !! Phase B rank-awareness: apply_bcs is a no-op on inner ranks. Only rank 0 !! writes the left-boundary halo, and only rank N-1 writes the right-boundary !! halo. Halo exchange (called from compute_resid) fills the interior halo !! cells of inner ranks with neighbour data. At -np 1 a single rank is both !! `my_rank == 0` and `my_rank == n_ranks - 1`, so both edges are written and !! behaviour matches the serial code. !! !! Supported BC types (set via bc_left / bc_right in the &initial_condition !! namelist group): !! !! 'dirichlet' — fixed prescribed state (halo cells preset by the driver, !! no update here) !! 'inflow' — same as Dirichlet; use for a driven supersonic inlet !! 'outflow' — linear extrapolation: ghost = 2*boundary - penultimate cell !! (applied to every halo cell with the same formula) !! 'reflecting' — solid wall: density and energy mirrored, momentum negated !! 'periodic' — at -np > 1 handled by halo_exchange (left/right neighbours !! wrap around for periodic decomp); at -np 1 apply_bcs !! performs the in-process periodic copy directly !! 'nonreflecting' — characteristic (LODI/Thompson) non-reflecting outflow. !! Two algorithm variants are selected by nrbc_mode: !! 'pressure' (default) — isentropic pressure-relaxation; !! zeroes the incoming acoustic wave by constructing the !! ghost via isentropic relation and dp = rho*c*du. !! 'characteristic' — full LODI characteristic targeting !! (Poinsot & Lele 1992, Sec. 3); decomposes the boundary !! state into f1/f2/f3 characteristic amplitudes and relaxes !! each independently. !! 'supersonic_inlet' — fully prescribed supersonic inflow (Ma > 1). !! Halo cells preset by the driver and never updated. !! 'supersonic_outlet' — all characteristics exit; zero-gradient ghost. !! 'subsonic_inlet' — characteristic-based subsonic inflow (0 < Ma < 1). !! 'subsonic_outlet' — characteristic-based subsonic outflow (0 < Ma < 1). !! 'neumann' — homogeneous Neumann BC: ghost = boundary cell. !! 'neumann_gradient' — prescribed-gradient Neumann: ghost = q_wall + (dq/dn)*dx*k. !! !! For the multi-cell halo (halo_width > 1), the "constant ghost" BCs !! (dirichlet, inflow, supersonic_inlet, supersonic_outlet, neumann, the !! nonreflecting and subsonic variants) fill every halo cell with the same !! ghost vector — equivalent to zero-gradient extrapolation past the first !! halo cell. 'outflow' applies the linear extrapolation to every halo cell !! using the same penultimate-cell baseline. 'reflecting' mirrors halo cell !! k about the boundary onto interior cell k (so halo cell 1-k = mirror of !! interior cell k). !! !! 'outflow' uses first-order linear extrapolation (LeVeque 2002, Sec. 7.3). !! A positivity guard falls back to zero-gradient if extrapolation produces a !! non-physical ghost state (negative density or pressure). !! !! 'nonreflecting' follows Thompson (1987) / Poinsot & Lele (1992). At !! supersonic outflow all waves leave and the ghost degenerates to zero-gradient. !! sigma_nrbc in [0,1] blends between fully non-reflecting (0) and a soft !! Dirichlet target (1) for both modes. module boundary_conditions use precision, only: wp use logger, only: log_warn use option_registry, only: bc_periodic, bc_dirichlet, bc_inflow, bc_outflow, & bc_reflecting, bc_nonreflecting, bc_supersonic_inlet, & bc_supersonic_outlet, bc_neumann, bc_neumann_gradient, & bc_subsonic_outlet, bc_subsonic_inlet, & nrbc_mode_pressure, method_fdm, method_fvm implicit none private public :: apply_bcs contains ! --------------------------------------------------------------------------- !> Method dispatcher for boundary-condition enforcement. !! !! Reads the (single) block's discretization method and routes to the matching !! per-method kernel. The public signature is unchanged, so all existing !! callers (spatial_discretization::compute_resid) keep working irrespective !! of the active method. FDM is the default; the FDM path is byte-identical to !! the historical single-path implementation (the dispatcher adds only a !! select-case over a small character token). !! !! @param[inout] state Solver instance state. ! --------------------------------------------------------------------------- subroutine apply_bcs(state) use solver_state, only: solver_state_t implicit none type(solver_state_t), intent(inout) :: state select case (trim(state % blocks(1) % method)) case (method_fdm) call apply_bcs_fdm(state) case (method_fvm) call apply_bcs_fvm_1d(state) case default error stop 'boundary_conditions: unknown block method "'// & trim(state % blocks(1) % method)//'"' end select end subroutine apply_bcs ! --------------------------------------------------------------------------- !> Update halo cells of state%ub from the current interior solution and BC !! types (NODAL / finite-difference path). Rank-aware: only edge ranks act. !! !! Called (via apply_bcs) at the top of !! spatial_discretization::compute_resid(state), after halo_exchange has !! populated the interior-side halos. !! !! On exit: !! - state%is_periodic is set consistently with state%bc_left / state%bc_right. !! - On rank 0: state%ub(:, 1-h : 0) is updated according to bc_left. !! - On rank n_ranks - 1: state%ub(:, n_local+1 : n_local+h) is updated !! according to bc_right. !! - On inner ranks: no writes (halo_exchange already filled the cells). !! !! @param[inout] state Solver instance state. ! --------------------------------------------------------------------------- subroutine apply_bcs_fdm(state) use solver_state, only: solver_state_t implicit none type(solver_state_t), intent(inout) :: state integer :: h, n_local, k h = state % decomp % halo_width n_local = state % decomp % n_local ! Keep is_periodic in sync with the configuration; consumers in ! spatial_discretization still read this flag. state % is_periodic = (trim(state % cfg % bc_left) == bc_periodic) ! Periodic at -np 1: perform the in-process wrap directly because there ! are no neighbours for halo_exchange to talk to. At -np > 1 the ! decomp%left_neighbour/right_neighbour wrap to the opposite rank and ! halo_exchange supplies the data, so apply_bcs is a no-op for periodic. if (state % is_periodic) then if (state % decomp % n_ranks == 1) then ! n_pt = n_cell + 1 with ub(:, 1) == ub(:, n_pt) (same physical cell). ! Wrap such that ub(:, 1-k) = ub(:, n_pt - k) and ! ub(:, n_pt + k) = ub(:, 1 + k). do k = 1, h state % ub(:, 1 - k) = state % ub(:, n_local - k) state % ub(:, n_local + k) = state % ub(:, 1 + k) end do end if return end if if (state % decomp % my_rank == 0) then call apply_left_bc_into_halo(state, h) end if if (state % decomp % my_rank == state % decomp % n_ranks - 1) then call apply_right_bc_into_halo(state, h, n_local) end if ! Inner ranks (0 < my_rank < n_ranks - 1): nothing to do; halo_exchange ! has already filled state%ub's halo cells with neighbour interior data. end subroutine apply_bcs_fdm ! --------------------------------------------------------------------------- !> Update ghost cells of state%ub for the CELL-CENTERED finite-volume path. !! !! Fills the ghost-cell averages ub(:, 1-h : 0) (left) and !! ub(:, n_local+1 : n_local+h) (right) for the full BC token set. !! Rank-aware exactly like apply_bcs_fdm: only edge ranks write physical BCs; !! periodic at -np > 1 is handled by halo_exchange. !! !! Cell-centered vs nodal anchoring — the only two places it matters: !! * periodic at -np 1: there is NO duplicated endpoint node (unlike the !! nodal FDM grid where ub(:,1) == ub(:,n_pt)), so the in-process wrap !! copies whole interior cell averages: ghost(1-k) <- cell(n_local+1-k) !! and ghost(n_local+k) <- cell(k) (the 1D collapse of !! boundary_2d::apply_periodic_halos_2d). !! * reflecting (solid wall): mirrored about the boundary FACE, i.e. !! ghost(1-k) = mirror(cell k) with normal momentum negated (the 1D !! collapse of boundary_2d::fill_x_edge). !! !! Every other BC operates on the first/second interior cell as q_wall / !! q_penult state vectors (anchor-independent), so apply_left_bc_into_halo / !! apply_right_bc_into_halo and compute_boundary_ghost (and the characteristic !! / extrapolation helpers they call) are SHARED with the FDM path unchanged: !! - prescribed (dirichlet/inflow/supersonic_inlet): ghosts preset, no write; !! - zero-gradient / extrapolation (outflow, neumann, supersonic_outlet, !! neumann_gradient): copied/extrapolated from the adjacent interior cell; !! - characteristic (nonreflecting LODI/Thompson, subsonic_inlet/outlet): !! the relaxation acts on the state vector + gamma/sound speed (not on node !! positions), so the math carries over verbatim with the first interior !! cell supplying the wall state. !! The shared reflecting fill in apply_left_bc_into_halo / apply_right_bc_into_halo !! is ALREADY face-anchored (ghost(1-k) = mirror(cell k)), matching the !! cell-centered requirement, so no FVM-specific reflecting code is needed. !! !! @param[inout] state Solver instance state. ! --------------------------------------------------------------------------- subroutine apply_bcs_fvm_1d(state) use solver_state, only: solver_state_t implicit none type(solver_state_t), intent(inout) :: state integer :: h, n_local, k h = state % decomp % halo_width n_local = state % decomp % n_local ! Keep is_periodic in sync with the configuration; consumers in ! spatial_discretization still read this flag. state % is_periodic = (trim(state % cfg % bc_left) == bc_periodic) ! Periodic at -np 1: in-process wrap of cell averages. Cell-centered grids ! have no shared endpoint cell, so (unlike the nodal FDM wrap) the left ! ghosts take the RIGHTMOST interior cells and the right ghosts take the ! LEFTMOST interior cells: ghost(1-k) = cell(n_local+1-k), ! ghost(n_local+k) = cell(k). At -np > 1 the neighbours wrap and ! halo_exchange supplies the data, so this is a no-op for periodic. if (state % is_periodic) then if (state % decomp % n_ranks == 1) then do k = 1, h state % ub(:, 1 - k) = state % ub(:, n_local + 1 - k) state % ub(:, n_local + k) = state % ub(:, k) end do end if return end if ! Physical edges: first interior cell is q_wall, second is q_penult — the ! same indices the FDM helpers use, so they are reused directly (see the ! routine header for why this is anchor-correct for cell-centered ghosts). if (state % decomp % my_rank == 0) then call apply_left_bc_into_halo(state, h) end if if (state % decomp % my_rank == state % decomp % n_ranks - 1) then call apply_right_bc_into_halo(state, h, n_local) end if ! Inner ranks (0 < my_rank < n_ranks - 1): nothing to do; halo_exchange ! has already filled state%ub's halo cells with neighbour interior data. end subroutine apply_bcs_fvm_1d ! --------------------------------------------------------------------------- !> Compute the left-boundary ghost state and broadcast it across the h !! halo cells `state%ub(:, 1-h .. 0)`. ! --------------------------------------------------------------------------- subroutine apply_left_bc_into_halo(state, h) use solver_state, only: solver_state_t type(solver_state_t), intent(inout) :: state integer, intent(in) :: h real(wp) :: q_wall(3), q_penult(3), q_ghost(3) integer :: k logical :: per_cell_pattern, no_write q_wall = state % ub(:, 1) q_penult = state % ub(:, 2) call compute_boundary_ghost(state, trim(state % cfg % bc_left), 'L', & & q_wall, q_penult, q_ghost, per_cell_pattern, no_write) if (no_write) return ! e.g. dirichlet/inflow/supersonic_inlet — halo preset if (per_cell_pattern) then ! 'reflecting': mirror cell k about the boundary => halo(1-k) = mirror(ub(:,k)) do k = 1, h state % ub(1, 1 - k) = state % ub(1, k) state % ub(2, 1 - k) = -state % ub(2, k) state % ub(3, 1 - k) = state % ub(3, k) end do else do k = 1, h state % ub(:, 1 - k) = q_ghost end do end if end subroutine apply_left_bc_into_halo ! --------------------------------------------------------------------------- !> Compute the right-boundary ghost state and broadcast it across the h !! halo cells `state%ub(:, n_local+1 .. n_local+h)`. ! --------------------------------------------------------------------------- subroutine apply_right_bc_into_halo(state, h, n_local) use solver_state, only: solver_state_t type(solver_state_t), intent(inout) :: state integer, intent(in) :: h, n_local real(wp) :: q_wall(3), q_penult(3), q_ghost(3) integer :: k logical :: per_cell_pattern, no_write q_wall = state % ub(:, n_local) q_penult = state % ub(:, n_local - 1) call compute_boundary_ghost(state, trim(state % cfg % bc_right), 'R', & & q_wall, q_penult, q_ghost, per_cell_pattern, no_write) if (no_write) return if (per_cell_pattern) then do k = 1, h state % ub(1, n_local + k) = state % ub(1, n_local - k + 1) state % ub(2, n_local + k) = -state % ub(2, n_local - k + 1) state % ub(3, n_local + k) = state % ub(3, n_local - k + 1) end do else do k = 1, h state % ub(:, n_local + k) = q_ghost end do end if end subroutine apply_right_bc_into_halo ! --------------------------------------------------------------------------- !> Compute one boundary ghost-state vector from a boundary cell. Caller !! decides whether to broadcast the value across all halo cells or to use !! the per-cell mirror pattern (signalled by @p per_cell_pattern). !! !! @param[inout] state Solver instance (reads gam, BC config, etc.) !! @param[in] bc_type Boundary condition type string. !! @param[in] side Which boundary: 'L' (left) or 'R' (right). !! @param[in] q_wall Conserved state at the interior boundary cell. !! @param[in] q_penult Conserved state at the second-from-boundary cell. !! @param[out] q_ghost Computed ghost conserved state (when applicable). !! @param[out] per_cell_pattern .true. for 'reflecting': caller must mirror !! per halo index. !! @param[out] no_write .true. when the halo should not be touched !! (dirichlet/inflow/supersonic_inlet, or !! supersonic-inflow guards in nrbc paths). ! --------------------------------------------------------------------------- subroutine compute_boundary_ghost(state, bc_type, side, q_wall, q_penult, & q_ghost, per_cell_pattern, no_write) use solver_state, only: solver_state_t implicit none type(solver_state_t), intent(inout) :: state character(len=*), intent(in) :: bc_type character(len=1), intent(in) :: side real(wp), intent(in) :: q_wall(3), q_penult(3) real(wp), intent(out) :: q_ghost(3) logical, intent(out) :: per_cell_pattern, no_write per_cell_pattern = .false. no_write = .false. q_ghost = q_wall ! safe default select case (trim(bc_type)) case (bc_dirichlet, bc_inflow, bc_supersonic_inlet) ! Halo state was pre-filled by the driver/IC code; do not overwrite. no_write = .true. return case (bc_outflow) call outflow_ghost(state, q_wall, q_penult, q_ghost) case (bc_reflecting) ! Per-halo-cell mirror about the boundary: caller mirrors ub(:, k). per_cell_pattern = .true. case (bc_nonreflecting) call nonreflecting_ghost(state, side, q_wall, q_penult, q_ghost, no_write) if (no_write) return case (bc_supersonic_outlet) ! Zero-gradient (all characteristics leave). q_ghost = q_wall case (bc_neumann) q_ghost = q_wall case (bc_neumann_gradient) call neumann_gradient_ghost(state, side, q_wall, q_ghost) case (bc_subsonic_outlet) call subsonic_outlet_ghost(state, side, q_wall, q_ghost) case (bc_subsonic_inlet) call subsonic_inlet_ghost(state, side, q_wall, q_ghost) case (bc_periodic) ! Handled at the apply_bcs level (n_ranks == 1 in-process copy, or ! halo_exchange for multi-rank). Should not reach here. no_write = .true. return case default no_write = .true. return end select end subroutine compute_boundary_ghost ! --------------------------------------------------------------------------- !> 'outflow' ghost: first-order linear extrapolation with a positivity guard. !! !! ghost = 2*boundary - penultimate (LeVeque 2002, Sec. 7.3). Falls back to !! zero-gradient (q_wall) when the extrapolation yields a non-physical ghost !! state (negative density or pressure). ! --------------------------------------------------------------------------- subroutine outflow_ghost(state, q_wall, q_penult, q_ghost) use solver_state, only: solver_state_t type(solver_state_t), intent(in) :: state real(wp), intent(in) :: q_wall(3), q_penult(3) real(wp), intent(out) :: q_ghost(3) ! Linear extrapolation (first-order): ghost = 2*boundary - penultimate. ! LeVeque (2002), Sec. 7.3. Positivity guard falls back to zero-gradient. q_ghost = 2.0_wp * q_wall - q_penult block real(wp) :: rho_g, p_g rho_g = q_ghost(1) if (rho_g > 0.0_wp) then p_g = (q_ghost(3) - 0.5_wp * q_ghost(2)**2 / rho_g) * (state % cfg % gam - 1.0_wp) else p_g = -1.0_wp end if if (rho_g <= 0.0_wp .or. p_g <= 0.0_wp) q_ghost = q_wall end block end subroutine outflow_ghost ! --------------------------------------------------------------------------- !> 'nonreflecting' ghost: characteristic (LODI/Thompson) non-reflecting BC. !! !! Two variants are selected by nrbc_mode: 'pressure' (isentropic !! pressure-relaxation) and 'characteristic' (full LODI characteristic !! targeting). Each retains its exact per-side ('L'/'R') and per-Mach-regime !! branches. @p no_write is set .true. for the supersonic-inflow guards that !! must leave the halo untouched. !! !! See Thompson (1987), J. Comput. Phys. 68, 1-24; !! Poinsot & Lele (1992), J. Comput. Phys. 101, 104-129. ! --------------------------------------------------------------------------- subroutine nonreflecting_ghost(state, side, q_wall, q_penult, q_ghost, no_write) use solver_state, only: solver_state_t type(solver_state_t), intent(in) :: state character(len=1), intent(in) :: side real(wp), intent(in) :: q_wall(3), q_penult(3) real(wp), intent(out) :: q_ghost(3) logical, intent(out) :: no_write no_write = .false. q_ghost = q_wall ! safe default ! Non-reflecting BC; two variants selected by nrbc_mode. ! See Thompson (1987), J. Comput. Phys. 68, 1-24; ! Poinsot & Lele (1992), J. Comput. Phys. 101, 104-129. block real(wp) :: rho_b, u_b, p_b, c_b, p_ref real(wp) :: p_ghost, rho_ghost, u_ghost real(wp) :: rho0, c0, f1_b, f2_b, f3_b real(wp) :: f1_ref, f2_ref, f1_g, f2_g, f3_g real(wp) :: u_ref, rho_ref rho_b = q_wall(1) u_b = q_wall(2) / rho_b p_b = (q_wall(3) - 0.5_wp * rho_b * u_b**2) * (state % cfg % gam - 1.0_wp) c_b = sqrt(state % cfg % gam * p_b / rho_b) if (state % cfg % nrbc_mode == nrbc_mode_pressure) then if (side == 'R') then p_ref = state % cfg % p_ref_right if (u_b / c_b >= 1.0_wp) then q_ghost = q_wall else if (u_b > 0.0_wp) then p_ghost = p_b - state % cfg % sigma_nrbc * (p_b - p_ref) rho_ghost = rho_b * (p_ghost / p_b)**(1.0_wp / state % cfg % gam) u_ghost = u_b + (p_b - p_ghost) / (rho_b * c_b) q_ghost(1) = rho_ghost q_ghost(2) = rho_ghost * u_ghost q_ghost(3) = p_ghost / (state % cfg % gam - 1.0_wp) & & + 0.5_wp * rho_ghost * u_ghost**2 else no_write = .true. return end if else ! side == 'L' p_ref = state % cfg % p_ref_left if (u_b / c_b <= -1.0_wp) then no_write = .true. return else if (u_b >= 0.0_wp) then q_ghost = 2.0_wp * q_wall - q_penult else p_ghost = p_b - state % cfg % sigma_nrbc * (p_b - p_ref) rho_ghost = rho_b * (p_ghost / p_b)**(1.0_wp / state % cfg % gam) u_ghost = u_b - (p_b - p_ghost) / (rho_b * c_b) q_ghost(1) = rho_ghost q_ghost(2) = rho_ghost * u_ghost q_ghost(3) = p_ghost / (state % cfg % gam - 1.0_wp) & & + 0.5_wp * rho_ghost * u_ghost**2 end if end if else ! nrbc_mode == 'characteristic': full LODI rho0 = rho_b c0 = c_b f1_b = p_b - rho0 * c0 * u_b f2_b = rho_b - p_b / c0**2 f3_b = p_b + rho0 * c0 * u_b if (side == 'R') then p_ref = state % cfg % p_ref_right u_ref = state % cfg % u_ref_right rho_ref = state % cfg % rho_ref_right if (u_b / c_b >= 1.0_wp) then q_ghost = q_wall else if (u_b > 0.0_wp) then f1_ref = p_ref - rho0 * c0 * u_ref f2_ref = rho_ref - p_ref / c0**2 f3_g = f3_b f2_g = f2_b - state % cfg % sigma_nrbc_entropy * (f2_b - f2_ref) f1_g = f1_b - state % cfg % sigma_nrbc * (f1_b - f1_ref) p_ghost = 0.5_wp * (f1_g + f3_g) u_ghost = (f3_g - f1_g) / (2.0_wp * rho0 * c0) rho_ghost = f2_g + p_ghost / c0**2 if (rho_ghost <= 0.0_wp .or. p_ghost <= 0.0_wp) then q_ghost = q_wall else q_ghost(1) = rho_ghost q_ghost(2) = rho_ghost * u_ghost q_ghost(3) = p_ghost / (state % cfg % gam - 1.0_wp) & & + 0.5_wp * rho_ghost * u_ghost**2 end if else no_write = .true. return end if else ! side == 'L' p_ref = state % cfg % p_ref_left u_ref = state % cfg % u_ref_left rho_ref = state % cfg % rho_ref_left if (u_b / c_b <= -1.0_wp) then no_write = .true. return else if (u_b >= 0.0_wp) then q_ghost = 2.0_wp * q_wall - q_penult else f1_ref = p_ref - rho0 * c0 * u_ref f2_ref = rho_ref - p_ref / c0**2 f1_g = f1_b f2_g = f2_b - state % cfg % sigma_nrbc_entropy * (f2_b - f2_ref) f3_g = f3_b - state % cfg % sigma_nrbc * (f3_b - (p_ref + rho0 * c0 * u_ref)) p_ghost = 0.5_wp * (f1_g + f3_g) u_ghost = (f3_g - f1_g) / (2.0_wp * rho0 * c0) rho_ghost = f2_g + p_ghost / c0**2 if (rho_ghost <= 0.0_wp .or. p_ghost <= 0.0_wp) then q_ghost = q_wall else q_ghost(1) = rho_ghost q_ghost(2) = rho_ghost * u_ghost q_ghost(3) = p_ghost / (state % cfg % gam - 1.0_wp) & & + 0.5_wp * rho_ghost * u_ghost**2 end if end if end if end if end block end subroutine nonreflecting_ghost ! --------------------------------------------------------------------------- !> 'neumann_gradient' ghost: prescribed-gradient Neumann BC. !! !! ghost = q_wall + (dq/dn) * dx, using the per-side prescribed gradient and !! mesh spacing. ! --------------------------------------------------------------------------- subroutine neumann_gradient_ghost(state, side, q_wall, q_ghost) use solver_state, only: solver_state_t type(solver_state_t), intent(in) :: state character(len=1), intent(in) :: side real(wp), intent(in) :: q_wall(3) real(wp), intent(out) :: q_ghost(3) if (side == 'L') then q_ghost = q_wall + state % cfg % neumann_grad_left * state % mesh % h_left else q_ghost = q_wall + state % cfg % neumann_grad_right * state % mesh % h_right end if end subroutine neumann_gradient_ghost ! --------------------------------------------------------------------------- !> 'subsonic_outlet' ghost: characteristic-based subsonic outflow. !! !! The outgoing Riemann invariant is extrapolated from the interior; density !! is recovered isentropically from the prescribed back pressure. Falls back !! to zero-gradient on a non-physical ghost. ! --------------------------------------------------------------------------- subroutine subsonic_outlet_ghost(state, side, q_wall, q_ghost) use solver_state, only: solver_state_t type(solver_state_t), intent(in) :: state character(len=1), intent(in) :: side real(wp), intent(in) :: q_wall(3) real(wp), intent(out) :: q_ghost(3) block real(wp) :: rho_b, u_b, p_b, c_b, gm1 real(wp) :: r_out, p_back, rho_g, p_g, c_g, u_g gm1 = state % cfg % gam - 1.0_wp rho_b = q_wall(1) u_b = q_wall(2) / rho_b p_b = (q_wall(3) - 0.5_wp * rho_b * u_b**2) * gm1 c_b = sqrt(state % cfg % gam * p_b / rho_b) if (side == 'R') then r_out = u_b + 2.0_wp * c_b / gm1 p_back = state % cfg % p_back_right else r_out = u_b - 2.0_wp * c_b / gm1 p_back = state % cfg % p_back_left end if rho_g = rho_b * (p_back / p_b)**(1.0_wp / state % cfg % gam) p_g = p_back c_g = sqrt(state % cfg % gam * p_g / rho_g) if (side == 'R') then u_g = r_out - 2.0_wp * c_g / gm1 else u_g = r_out + 2.0_wp * c_g / gm1 end if if (rho_g <= 0.0_wp .or. p_g <= 0.0_wp) then call log_warn('boundary_conditions: subsonic_outlet: non-physical ghost, ' & & //'falling back to zero-gradient') q_ghost = q_wall else q_ghost(1) = rho_g q_ghost(2) = rho_g * u_g q_ghost(3) = p_g / gm1 + 0.5_wp * rho_g * u_g**2 end if end block end subroutine subsonic_outlet_ghost ! --------------------------------------------------------------------------- !> 'subsonic_inlet' ghost: characteristic-based subsonic inflow. !! !! The incoming state is reconstructed from the prescribed stagnation !! conditions and the outgoing Riemann invariant extrapolated from the !! interior. Falls back to zero-gradient on a negative discriminant or a !! non-physical ghost. ! --------------------------------------------------------------------------- subroutine subsonic_inlet_ghost(state, side, q_wall, q_ghost) use solver_state, only: solver_state_t type(solver_state_t), intent(in) :: state character(len=1), intent(in) :: side real(wp), intent(in) :: q_wall(3) real(wp), intent(out) :: q_ghost(3) block real(wp) :: rho_b, u_b, p_b, c_b, gm1, gp1 real(wp) :: r_inv, h0, c_stag real(wp) :: disc, x_root, c_in, u_in, rho_in, p_in real(wp) :: p_stag, rho_stag gm1 = state % cfg % gam - 1.0_wp gp1 = state % cfg % gam + 1.0_wp rho_b = q_wall(1) u_b = q_wall(2) / rho_b p_b = (q_wall(3) - 0.5_wp * rho_b * u_b**2) * gm1 c_b = sqrt(state % cfg % gam * p_b / rho_b) if (side == 'L') then p_stag = state % cfg % p_stag_left rho_stag = state % cfg % rho_stag_left r_inv = u_b - 2.0_wp * c_b / gm1 else p_stag = state % cfg % p_stag_right rho_stag = state % cfg % rho_stag_right r_inv = u_b + 2.0_wp * c_b / gm1 end if h0 = state % cfg % gam * p_stag / (gm1 * rho_stag) c_stag = sqrt(gm1 * h0) disc = r_inv**2 * (1.0_wp - 0.5_wp * gp1) + gp1 * h0 if (disc < 0.0_wp) then call log_warn('boundary_conditions: subsonic_inlet: negative discriminant, ' & & //'falling back to zero-gradient (check p_stag / rho_stag)') q_ghost = q_wall else x_root = (-r_inv + sqrt(disc)) / gp1 c_in = gm1 * x_root if (side == 'L') then u_in = r_inv + 2.0_wp * x_root else u_in = r_inv - 2.0_wp * x_root end if rho_in = rho_stag * (c_in / c_stag)**(2.0_wp / gm1) p_in = rho_in * c_in**2 / state % cfg % gam if (rho_in <= 0.0_wp .or. p_in <= 0.0_wp .or. c_in <= 0.0_wp) then call log_warn('boundary_conditions: subsonic_inlet: non-physical ghost, ' & & //'falling back to zero-gradient') q_ghost = q_wall else q_ghost(1) = rho_in q_ghost(2) = rho_in * u_in q_ghost(3) = p_in / gm1 + 0.5_wp * rho_in * u_in**2 end if end if end block end subroutine subsonic_inlet_ghost end module boundary_conditions