!> @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 states back into the !! state ghost-state vectors and sets state%is_periodic. !! !! Supported BC types (set via bc_left / bc_right in the &initial_condition !! namelist group): !! !! 'dirichlet' — fixed prescribed state (set once by the driver, no update) !! 'inflow' — same as Dirichlet; use for a driven supersonic inlet !! 'outflow' — linear extrapolation: ghost = 2*boundary - penultimate cell !! 'reflecting' — solid wall: density and energy copied, momentum negated !! 'periodic' — index wrapping handled by spatial_discretization !! '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. !! Requires p_ref_left/right and sigma_nrbc. !! 'characteristic' — full LODI characteristic targeting !! (Poinsot & Lele 1992, Sec. 3); decomposes the boundary !! state into f1/f2/f3 characteristic amplitudes and relaxes !! each independently. Requires nrbc_mode, p_ref_*/u_ref_*/ !! rho_ref_*, sigma_nrbc, sigma_nrbc_entropy. !! 'supersonic_inlet' — fully prescribed supersonic inflow (Ma > 1). !! Semantically equivalent to 'inflow'/'dirichlet'; ghost set !! once by the driver and never updated. !! 'subsonic_inlet' — characteristic-based subsonic inflow (0 < Ma < 1). !! Extrapolates the outgoing Riemann invariant from the !! interior and reconstructs the inlet state from stagnation !! conditions (p_stag, rho_stag) in the config. !! See Poinsot & Lele (1992), Sec. 3.3; Anderson (2003), Sec. 11.3. !! 'supersonic_outlet' — all characteristics exit; zero-gradient ghost (q_ghost = q_wall). !! 'subsonic_outlet' — characteristic-based subsonic outflow (0 < Ma < 1). !! Extrapolates the outgoing Riemann invariant and prescribes !! back pressure (p_back) from the config; density recovered via !! isentropic relation. !! 'neumann' — homogeneous Neumann BC: dq/dn = 0; ghost = boundary cell. !! 'neumann_gradient' — prescribed-gradient Neumann: ghost = q_wall + (dq/dn)*dx. !! Gradient specified per-side via neumann_grad_left/right in config. !! !! 'dirichlet', 'inflow', and 'supersonic_inlet' are identical in implementation; they !! differ only in the user's semantic intent. !! !! '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 implicit none private public :: apply_bcs contains ! --------------------------------------------------------------------------- !> Update ghost-state vectors from the current solution and BC types. !! !! Called at the top of spatial_discretization::compute_resid(state). For !! 'dirichlet' and 'inflow' BCs the ghost state is already correct (set once !! by the driver) and no action is taken. For 'outflow', 'reflecting', and !! 'nonreflecting' the ghost state is recomputed every call. !! !! On exit: !! - state%is_periodic is set consistently with state%bc_left / state%bc_right. !! - state%q_left, state%fl_pos, state%fl_neg are updated for the left ghost. !! - state%q_right, state%fr_pos, state%fr_neg are updated for the right ghost. !! !! @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 ! Periodic: index wrapping in spatial_discretization handles all stencil ! boundary treatment; no ghost-state vectors are used. if (trim(state % cfg % bc_left) == bc_periodic) then state % is_periodic = .true. return end if state % is_periodic = .false. ! --- Left boundary --- call update_ghost(state, state % cfg % bc_left, 'L', & & state % ub(:, 1), state % ub(:, 2), & & state % q_left, state % fl_pos, state % fl_neg) ! --- Right boundary --- call update_ghost(state, state % cfg % bc_right, 'R', & & state % ub(:, state % n_pt), state % ub(:, state % n_pt - 1), & & state % q_right, state % fr_pos, state % fr_neg) end subroutine apply_bcs ! --------------------------------------------------------------------------- !> Compute one ghost-state vector and its split fluxes from a boundary cell. !! !! @param[inout] state Solver instance (reads use_fds, flux_split, gam, !! p_ref_left, p_ref_right, sigma_nrbc). !! @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. !! Used by 'outflow' (linear extrapolation) and !! 'nonreflecting' (subsonic left outflow fallback). !! @param[out] q_ghost Ghost conserved state. !! @param[out] f_pos F^+(q_ghost) — used by FVS stencil assembly. !! @param[out] f_neg F^-(q_ghost) — used by FVS stencil assembly. ! --------------------------------------------------------------------------- subroutine update_ghost(state, bc_type, side, q_wall, q_penult, q_ghost, f_pos, f_neg) 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 !< 'L' or 'R' real(wp), intent(in) :: q_wall(3) !< boundary cell real(wp), intent(in) :: q_penult(3) !< second cell from boundary real(wp), intent(out) :: q_ghost(3), f_pos(3), f_neg(3) select case (trim(bc_type)) case (bc_dirichlet, bc_inflow) ! Ghost state was set once by the driver; nothing to update here. ! Split fluxes were likewise pre-computed by the driver. return case (bc_outflow) ! Linear extrapolation (first-order): ghost = 2*boundary - penultimate. ! Reduces the spurious gradient discontinuity at the boundary and damps ! acoustic wave reflections compared with a zero-order copy. ! See LeVeque (2002), Sec. 7.3. q_ghost = 2.0_wp * q_wall - q_penult ! Positivity guard: fall back to zero-gradient if extrapolation yields ! non-physical density or pressure (can occur near strong rarefactions). 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 case (bc_reflecting) ! Solid wall: copy density and total energy, negate the momentum component. ! See LeVeque (2002), Sec. 7.3. q_ghost(1) = q_wall(1) ! rho q_ghost(2) = -q_wall(2) ! -(rho*u) [mirror velocity] q_ghost(3) = q_wall(3) ! E case (bc_nonreflecting) ! 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. ! ! Primitive extraction is inlined here to avoid a circular USE ! dependency between boundary_conditions and euler_physics. block real(wp) :: rho_b, u_b, p_b, c_b, p_ref real(wp) :: p_ghost, rho_ghost, u_ghost ! Characteristic-mode locals. 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 ! ---- Pressure-relaxation path (Thompson 1987) ---- ! At the right subsonic outflow boundary the u-c eigenvalue is ! negative (wave enters from outside); its amplitude is zeroed by ! constructing the ghost from the isentropic relation and the ! acoustic jump dp = rho*c*du. sigma_nrbc in [0,1] blends between ! fully non-reflecting (0) and a soft Dirichlet target pressure (1). if (side == 'R') then p_ref = state % cfg % p_ref_right if (u_b / c_b >= 1.0_wp) then ! Supersonic outflow: all waves leave; zero-gradient is correct. q_ghost = q_wall else if (u_b > 0.0_wp) then ! Subsonic outflow: zero (or soft-target) incoming acoustic wave. 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 ! Inflow at right boundary: leave Dirichlet ghost as-is. return end if else ! side == 'L' p_ref = state % cfg % p_ref_left if (u_b / c_b <= -1.0_wp) then ! Supersonic inflow at left boundary: Dirichlet; no change. return else if (u_b >= 0.0_wp) then ! Subsonic outflow at left: linear extrapolation. q_ghost = 2.0_wp * q_wall - q_penult else ! Subsonic inflow from right at left boundary: zero incoming wave. 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 ! ---- Full characteristic LODI path (Poinsot & Lele 1992, Sec. 3) ---- ! Characteristic amplitudes, linearised about local rho0, c0: ! f1 = p - rho0*c0*u (u-c wave; incoming at right subsonic outflow) ! f2 = rho - p/c0^2 (entropy wave) ! f3 = p + rho0*c0*u (u+c wave; outgoing at right subsonic outflow) ! Inverse: p = (f1+f3)/2, u = (f3-f1)/(2*rho0*c0), rho = f2 + p/c0^2. ! Outgoing amplitudes are extrapolated; incoming ones are relaxed ! toward far-field reference values using sigma_nrbc (acoustic) and ! sigma_nrbc_entropy (entropy wave) independently. 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 ! Subsonic outflow: f1 incoming; f2, f3 outgoing. 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 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 return else if (u_b >= 0.0_wp) then q_ghost = 2.0_wp * q_wall - q_penult else ! Subsonic inflow from right at left boundary: f3 incoming. 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 case (bc_supersonic_inlet) ! Fully prescribed supersonic inflow; all characteristics enter from outside. ! Ghost was set once by the driver (identical to 'inflow'/'dirichlet'). ! See Anderson (2003), Modern Compressible Flow, Sec. 11.3. return case (bc_supersonic_outlet) ! All characteristics exit the domain; zero-gradient ghost. ! See Anderson (2003), Modern Compressible Flow, Sec. 11.3. q_ghost = q_wall case (bc_neumann) ! Homogeneous Neumann BC: dq/dn = 0. Ghost = boundary cell (zeroth-order). ! Imposes zero first-derivative at the boundary; simpler than 'outflow' ! (which imposes zero second-derivative). See LeVeque (2002), Sec. 7.3. q_ghost = q_wall case (bc_neumann_gradient) ! Prescribed-gradient Neumann: q_ghost = q_wall + (dq/dn) * dx. ! Gradient is specified independently per side via config neumann_grad_left/right. if (side == 'L') then q_ghost = q_wall + state % cfg % neumann_grad_left * state % dx else q_ghost = q_wall + state % cfg % neumann_grad_right * state % dx end if case (bc_subsonic_outlet) ! Characteristic-based subsonic outflow. ! Extrapolates the outgoing Riemann invariant from the interior and ! prescribes back pressure from config; density is recovered isentropically. ! ! Right boundary (rightward subsonic flow, 0 < Ma < 1): ! Outgoing: R+ = u + 2c/(gam-1); prescribe p = p_back_right. ! Left boundary (leftward subsonic outflow): ! Outgoing: R- = u - 2c/(gam-1); prescribe p = p_back_left. ! ! See Anderson (2003), Sec. 11.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 ! outgoing R+ at right boundary p_back = state % cfg % p_back_right else r_out = u_b - 2.0_wp * c_b / gm1 ! outgoing R- at left boundary p_back = state % cfg % p_back_left end if ! Isentropic density from back pressure. 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 case (bc_subsonic_inlet) ! Characteristic-based subsonic inflow. ! One outgoing Riemann invariant is extrapolated from the interior; ! the inlet state is reconstructed from stagnation conditions ! (p_stag, rho_stag) using isentropic relations. ! ! Left boundary (rightward inflow, 0 < Ma < 1): ! Outgoing: R- = u - 2c/(gam-1) (left-running acoustic; extrapolated). ! Stagnation enthalpy: h0 = gam*p_stag/((gam-1)*rho_stag). ! Quadratic: (gam+1)*x^2 + 2*R-*x + (R-^2/2 - h0) = 0, x = c/(gam-1). ! c_in = (gam-1)*x; u_in = R- + 2*x (left inlet). ! rho_in = rho_stag*(c_in/c_stag)^(2/(gam-1)); p_in = rho_in*c_in^2/gam. ! Right boundary: same quadratic with R+ = u + 2c/(gam-1); u_in = R+ - 2*x. ! ! See Poinsot & Lele (1992), J. Comput. Phys. 101, 104-129, Sec. 3.3; ! Anderson (2003), Sec. 11.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 ! Outgoing left-running Riemann invariant (exits through left boundary). 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 ! Outgoing right-running Riemann invariant (exits through right boundary). r_inv = u_b + 2.0_wp * c_b / gm1 end if ! Stagnation enthalpy and speed of sound from supplied conditions. h0 = state % cfg % gam * p_stag / (gm1 * rho_stag) c_stag = sqrt(gm1 * h0) ! Quadratic (gam+1)*x^2 + 2*r_inv*x + (r_inv^2/2 - h0) = 0, x = c/(gam-1). ! Discriminant: disc = r_inv^2*(1-(gam+1)/2) + (gam+1)*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 ! positive root: c > 0 c_in = gm1 * x_root if (side == 'L') then u_in = r_inv + 2.0_wp * x_root ! u = R- + 2*x (left inlet) else u_in = r_inv - 2.0_wp * x_root ! u = R+ - 2*x (right inlet) end if ! Isentropic density: rho = rho_stag*(c/c_stag)^(2/(gam-1)) 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 select ! Recompute split fluxes at the ghost state (FVS path only). if (.not. state % use_fds) then call state % flux_split(q_ghost, f_pos, 1, state % cfg % gam) call state % flux_split(q_ghost, f_neg, -1, state % cfg % gam) end if end subroutine update_ghost end module boundary_conditions