!> @file boundary_2d.f90 !> @brief 2D boundary halo filling: doubly-periodic wrap plus per-edge physical !! BCs (reflecting wall / outflow / supersonic inlet), normal-aware on !! curvilinear grids, with an MPI dual-path for inter-rank edges. module boundary_2d use precision, only: wp use solver_state_2d, only: solver_state_2d_t use option_registry, only: bc_periodic, bc_outflow, bc_reflecting, bc_supersonic_inlet, & bc_farfield, method_fdm use euler_physics_2d, only: primitives_2d implicit none private public :: apply_periodic_halos_2d, apply_halos_2d, wall_edges_2d contains !> Report which of the four block edges are PHYSICAL solid (reflecting) walls, !! for the curvilinear residual's weak pressure-wall flux. Mirrors !! apply_halos_2d's serial/MPI dual-path edge selection exactly (single source !! of truth): serial (np=1 or no Cart comm) => all four edges are physical; !! MPI => an edge is physical only if its Cartesian neighbour is MPI_PROC_NULL. !! An edge is a wall iff it is physical AND its cfg % bc_* is 'reflecting'. subroutine wall_edges_2d(state, w_xlo, w_xhi, w_ylo, w_yhi) use mpi_cart_2d, only: cart_is_setup use mpi_runtime, only: mpi_proc_null_value, n_ranks type(solver_state_2d_t), intent(in) :: state logical, intent(out) :: w_xlo, w_xhi, w_ylo, w_yhi logical :: phys_xlo, phys_xhi, phys_ylo, phys_yhi integer :: null_rank if (n_ranks() == 1 .or. .not. cart_is_setup()) then phys_xlo = .true.; phys_xhi = .true.; phys_ylo = .true.; phys_yhi = .true. else null_rank = mpi_proc_null_value() phys_xlo = (state % decomp_2d % x_lo_neighbour == null_rank) phys_xhi = (state % decomp_2d % x_hi_neighbour == null_rank) phys_ylo = (state % decomp_2d % y_lo_neighbour == null_rank) phys_yhi = (state % decomp_2d % y_hi_neighbour == null_rank) end if w_xlo = phys_xlo .and. trim(state % cfg % bc_left) == bc_reflecting w_xhi = phys_xhi .and. trim(state % cfg % bc_right) == bc_reflecting w_ylo = phys_ylo .and. trim(state % cfg % bc_bottom) == bc_reflecting w_yhi = phys_yhi .and. trim(state % cfg % bc_top) == bc_reflecting end subroutine wall_edges_2d !> Fill all four halo bands by periodic wrap of the opposite interior edge. subroutine apply_periodic_halos_2d(state) type(solver_state_2d_t), intent(inout) :: state integer :: i, j, k, h, nx, ny h = state % halo_width; nx = state % nx_local; ny = state % ny_local ! x halos (wrap columns), full y interior do j = 1, ny do k = 1, h state % ub(:, 1 - k, j) = state % ub(:, nx + 1 - k, j) ! left <- right interior state % ub(:, nx + k, j) = state % ub(:, k, j) ! right <- left interior end do end do ! y halos (wrap rows), full x range INCLUDING the x-halos just filled (corner consistency) do k = 1, h do i = 1 - h, nx + h state % ub(:, i, 1 - k) = state % ub(:, i, ny + 1 - k) ! bottom <- top interior state % ub(:, i, ny + k) = state % ub(:, i, k) ! top <- bottom interior end do end do end subroutine apply_periodic_halos_2d !> Fill all four halo bands per cfg % bc_* edge type. !! !! Dual-path: if no MPI Cartesian communicator has been set up (serial unit !! tests), falls back to the Phase-2C pure per-edge fills for all four edges. !! When the cart comm exists (MPI app / multi-rank tests), calls !! exchange_halos_2d to fill from real neighbours and then applies the !! physical BC fill only on edges whose Cart neighbour is MPI_PROC_NULL. !! !! Diagonal corner ghosts are intentionally left unfilled (the !! dimension-by-dimension residual never reads them). subroutine apply_halos_2d(state) use halo_exchange_2d, only: exchange_halos_2d use mpi_cart_2d, only: cart_is_setup use mpi_runtime, only: mpi_proc_null_value, n_ranks type(solver_state_2d_t), intent(inout) :: state integer :: null_rank ! NODAL finite-difference path: route to the node-anchored BC kernel. The ! FDM grid reflects/wraps about the boundary NODE (not the cell face), so the ! cell-centred FVM fills below would be off by half a cell. Keep the FVM body ! byte-identical for the default (FVM) path. if (trim(state % blocks(1) % method) == method_fdm) then call apply_bcs_fdm_2d(state) return end if if (n_ranks() == 1 .or. .not. cart_is_setup()) then ! Serial (single rank, or no MPI cart): pure per-edge physical/periodic fills. call fill_x_edge(state, .true., trim(state % cfg % bc_left)) call fill_x_edge(state, .false., trim(state % cfg % bc_right)) call fill_y_edge(state, .true., trim(state % cfg % bc_bottom)) call fill_y_edge(state, .false., trim(state % cfg % bc_top)) return end if ! MPI mode: exchange with real neighbours, then physical fills at PROC_NULL edges. null_rank = mpi_proc_null_value() call exchange_halos_2d(state, state % decomp_2d) if (state % decomp_2d % x_lo_neighbour == null_rank) & call fill_x_edge(state, .true., trim(state % cfg % bc_left)) if (state % decomp_2d % x_hi_neighbour == null_rank) & call fill_x_edge(state, .false., trim(state % cfg % bc_right)) if (state % decomp_2d % y_lo_neighbour == null_rank) & call fill_y_edge(state, .true., trim(state % cfg % bc_bottom)) if (state % decomp_2d % y_hi_neighbour == null_rank) & call fill_y_edge(state, .false., trim(state % cfg % bc_top)) end subroutine apply_halos_2d !> Fill the four halo bands for the NODAL finite-difference (FDM) path. !! !! Node anchoring: the physical boundary lies ON the boundary node (column 1 / !! nx_local, row 1 / ny_local), not on a cell face, so the ghost fills differ !! from the cell-centred FVM kernel: !! * reflecting: mirror about the boundary NODE — ghost(1-k)=node(1+k) with the !! wall-normal momentum negated (vs. the FVM face mirror ghost(1-k)=cell(k)). !! * periodic: duplicate-endpoint-aware nodal wrap. The nodal mesh has nodes !! ON both boundaries with node 1 == node nx_local the SAME physical point !! (period = (nx_local-1)*dx), so the wrap SKIPS the duplicate endpoint: !! ub(1-k)=ub(nx_local-k), ub(nx_local+k)=ub(1+k). This is the 2D analogue !! of the verified-correct 1D convention (boundary_conditions::apply_bcs_fdm, !! n_pt = n_cell+1 with ub(:,1)==ub(:,n_pt)). It differs from the FVM !! cell-centred wrap at the boundary (FVM cell averages have no duplicate !! endpoint), so the FDM operator equals FVM only in the INTERIOR away from !! the boundary ghosts — which is what the equivalence oracle in !! test_spatial_discretization_fdm_2d now pins. !! * outflow: zeroth-order extrapolation (copy the boundary node). !! * supersonic_inlet: free-stream conserved state. !! !! Dual-path identical to apply_halos_2d: serial (np=1 / no cart comm) fills all !! four physical edges; under MPI the interior partition edges are filled by !! exchange_halos_2d and only the global-boundary edges (Cart neighbour == !! MPI_PROC_NULL) get a physical-BC fill, so no interior edge is double-filled. !! !! NOTE (np>1 periodic): the MPI cartesian wrap in exchange_halos_2d performs !! the neighbour copy across the periodic boundary; verifying it stays !! consistent with the nodal grid under decomposition is deferred to the !! Phase 5 np>1 full-solve task. subroutine apply_bcs_fdm_2d(state) use halo_exchange_2d, only: exchange_halos_2d use mpi_cart_2d, only: cart_is_setup use mpi_runtime, only: mpi_proc_null_value, n_ranks type(solver_state_2d_t), intent(inout) :: state integer :: null_rank if (n_ranks() == 1 .or. .not. cart_is_setup()) then call fill_x_edge_fdm(state, .true., trim(state % cfg % bc_left)) call fill_x_edge_fdm(state, .false., trim(state % cfg % bc_right)) call fill_y_edge_fdm(state, .true., trim(state % cfg % bc_bottom)) call fill_y_edge_fdm(state, .false., trim(state % cfg % bc_top)) return end if null_rank = mpi_proc_null_value() call exchange_halos_2d(state, state % decomp_2d) if (state % decomp_2d % x_lo_neighbour == null_rank) & call fill_x_edge_fdm(state, .true., trim(state % cfg % bc_left)) if (state % decomp_2d % x_hi_neighbour == null_rank) & call fill_x_edge_fdm(state, .false., trim(state % cfg % bc_right)) if (state % decomp_2d % y_lo_neighbour == null_rank) & call fill_y_edge_fdm(state, .true., trim(state % cfg % bc_bottom)) if (state % decomp_2d % y_hi_neighbour == null_rank) & call fill_y_edge_fdm(state, .false., trim(state % cfg % bc_top)) end subroutine apply_bcs_fdm_2d !> Fill a left (is_min=.true.) or right x-edge band over rows j=1..ny, NODAL. subroutine fill_x_edge_fdm(state, is_min, bc) type(solver_state_2d_t), intent(inout) :: state logical, intent(in) :: is_min character(len=*), intent(in) :: bc integer :: j, k, h, nx, gc, src real(wp) :: q_inf(size(state % ub, 1)), q_ff_ref(size(state % ub, 1)) h = state % halo_width; nx = state % nx_local if (bc == bc_supersonic_inlet) q_inf = ref_state(state) if (bc == bc_farfield) q_ff_ref = ref_state(state) do j = 1, state % ny_local do k = 1, h if (is_min) then; gc = 1 - k; else; gc = nx + k; end if select case (bc) case (bc_periodic) ! Duplicate-endpoint-aware nodal wrap (period (nx-1)*dx; node 1 == node ! nx are the SAME physical point). Skip the duplicate endpoint so the ! periodic stencil has no repeated point: ub(1-k)=ub(nx-k), ! ub(nx+k)=ub(1+k). This is the 2D analogue of the verified-correct 1D ! convention (boundary_conditions::apply_bcs_fdm). if (is_min) then; src = nx - k; else; src = 1 + k; end if state % ub(:, gc, j) = state % ub(:, src, j) case (bc_outflow) if (is_min) then; src = 1; else; src = nx; end if state % ub(:, gc, j) = state % ub(:, src, j) case (bc_reflecting) ! Node mirror about the boundary node; wall-normal momentum negated. if (is_min) then; src = 1 + k; else; src = nx - k; end if if (state % mesh % fdm_curvilinear) then ! Curvilinear nodal slip-wall: reflect the interior partner's ! velocity about the boundary-NODE metric normal sx_xi(:,fcol,j) ! (S_xi = (y_eta,-x_eta)). Every ghost layer uses the single ! boundary-node normal, matching the FVM curvilinear branch in ! fill_x_edge (node-anchored here rather than face-anchored). block real(wp) :: nrm(2), smag integer :: fcol if (is_min) then; fcol = 1; else; fcol = nx + 1; end if smag = sqrt(state % mesh % sx_xi(1, fcol, j)**2 + state % mesh % sx_xi(2, fcol, j)**2) if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate (zero-length) x-boundary node metric' nrm = state % mesh % sx_xi(:, fcol, j) / smag state % ub(:, gc, j) = reflect_velocity_2d(state % ub(:, src, j), nrm) end block else state % ub(:, gc, j) = state % ub(:, src, j) state % ub(2, gc, j) = -state % ub(2, src, j) end if case (bc_supersonic_inlet) state % ub(:, gc, j) = q_inf case (bc_farfield) ! Characteristic far-field about the OUTWARD boundary-NODE metric ! normal (sx_xi points +xi; left-edge outward normal is -sx_xi). ! Uniform nodal FDM has no nodal metric -> use the axis normal. block real(wp) :: nrm(2), smag integer :: fcol, isrc if (is_min) then; fcol = 1; isrc = 1; else; fcol = nx + 1; isrc = nx; end if if (state % mesh % fdm_curvilinear) then smag = sqrt(state % mesh % sx_xi(1, fcol, j)**2 + state % mesh % sx_xi(2, fcol, j)**2) if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate x-boundary node metric (farfield)' nrm = state % mesh % sx_xi(:, fcol, j) / smag if (is_min) nrm = -nrm else nrm = 0.0_wp if (is_min) then; nrm(1) = -1.0_wp; else; nrm(1) = 1.0_wp; end if end if state % ub(:, gc, j) = set_farfield_ghost(state % ub(:, isrc, j), q_ff_ref, nrm, state % cfg % gam) end block case default error stop 'boundary_2d: unknown x-edge FDM BC "'//bc//'"' end select end do end do end subroutine fill_x_edge_fdm !> Fill a bottom (is_min=.true.) or top y-edge band over columns i=1..nx, NODAL. subroutine fill_y_edge_fdm(state, is_min, bc) type(solver_state_2d_t), intent(inout) :: state logical, intent(in) :: is_min character(len=*), intent(in) :: bc integer :: i, k, h, ny, gc, src real(wp) :: q_inf(size(state % ub, 1)), q_ff_ref(size(state % ub, 1)) h = state % halo_width; ny = state % ny_local if (bc == bc_supersonic_inlet) q_inf = ref_state(state) if (bc == bc_farfield) q_ff_ref = ref_state(state) do k = 1, h do i = 1, state % nx_local if (is_min) then; gc = 1 - k; else; gc = ny + k; end if select case (bc) case (bc_periodic) ! Duplicate-endpoint-aware nodal wrap (period (ny-1)*dy; row 1 == row ny ! are the SAME physical point). Skip the duplicate endpoint: ! ub(1-k)=ub(ny-k), ub(ny+k)=ub(1+k). Matches the 1D convention. if (is_min) then; src = ny - k; else; src = 1 + k; end if state % ub(:, i, gc) = state % ub(:, i, src) case (bc_outflow) if (is_min) then; src = 1; else; src = ny; end if state % ub(:, i, gc) = state % ub(:, i, src) case (bc_reflecting) ! Node mirror about the boundary node; wall-normal momentum negated. if (is_min) then; src = 1 + k; else; src = ny - k; end if if (state % mesh % fdm_curvilinear) then ! Curvilinear nodal slip-wall: reflect about the boundary-NODE metric ! normal sx_eta(:,i,frow) (S_eta = (-y_xi,x_xi)); see fill_x_edge_fdm. block real(wp) :: nrm(2), smag integer :: frow if (is_min) then; frow = 1; else; frow = ny + 1; end if smag = sqrt(state % mesh % sx_eta(1, i, frow)**2 + state % mesh % sx_eta(2, i, frow)**2) if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate (zero-length) y-boundary node metric' nrm = state % mesh % sx_eta(:, i, frow) / smag state % ub(:, i, gc) = reflect_velocity_2d(state % ub(:, i, src), nrm) end block else state % ub(:, i, gc) = state % ub(:, i, src) state % ub(3, i, gc) = -state % ub(3, i, src) end if case (bc_supersonic_inlet) state % ub(:, i, gc) = q_inf case (bc_farfield) ! Characteristic far-field about the OUTWARD boundary-NODE metric ! normal (sx_eta points +eta; bottom-edge outward normal is -sx_eta). block real(wp) :: nrm(2), smag integer :: frow, jsrc if (is_min) then; frow = 1; jsrc = 1; else; frow = ny + 1; jsrc = ny; end if if (state % mesh % fdm_curvilinear) then smag = sqrt(state % mesh % sx_eta(1, i, frow)**2 + state % mesh % sx_eta(2, i, frow)**2) if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate y-boundary node metric (farfield)' nrm = state % mesh % sx_eta(:, i, frow) / smag if (is_min) nrm = -nrm else nrm = 0.0_wp if (is_min) then; nrm(2) = -1.0_wp; else; nrm(2) = 1.0_wp; end if end if state % ub(:, i, gc) = set_farfield_ghost(state % ub(:, i, jsrc), q_ff_ref, nrm, state % cfg % gam) end block case default error stop 'boundary_2d: unknown y-edge FDM BC "'//bc//'"' end select end do end do end subroutine fill_y_edge_fdm !> Reflect the velocity of a conserved state about unit normal nrm, preserving !! density and total energy (|v| is unchanged). Inviscid no-penetration wall. !! The density is floored to eps (matching limit_pos_2d): the halo fill runs !! BEFORE any residual-side positivity guard can fire, so a transient !! non-positive boundary-adjacent state (near-vacuum off a strong expansion) !! must not divide-by-zero / SIGFPE here (audit 2026-07-06 N2); the residual !! guards still abort cleanly on the interior state itself. pure function reflect_velocity_2d(q, nrm) result(qr) real(wp), intent(in) :: q(:), nrm(2) real(wp) :: qr(size(q)) real(wp), parameter :: eps = 1.0e-13_wp real(wp) :: rho, u, v, vn rho = max(q(1), eps) u = q(2) / rho v = q(3) / rho vn = u * nrm(1) + v * nrm(2) qr(1) = rho qr(2) = rho * (u - 2.0_wp * vn * nrm(1)) qr(3) = rho * (v - 2.0_wp * vn * nrm(2)) qr(4) = q(4) end function reflect_velocity_2d !> Conserved free-stream reference state from config. pure function ref_state(state) result(q) type(solver_state_2d_t), intent(in) :: state real(wp) :: q(size(state % ub, 1)) associate (c => state % cfg) q(1) = c % rho_inf q(2) = c % rho_inf * c % u_inf q(3) = c % rho_inf * c % v_inf q(4) = c % p_inf / (c % gam - 1.0_wp) + 0.5_wp * c % rho_inf * (c % u_inf**2 + c % v_inf**2) end associate end function ref_state !> Characteristic (Riemann-invariant) far-field ghost state for an OUTWARD unit !! normal nrm. q_int is the boundary-adjacent interior conserved state, q_inf !! the free-stream. With interior/free-stream normal velocities Un_i, Un_inf !! and sound speeds c_i, c_inf: !! * supersonic in/out (|Un_i| >= c_i) short-circuit to pure free-stream / !! pure interior extrapolation; !! * otherwise combine the outgoing (interior) invariant !! R+ = Un_i + 2 c_i/(gam-1) with the incoming (free-stream) !! R- = Un_inf - 2 c_inf/(gam-1) to get Un_b = (R+ + R-)/2, !! c_b = (gam-1)/4 (R+ - R-). Entropy and tangential velocity are taken !! from the free-stream on inflow (Un_b <= 0) and from the interior on !! outflow; rho_b, p_b follow from entropy + c_b. !! gam is passed explicitly (no module-level thermodynamic state here). pure function set_farfield_ghost(q_int, q_inf, nrm, gam) result(q_ghost) real(wp), intent(in) :: q_int(:), q_inf(:), nrm(2), gam real(wp) :: q_ghost(size(q_int)) real(wp), parameter :: eps = 1.0e-13_wp real(wp) :: qi(size(q_int)) real(wp) :: gm1, rho_i, u_i, v_i, p_i, c_i, un_i real(wp) :: rho_o, u_o, v_o, p_o, c_o, un_o real(wp) :: rp, rm, un_b, c_b, s_b, ut, vt, rho_b, p_b, ux, uy gm1 = gam - 1.0_wp ! Positivity guard (audit 2026-07-06 N2): the halo fill runs before the ! residual guards, so a transient non-positive interior density/pressure ! must not reach the divide / sqrt below (NaN ghosts, or SIGFPE under ! -ffpe-trap, with no clean abort). Floor to eps like limit_pos_2d; with ! rho_i, p_i > 0 every downstream expression (c_i, s_b, rho_b) stays ! finite. The free-stream reference is config-validated positive. qi = q_int qi(1) = max(qi(1), eps) call primitives_2d(qi, gam, rho_i, u_i, v_i, p_i) p_i = max(p_i, eps) call primitives_2d(q_inf, gam, rho_o, u_o, v_o, p_o) c_i = sqrt(gam * p_i / rho_i) c_o = sqrt(gam * p_o / rho_o) un_i = u_i * nrm(1) + v_i * nrm(2) un_o = u_o * nrm(1) + v_o * nrm(2) if (un_i <= -c_i) then q_ghost = q_inf; return ! supersonic inflow -> pure free-stream else if (un_i >= c_i) then q_ghost = q_int; return ! supersonic outflow -> pure extrapolation end if rp = un_i + 2.0_wp * c_i / gm1 ! outgoing invariant (from interior) rm = un_o - 2.0_wp * c_o / gm1 ! incoming invariant (from free-stream) un_b = 0.5_wp * (rp + rm) c_b = 0.25_wp * gm1 * (rp - rm) if (un_b <= 0.0_wp) then ! inflow: entropy + tangential from free-stream s_b = p_o / rho_o**gam ut = u_o - un_o * nrm(1) vt = v_o - un_o * nrm(2) else ! outflow: entropy + tangential from interior s_b = p_i / rho_i**gam ut = u_i - un_i * nrm(1) vt = v_i - un_i * nrm(2) end if rho_b = (c_b * c_b / (gam * s_b))**(1.0_wp / gm1) p_b = rho_b * c_b * c_b / gam ux = ut + un_b * nrm(1) uy = vt + un_b * nrm(2) q_ghost(1) = rho_b q_ghost(2) = rho_b * ux q_ghost(3) = rho_b * uy q_ghost(4) = p_b / gm1 + 0.5_wp * rho_b * (ux * ux + uy * uy) end function set_farfield_ghost !> Fill a left (is_min=.true.) or right x-edge band over rows j=1..ny. subroutine fill_x_edge(state, is_min, bc) type(solver_state_2d_t), intent(inout) :: state logical, intent(in) :: is_min character(len=*), intent(in) :: bc integer :: j, k, h, nx, gc, src real(wp) :: q_inf(size(state % ub, 1)), q_ff_ref(size(state % ub, 1)) h = state % halo_width; nx = state % nx_local if (bc == bc_supersonic_inlet) q_inf = ref_state(state) if (bc == bc_farfield) q_ff_ref = ref_state(state) do j = 1, state % ny_local do k = 1, h if (is_min) then; gc = 1 - k; else; gc = nx + k; end if select case (bc) case (bc_periodic) if (is_min) then; src = nx + 1 - k; else; src = k; end if state % ub(:, gc, j) = state % ub(:, src, j) case (bc_outflow) if (is_min) then; src = 1; else; src = nx; end if state % ub(:, gc, j) = state % ub(:, src, j) case (bc_reflecting) if (is_min) then; src = k; else; src = nx + 1 - k; end if if (state % mesh % uniform) then state % ub(:, gc, j) = state % ub(:, src, j) state % ub(2, gc, j) = -state % ub(2, src, j) else ! Curvilinear reflecting wall: all halo layers (gc = 1-k) reflect ! their interior partner about the SINGLE boundary-face normal ! s_xi(:,fcol,j). The ghost layers march along the wall-normal ! direction (into the domain), not along the wall, so there is no ! distinct per-layer wall normal to use; the wall face is the ! correct reflection plane for every layer. No-penetration is ! enforced exactly at the boundary face. This single-normal ! treatment is standard practice for curvilinear FV walls. See ! doc/theory-guide.org (curvilinear boundary conditions). block real(wp) :: nrm(2), smag integer :: fcol if (is_min) then; fcol = 1; else; fcol = nx + 1; end if smag = sqrt(state % mesh % s_xi(1, fcol, j)**2 + state % mesh % s_xi(2, fcol, j)**2) if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate (zero-length) x-boundary face' nrm = state % mesh % s_xi(:, fcol, j) / smag state % ub(:, gc, j) = reflect_velocity_2d(state % ub(:, src, j), nrm) end block end if case (bc_supersonic_inlet) state % ub(:, gc, j) = q_inf case (bc_farfield) ! Characteristic far-field about the OUTWARD face normal. s_xi points ! +xi (into the domain at the left face, out at the right), so the ! left-edge outward normal is -s_xi. Uniform FVM has no s_xi -> use ! the axis normal. block real(wp) :: nrm(2), smag integer :: fcol, isrc if (is_min) then; fcol = 1; isrc = 1; else; fcol = nx + 1; isrc = nx; end if if (state % mesh % uniform) then nrm = 0.0_wp if (is_min) then; nrm(1) = -1.0_wp; else; nrm(1) = 1.0_wp; end if else smag = sqrt(state % mesh % s_xi(1, fcol, j)**2 + state % mesh % s_xi(2, fcol, j)**2) if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate x-boundary face (farfield)' nrm = state % mesh % s_xi(:, fcol, j) / smag if (is_min) nrm = -nrm end if state % ub(:, gc, j) = set_farfield_ghost(state % ub(:, isrc, j), q_ff_ref, nrm, state % cfg % gam) end block case default error stop 'boundary_2d: unknown x-edge BC "'//bc//'"' end select end do end do end subroutine fill_x_edge !> Fill a bottom (is_min=.true.) or top y-edge band over columns i=1..nx. subroutine fill_y_edge(state, is_min, bc) type(solver_state_2d_t), intent(inout) :: state logical, intent(in) :: is_min character(len=*), intent(in) :: bc integer :: i, k, h, ny, gc, src real(wp) :: q_inf(size(state % ub, 1)), q_ff_ref(size(state % ub, 1)) h = state % halo_width; ny = state % ny_local if (bc == bc_supersonic_inlet) q_inf = ref_state(state) if (bc == bc_farfield) q_ff_ref = ref_state(state) do k = 1, h do i = 1, state % nx_local if (is_min) then; gc = 1 - k; else; gc = ny + k; end if select case (bc) case (bc_periodic) if (is_min) then; src = ny + 1 - k; else; src = k; end if state % ub(:, i, gc) = state % ub(:, i, src) case (bc_outflow) if (is_min) then; src = 1; else; src = ny; end if state % ub(:, i, gc) = state % ub(:, i, src) case (bc_reflecting) if (is_min) then; src = k; else; src = ny + 1 - k; end if if (state % mesh % uniform) then state % ub(:, i, gc) = state % ub(:, i, src) state % ub(3, i, gc) = -state % ub(3, i, src) else ! Curvilinear reflecting wall: see fill_x_edge — all halo layers ! reflect about the single boundary-face normal s_eta(:,i,frow); ! the wall face is the correct reflection plane for every layer. block real(wp) :: nrm(2), smag integer :: frow if (is_min) then; frow = 1; else; frow = ny + 1; end if smag = sqrt(state % mesh % s_eta(1, i, frow)**2 + state % mesh % s_eta(2, i, frow)**2) if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate (zero-length) y-boundary face' nrm = state % mesh % s_eta(:, i, frow) / smag state % ub(:, i, gc) = reflect_velocity_2d(state % ub(:, i, src), nrm) end block end if case (bc_supersonic_inlet) state % ub(:, i, gc) = q_inf case (bc_farfield) ! Characteristic far-field about the OUTWARD face normal (s_eta points ! +eta; bottom-edge outward normal is -s_eta). Uniform FVM -> axis. block real(wp) :: nrm(2), smag integer :: frow, jsrc if (is_min) then; frow = 1; jsrc = 1; else; frow = ny + 1; jsrc = ny; end if if (state % mesh % uniform) then nrm = 0.0_wp if (is_min) then; nrm(2) = -1.0_wp; else; nrm(2) = 1.0_wp; end if else smag = sqrt(state % mesh % s_eta(1, i, frow)**2 + state % mesh % s_eta(2, i, frow)**2) if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate y-boundary face (farfield)' nrm = state % mesh % s_eta(:, i, frow) / smag if (is_min) nrm = -nrm end if state % ub(:, i, gc) = set_farfield_ghost(state % ub(:, i, jsrc), q_ff_ref, nrm, state % cfg % gam) end block case default error stop 'boundary_2d: unknown y-edge BC "'//bc//'"' end select end do end do end subroutine fill_y_edge end module boundary_2d