!> @file spatial_discretization_2d.f90 !> @brief Unsplit method-of-lines spatial residual for 2D Euler. !! R(Q) = -( (Fx_{i+1/2,j}-Fx_{i-1/2,j})/dx + (Gy_{i,j+1/2}-Gy_{i,j-1/2})/dy ). !! Per-line component-wise reconstruction (reused via state % reconstruct) + 4-component !! face flux (euler_riemann_2d). Characteristic projection is a later refinement. !! !! Bounds note: state%ub has declared bounds (neq2d, 1-h:nx+h, 1-h:ny+h). Stencil !! gathers are inlined (not delegated to assumed-shape subroutine dummies) so that !! Fortran component-access retains those explicit lower bounds and halo indices like !! -2 map to the correct ghost cells. module spatial_discretization_2d use precision, only: wp use solver_state_2d, only: solver_state_2d_t, neq2d use boundary_2d, only: apply_halos_2d, wall_edges_2d use euler_riemann_2d, only: face_flux_2d, face_flux_2d_normal, wall_flux_2d_normal use euler_physics_2d, only: split_flux_2d, split_contravariant_2d use option_registry, only: method_fdm use mpi_runtime, only: parallel_fatal implicit none private public :: compute_resid_2d contains !> Allocate (or re-allocate on a dimension change) the reconstruction scratch !! arrays shared by the uniform and curvilinear residual paths. Consistent !! stat= handling: the two paths previously diverged (one omitted stat=). subroutine ensure_resid_scratch_2d(state, nx, ny, sw) type(solver_state_2d_t), intent(inout) :: state integer, intent(in) :: nx, ny, sw integer :: astat if (allocated(state % rs_stencil)) then if (size(state % rs_stencil, 2) /= sw) then deallocate (state % rs_stencil, stat=astat) if (astat /= 0) error stop 'spatial_discretization_2d: rs_stencil deallocation failed' end if end if if (.not. allocated(state % rs_stencil)) then allocate (state % rs_stencil(neq2d, sw), stat=astat) if (astat /= 0) error stop 'spatial_discretization_2d: rs_stencil allocation failed' end if if (allocated(state % rs_fx)) then if (size(state % rs_fx, 2) /= nx + 1) then deallocate (state % rs_fx, stat=astat) if (astat /= 0) error stop 'spatial_discretization_2d: rs_fx deallocation failed' end if end if if (.not. allocated(state % rs_fx)) then allocate (state % rs_fx(neq2d, nx + 1), stat=astat) if (astat /= 0) error stop 'spatial_discretization_2d: rs_fx allocation failed' end if if (allocated(state % rs_gy)) then if (size(state % rs_gy, 2) /= ny + 1) then deallocate (state % rs_gy, stat=astat) if (astat /= 0) error stop 'spatial_discretization_2d: rs_gy deallocation failed' end if end if if (.not. allocated(state % rs_gy)) then allocate (state % rs_gy(neq2d, ny + 1), stat=astat) if (astat /= 0) error stop 'spatial_discretization_2d: rs_gy allocation failed' end if end subroutine ensure_resid_scratch_2d !> Zhang-Shu positivity limiter for a 2D Euler face state (neq2d=4). !! Blends the reconstructed face state toward the cell average until !! density rho >= eps and pressure p >= eps. A no-op in smooth regions. pure subroutine limit_pos_2d(qf, qa, gam) real(wp), intent(inout) :: qf(neq2d) !< reconstructed face state real(wp), intent(in) :: qa(neq2d) !< cell-average conserved state real(wp), intent(in) :: gam real(wp), parameter :: eps = 1.0e-13_wp real(wp) :: rho_f, rho_a, p_f, p_a, ke_f, ke_a, theta, denom rho_f = qf(1); rho_a = qa(1) ! Step 1: enforce rho >= eps if (rho_f < eps) then ! Guard denominator: when rho_a == rho_f the blend is degenerate; θ → 1. denom = rho_a - rho_f ! if/else, not merge — merge evaluates both args, so the division would ! still trap under -ffpe-trap when denom == 0. if (abs(denom) > tiny(1.0_wp)) then theta = (rho_a - eps) / denom else theta = 1.0_wp end if theta = min(1.0_wp, max(0.0_wp, theta)) qf = qa + theta * (qf - qa) end if ! Step 2: enforce p >= eps (p = (E - 0.5*(m_x^2+m_y^2)/rho)*(gam-1)) ! Floor rho_f to eps before dividing (covers degenerate rho_a==rho_f<eps). rho_f = max(qf(1), eps) qf(1) = rho_f ke_f = 0.5_wp * (qf(2)**2 + qf(3)**2) / rho_f p_f = (qf(4) - ke_f) * (gam - 1.0_wp) if (p_f < eps) then ke_a = 0.5_wp * (qa(2)**2 + qa(3)**2) / rho_a p_a = (qa(4) - ke_a) * (gam - 1.0_wp) denom = p_a - p_f if (abs(denom) > tiny(1.0_wp)) then theta = (p_a - eps) / denom else theta = 1.0_wp end if theta = min(1.0_wp, max(0.0_wp, theta)) qf = qa + theta * (qf - qa) end if end subroutine limit_pos_2d subroutine compute_resid_2d(state) type(solver_state_2d_t), intent(inout) :: state integer :: i, j, k, nx, ny, sw, so, col, row real(wp) :: qL(neq2d), qR(neq2d) logical :: limit if (.not. associated(state % reconstruct)) & error stop 'compute_resid_2d: reconstruction scheme not initialised (call init_recon_scheme_2d first)' nx = state % nx_local; ny = state % ny_local ! ---- method-aware dispatch ---- ! FDM is uniform-Cartesian only. FDM+FVS routes to its own operator ! (Phase 6, stubbed below); FDM+FDS reuses the shared uniform body that ! follows -- the nodal FDS operator IS the FVM uniform operator. FVM ! (method never == fdm) keeps the existing uniform/curvilinear split. if (trim(state % blocks(1) % method) == method_fdm) then if (.not. state % mesh % uniform) then ! Curvilinear FDM (nodal transformation metrics). FDS routes to the ! contravariant nodal residual; FVS routes to the contravariant ! flux-vector-splitting nodal residual. if (state % use_fds) then call compute_resid_fdm_curv_2d(state) else call compute_resid_fdm_curv_fvs_2d(state) end if return end if if (.not. state % use_fds) then call compute_resid_fdm_fvs_2d(state) ! Phase 6 return end if ! fdm + fds (uniform): fall through to the shared uniform body below. else if (.not. state % mesh % uniform) then call compute_resid_2d_curvilinear(state) return end if sw = state % stencil_width; so = state % stencil_start_offset limit = state % cfg % use_positivity_limiter call ensure_resid_scratch_2d(state, nx, ny, sw) associate (fx => state % rs_fx, gy => state % rs_gy, stencil => state % rs_stencil) call apply_halos_2d(state) state % resid(:, 1:nx, 1:ny) = 0.0_wp ! ---- x sweep: faces i = 1..nx+1 (face i is between cell i-1 and cell i) ---- ! Left-biased stencil for qL at face i: cells i+so, i+so+1, ..., i+so+sw-1 ! (so < 0 for standard schemes, e.g. weno5: so=-3 => i-3..i+1) ! Right-biased stencil for qR at face i: cells i-so-1, i-so-2, ..., i-so-sw ! (mirror, e.g. weno5: so=-3 => i+2, i+1, ..., i-2) ! Inlined to preserve state%ub's explicit bounds (1-h:nx_local+h). do j = 1, ny do i = 1, nx + 1 ! Left-biased stencil: stencil(:,k) = ub(:, i+so + (k-1), j) do k = 1, sw col = i + so + (k - 1) stencil(:, k) = state % ub(:, col, j) end do call state % reconstruct(stencil, qL) ! Right-biased stencil: stencil(:,k) = ub(:, i-so-1 - (k-1), j) do k = 1, sw col = i - so - 1 - (k - 1) stencil(:, k) = state % ub(:, col, j) end do call state % reconstruct(stencil, qR) ! Zhang-Shu positivity: blend toward the left/right cell averages. ! Guards removed: at np>1 face 1 and face nx_local+1 are interior ! global faces; halos fill ub(:,0,j) and ub(:,nx_local+1,j) so both ! cell-average references are always valid. if (limit) call limit_pos_2d(qL, state % ub(:, i - 1, j), state % cfg % gam) if (limit) call limit_pos_2d(qR, state % ub(:, i, j), state % cfg % gam) ! Non-positive state guard (mirrors 1D FDS path and 2D FVS path). ! Check both qL and qR for non-positive density or pressure before ! calling face_flux_2d, so a reconstruction overshoot aborts cleanly ! instead of producing a nonphysical flux silently. block real(wp) :: rho_l, p_l, rho_r, p_r ! Check density BEFORE computing pressure — the pressure formula ! divides by rho, so a non-positive rho must abort first (else the ! division itself would trap under -ffpe-trap). rho_l = qL(1) rho_r = qR(1) if (rho_l <= 0.0_wp .or. rho_r <= 0.0_wp) & call parallel_fatal('compute_resid_2d: non-positive density in FDS ' & //'reconstructed face (x-sweep); enable use_positivity_limiter') p_l = (qL(4) - 0.5_wp * (qL(2)**2 + qL(3)**2) / rho_l) * (state % cfg % gam - 1.0_wp) p_r = (qR(4) - 0.5_wp * (qR(2)**2 + qR(3)**2) / rho_r) * (state % cfg % gam - 1.0_wp) if (p_l <= 0.0_wp .or. p_r <= 0.0_wp) & call parallel_fatal('compute_resid_2d: non-positive pressure in FDS ' & //'reconstructed face (x-sweep); enable use_positivity_limiter') end block call face_flux_2d(qL, qR, 1, state % cfg % flux_scheme, state % cfg % gam, fx(:, i)) end do do i = 1, nx state % resid(:, i, j) = state % resid(:, i, j) & & - (fx(:, i + 1) - fx(:, i)) / state % dx end do end do ! ---- y sweep: faces j = 1..ny+1 (face j is between cell j-1 and cell j) ---- do i = 1, nx do j = 1, ny + 1 ! Left-biased stencil along y do k = 1, sw row = j + so + (k - 1) stencil(:, k) = state % ub(:, i, row) end do call state % reconstruct(stencil, qL) ! Right-biased stencil along y do k = 1, sw row = j - so - 1 - (k - 1) stencil(:, k) = state % ub(:, i, row) end do call state % reconstruct(stencil, qR) ! Zhang-Shu positivity: blend toward the bottom/top cell averages. ! Guards removed: at np>1 face 1 and face ny_local+1 are interior ! global faces; halos fill ub(:,i,0) and ub(:,i,ny_local+1) so both ! cell-average references are always valid. if (limit) call limit_pos_2d(qL, state % ub(:, i, j - 1), state % cfg % gam) if (limit) call limit_pos_2d(qR, state % ub(:, i, j), state % cfg % gam) ! Non-positive state guard (mirrors 1D FDS path and 2D FVS path). block real(wp) :: rho_l, p_l, rho_r, p_r ! Check density BEFORE computing pressure (see x-sweep note). rho_l = qL(1) rho_r = qR(1) if (rho_l <= 0.0_wp .or. rho_r <= 0.0_wp) & call parallel_fatal('compute_resid_2d: non-positive density in FDS ' & //'reconstructed face (y-sweep); enable use_positivity_limiter') p_l = (qL(4) - 0.5_wp * (qL(2)**2 + qL(3)**2) / rho_l) * (state % cfg % gam - 1.0_wp) p_r = (qR(4) - 0.5_wp * (qR(2)**2 + qR(3)**2) / rho_r) * (state % cfg % gam - 1.0_wp) if (p_l <= 0.0_wp .or. p_r <= 0.0_wp) & call parallel_fatal('compute_resid_2d: non-positive pressure in FDS ' & //'reconstructed face (y-sweep); enable use_positivity_limiter') end block call face_flux_2d(qL, qR, 2, state % cfg % flux_scheme, state % cfg % gam, gy(:, j)) end do do j = 1, ny state % resid(:, i, j) = state % resid(:, i, j) & & - (gy(:, j + 1) - gy(:, j)) / state % dy end do end do end associate end subroutine compute_resid_2d !> Curvilinear residual: R = -(1/V)[ (Fhat_{i+1}-Fhat_i) + (Ghat_{j+1}-Ghat_j) ] !! with contravariant face fluxes Fhat = |S_xi| * f_n(qL,qR; n_xi), !! Ghat = |S_eta| * f_n(qL,qR; n_eta). Reconstruction is identical to the !! uniform path (uniform computational stencil); geometry enters via s_xi/s_eta/vol. subroutine compute_resid_2d_curvilinear(state) type(solver_state_2d_t), intent(inout) :: state integer :: i, j, k, nx, ny, sw, so, col, row real(wp) :: qL(neq2d), qR(neq2d), nrm(2), smag, fn(neq2d) logical :: limit, w_xlo, w_xhi, w_ylo, w_yhi nx = state % nx_local; ny = state % ny_local sw = state % stencil_width; so = state % stencil_start_offset limit = state % cfg % use_positivity_limiter call ensure_resid_scratch_2d(state, nx, ny, sw) associate (fx => state % rs_fx, gy => state % rs_gy, stencil => state % rs_stencil, & mesh => state % mesh) call apply_halos_2d(state) ! Physical reflecting-wall edges get a weak pressure-only wall flux at the ! boundary face (zero mass/energy flux, momentum = p*n) instead of a Riemann ! solve on the reflected ghost; see euler_riemann_2d % wall_flux_2d_normal. call wall_edges_2d(state, w_xlo, w_xhi, w_ylo, w_yhi) state % resid(:, 1:nx, 1:ny) = 0.0_wp ! ---- xi sweep: faces i = 1..nx+1 (ξ-face at node column i) ---- do j = 1, ny do i = 1, nx + 1 do k = 1, sw col = i + so + (k - 1) stencil(:, k) = state % ub(:, col, j) end do call state % reconstruct(stencil, qL) do k = 1, sw col = i - so - 1 - (k - 1) stencil(:, k) = state % ub(:, col, j) end do call state % reconstruct(stencil, qR) if (limit) call limit_pos_2d(qL, state % ub(:, i - 1, j), state % cfg % gam) if (limit) call limit_pos_2d(qR, state % ub(:, i, j), state % cfg % gam) smag = sqrt(mesh % s_xi(1, i, j)**2 + mesh % s_xi(2, i, j)**2) ! A degenerate face metric is rank-local data; error stop here would ! strand the peers in the next collective (halo exchange / par_sum) — ! tear down collectively like the rest of this module (audit 2026-07-06 N4). if (smag <= 0.0_wp) call parallel_fatal('compute_resid_2d_curvilinear: degenerate xi-face') nrm = mesh % s_xi(:, i, j) / smag if (w_xlo .and. i == 1) then call wall_flux_2d_normal(qR, nrm, state % cfg % gam, fn) ! interior side = qR else if (w_xhi .and. i == nx + 1) then call wall_flux_2d_normal(qL, nrm, state % cfg % gam, fn) ! interior side = qL else call face_flux_2d_normal(qL, qR, nrm, state % cfg % flux_scheme, state % cfg % gam, fn) end if fx(:, i) = smag * fn end do do i = 1, nx state % resid(:, i, j) = state % resid(:, i, j) & & - (fx(:, i + 1) - fx(:, i)) / mesh % vol(i, j) end do end do ! ---- eta sweep: faces j = 1..ny+1 (η-face at node row j) ---- do i = 1, nx do j = 1, ny + 1 do k = 1, sw row = j + so + (k - 1) stencil(:, k) = state % ub(:, i, row) end do call state % reconstruct(stencil, qL) do k = 1, sw row = j - so - 1 - (k - 1) stencil(:, k) = state % ub(:, i, row) end do call state % reconstruct(stencil, qR) if (limit) call limit_pos_2d(qL, state % ub(:, i, j - 1), state % cfg % gam) if (limit) call limit_pos_2d(qR, state % ub(:, i, j), state % cfg % gam) smag = sqrt(mesh % s_eta(1, i, j)**2 + mesh % s_eta(2, i, j)**2) ! A degenerate face metric is rank-local data; error stop here would ! strand the peers in the next collective (halo exchange / par_sum) — ! tear down collectively like the rest of this module (audit 2026-07-06 N4). if (smag <= 0.0_wp) call parallel_fatal('compute_resid_2d_curvilinear: degenerate eta-face') nrm = mesh % s_eta(:, i, j) / smag if (w_ylo .and. j == 1) then call wall_flux_2d_normal(qR, nrm, state % cfg % gam, fn) ! interior side = qR else if (w_yhi .and. j == ny + 1) then call wall_flux_2d_normal(qL, nrm, state % cfg % gam, fn) ! interior side = qL else call face_flux_2d_normal(qL, qR, nrm, state % cfg % flux_scheme, state % cfg % gam, fn) end if gy(:, j) = smag * fn end do do j = 1, ny state % resid(:, i, j) = state % resid(:, i, j) & & - (gy(:, j + 1) - gy(:, j)) / mesh % vol(i, j) end do end do end associate end subroutine compute_resid_2d_curvilinear !> Curvilinear-FDM (FDS) residual on NODAL transformation metrics: !! R = -(1/J)[ (Fhat_{i+1/2}-Fhat_{i-1/2}) + (Ghat_{j+1/2}-Ghat_{j-1/2}) ] !! with contravariant face fluxes Fhat = |S_xi| * f_n(qL,qR; n_xi), !! Ghat = |S_eta| * f_n(qL,qR; n_eta) and S_xi=mesh%sx_xi, S_eta=mesh%sx_eta, !! J=mesh%jac. Structurally identical to compute_resid_2d_curvilinear (the FVM !! curvilinear residual); only the metric SOURCE differs — FVM face-area vectors !! (s_xi/s_eta/vol) become the nodal transformation metrics here. Reconstruction !! and the contravariant flux machinery are reused verbatim. Freestream/GCL !! preservation (the Phase 1 gate) relies on the consistent metric evaluation in !! compute_metrics_fdm_2d: a constant flow telescopes the face fluxes to ~0. !! !! BCs: apply_halos_2d routes the FDM method to the nodal kernel apply_bcs_fdm_2d !! (duplicate-endpoint-aware periodic wrap, node-mirror reflecting). Reflecting !! walls reuse the curvilinear weak pressure-wall flux with the nodal metric !! normal (Phase 4 refines the curvilinear nodal wall); periodic BCs (the !! Phase 1 gate) need no wall handling. subroutine compute_resid_fdm_curv_2d(state) type(solver_state_2d_t), intent(inout) :: state integer :: i, j, k, nx, ny, sw, so, col, row real(wp) :: qL(neq2d), qR(neq2d), nrm(2), smag, fn(neq2d) logical :: limit, w_xlo, w_xhi, w_ylo, w_yhi nx = state % nx_local; ny = state % ny_local sw = state % stencil_width; so = state % stencil_start_offset limit = state % cfg % use_positivity_limiter call ensure_resid_scratch_2d(state, nx, ny, sw) associate (fx => state % rs_fx, gy => state % rs_gy, stencil => state % rs_stencil, & mesh => state % mesh) call apply_halos_2d(state) call wall_edges_2d(state, w_xlo, w_xhi, w_ylo, w_yhi) state % resid(:, 1:nx, 1:ny) = 0.0_wp ! ---- xi sweep: faces i = 1..nx+1 (xi-face between node i-1 and node i) ---- do j = 1, ny do i = 1, nx + 1 do k = 1, sw col = i + so + (k - 1) stencil(:, k) = state % ub(:, col, j) end do call state % reconstruct(stencil, qL) do k = 1, sw col = i - so - 1 - (k - 1) stencil(:, k) = state % ub(:, col, j) end do call state % reconstruct(stencil, qR) if (limit) call limit_pos_2d(qL, state % ub(:, i - 1, j), state % cfg % gam) if (limit) call limit_pos_2d(qR, state % ub(:, i, j), state % cfg % gam) smag = sqrt(mesh % sx_xi(1, i, j)**2 + mesh % sx_xi(2, i, j)**2) ! A degenerate face metric is rank-local data; error stop here would ! strand the peers in the next collective (halo exchange / par_sum) — ! tear down collectively like the rest of this module (audit 2026-07-06 N4). if (smag <= 0.0_wp) call parallel_fatal('compute_resid_fdm_curv_2d: degenerate xi-face') nrm = mesh % sx_xi(:, i, j) / smag if (w_xlo .and. i == 1) then call wall_flux_2d_normal(qR, nrm, state % cfg % gam, fn) ! interior side = qR else if (w_xhi .and. i == nx + 1) then call wall_flux_2d_normal(qL, nrm, state % cfg % gam, fn) ! interior side = qL else call face_flux_2d_normal(qL, qR, nrm, state % cfg % flux_scheme, state % cfg % gam, fn) end if fx(:, i) = smag * fn end do do i = 1, nx state % resid(:, i, j) = state % resid(:, i, j) & & - (fx(:, i + 1) - fx(:, i)) / mesh % jac(i, j) end do end do ! ---- eta sweep: faces j = 1..ny+1 (eta-face between node j-1 and node j) ---- do i = 1, nx do j = 1, ny + 1 do k = 1, sw row = j + so + (k - 1) stencil(:, k) = state % ub(:, i, row) end do call state % reconstruct(stencil, qL) do k = 1, sw row = j - so - 1 - (k - 1) stencil(:, k) = state % ub(:, i, row) end do call state % reconstruct(stencil, qR) if (limit) call limit_pos_2d(qL, state % ub(:, i, j - 1), state % cfg % gam) if (limit) call limit_pos_2d(qR, state % ub(:, i, j), state % cfg % gam) smag = sqrt(mesh % sx_eta(1, i, j)**2 + mesh % sx_eta(2, i, j)**2) ! A degenerate face metric is rank-local data; error stop here would ! strand the peers in the next collective (halo exchange / par_sum) — ! tear down collectively like the rest of this module (audit 2026-07-06 N4). if (smag <= 0.0_wp) call parallel_fatal('compute_resid_fdm_curv_2d: degenerate eta-face') nrm = mesh % sx_eta(:, i, j) / smag if (w_ylo .and. j == 1) then call wall_flux_2d_normal(qR, nrm, state % cfg % gam, fn) ! interior side = qR else if (w_yhi .and. j == ny + 1) then call wall_flux_2d_normal(qL, nrm, state % cfg % gam, fn) ! interior side = qL else call face_flux_2d_normal(qL, qR, nrm, state % cfg % flux_scheme, state % cfg % gam, fn) end if gy(:, j) = smag * fn end do do j = 1, ny state % resid(:, i, j) = state % resid(:, i, j) & & - (gy(:, j + 1) - gy(:, j)) / mesh % jac(i, j) end do end do end associate end subroutine compute_resid_fdm_curv_2d !> Curvilinear-FDM FLUX-VECTOR-SPLITTING residual on NODAL transformation !! metrics (the FVS counterpart of compute_resid_fdm_curv_2d): !! R = -(1/J)[ (Fhat_{i+1/2}-Fhat_{i-1/2}) + (Ghat_{j+1/2}-Ghat_{j-1/2}) ] !! with the contravariant split-flux reconstruction !! Fhat_{i+1/2} = recon_L[ Fhat^+ ] + recon_R[ Fhat^- ], !! mirroring the uniform FVS residual compute_resid_fdm_fvs_2d but with !! CONTRAVARIANT split fluxes and the 1/J divisor. !! !! FREESTREAM (the headline gate): the metric used to SPLIT each stencil node is !! FROZEN at the FACE metric (mesh%sx_xi(:,i,j) for face i / mesh%sx_eta(:,i,j) !! for face j) — the same GCL-telescoping face metric the FDS branch uses — NOT !! a per-node metric. Then for a constant free-stream the split is identical at !! every stencil node, the (constant-reproducing) reconstruction returns it !! unchanged, and the face flux collapses to Fhat_face = S_face . (F,G). The !! conservative difference therefore telescopes through the discrete GCL identity !! (Delta_xi S_xi + Delta_eta S_eta = 0) to machine zero, exactly as the FDS !! branch does. A naive per-node-metric split would instead leak the !! reconstruction's truncation error into the free-stream residual. !! !! Cost note: because the metric is per-face, the split is recomputed for each !! stencil node at each face (2*sw split_contravariant_2d calls per face) rather !! than precomputed once per node as in the uniform FVS path. This keeps the !! formulation freestream-consistent and high-order; the LF split (no LAPACK) !! used by the accuracy/positivity gates keeps the cost modest. !! !! BCs: the full nodal BC set is applied via apply_halos_2d's nodal kernel !! (apply_bcs_fdm_2d): periodic, outflow, supersonic_inlet, characteristic !! farfield, and the curvilinear slip wall (reflection about the boundary-node !! metric normal). This is the production FVS residual for the curvilinear !! showcases (cylinder/airfoil/blunt body) and the bump/ramp validation. Unlike !! the FDS path it applies no weak wall-flux at the boundary face (the weak !! pressure wall is FDS-only by design); the nodal ghost reflection is its wall !! treatment. No positivity limiter is applied (FVS reconstructs fluxes, not !! states), mirroring compute_resid_fdm_fvs_2d. subroutine compute_resid_fdm_curv_fvs_2d(state) type(solver_state_2d_t), intent(inout) :: state integer :: i, j, k, nx, ny, sw, so, col, row real(wp) :: fpL(neq2d), fmR(neq2d), gam, sxi(2), seta(2), fp_k(neq2d), fm_k(neq2d) nx = state % nx_local; ny = state % ny_local sw = state % stencil_width; so = state % stencil_start_offset gam = state % cfg % gam call ensure_resid_scratch_2d(state, nx, ny, sw) associate (fx => state % rs_fx, gy => state % rs_gy, stencil => state % rs_stencil, & ub => state % ub, mesh => state % mesh) call apply_halos_2d(state) state % resid(:, 1:nx, 1:ny) = 0.0_wp ! ---- xi sweep: faces i = 1..nx+1, face metric S_xi = mesh%sx_xi(:,i,j) ---- do j = 1, ny do i = 1, nx + 1 sxi = mesh % sx_xi(:, i, j) ! Left-biased stencil on Fhat^+ (split each stencil node with the FACE metric). do k = 1, sw col = i + so + (k - 1) call fvs_vacuum_guard_2d(ub(:, col, j), gam, 'fdm-curv-fvs xi-sweep') call split_contravariant_2d(ub(:, col, j), sxi, state % cfg % flux_scheme, gam, fp_k, fm_k) stencil(:, k) = fp_k end do call state % reconstruct(stencil, fpL) ! Right-biased (mirror) stencil on Fhat^-. do k = 1, sw col = i - so - 1 - (k - 1) call fvs_vacuum_guard_2d(ub(:, col, j), gam, 'fdm-curv-fvs xi-sweep') call split_contravariant_2d(ub(:, col, j), sxi, state % cfg % flux_scheme, gam, fp_k, fm_k) stencil(:, k) = fm_k end do call state % reconstruct(stencil, fmR) fx(:, i) = fpL + fmR end do do i = 1, nx state % resid(:, i, j) = state % resid(:, i, j) & & - (fx(:, i + 1) - fx(:, i)) / mesh % jac(i, j) end do end do ! ---- eta sweep: faces j = 1..ny+1, face metric S_eta = mesh%sx_eta(:,i,j) ---- do i = 1, nx do j = 1, ny + 1 seta = mesh % sx_eta(:, i, j) do k = 1, sw row = j + so + (k - 1) call fvs_vacuum_guard_2d(ub(:, i, row), gam, 'fdm-curv-fvs eta-sweep') call split_contravariant_2d(ub(:, i, row), seta, state % cfg % flux_scheme, gam, fp_k, fm_k) stencil(:, k) = fp_k end do call state % reconstruct(stencil, fpL) do k = 1, sw row = j - so - 1 - (k - 1) call fvs_vacuum_guard_2d(ub(:, i, row), gam, 'fdm-curv-fvs eta-sweep') call split_contravariant_2d(ub(:, i, row), seta, state % cfg % flux_scheme, gam, fp_k, fm_k) stencil(:, k) = fm_k end do call state % reconstruct(stencil, fmR) gy(:, j) = fpL + fmR end do do j = 1, ny state % resid(:, i, j) = state % resid(:, i, j) & & - (gy(:, j + 1) - gy(:, j)) / mesh % jac(i, j) end do end do end associate end subroutine compute_resid_fdm_curv_fvs_2d !> Allocate (or re-allocate on a dimension change) the FVS split-flux scratch !! fp_all/fm_all, halo-padded like ub (1-h:nx+h, 1-h:ny+h). Allocated lazily !! from compute_resid_fdm_fvs_2d so non-FVS paths never pay the memory cost. subroutine ensure_fvs_scratch_2d(state, nx, ny, h) type(solver_state_2d_t), intent(inout) :: state integer, intent(in) :: nx, ny, h integer :: astat logical :: need need = .false. if (allocated(state % fp_all)) then if (size(state % fp_all, 2) /= nx + 2 * h .or. & size(state % fp_all, 3) /= ny + 2 * h) then deallocate (state % fp_all, state % fm_all, stat=astat) if (astat /= 0) error stop 'spatial_discretization_2d: fvs scratch deallocation failed' need = .true. end if else need = .true. end if if (need) then allocate (state % fp_all(neq2d, 1 - h:nx + h, 1 - h:ny + h), stat=astat) if (astat /= 0) error stop 'spatial_discretization_2d: fp_all allocation failed' allocate (state % fm_all(neq2d, 1 - h:nx + h, 1 - h:ny + h), stat=astat) if (astat /= 0) error stop 'spatial_discretization_2d: fm_all allocation failed' end if end subroutine ensure_fvs_scratch_2d !> Vacuum/positivity guard for the FVS split-flux precompute. The FVS split !! routines evaluate sqrt(gam*p/rho); a non-positive density or pressure would !! produce NaN that then propagates silently through the reconstructed face !! fluxes and the residual. Detect it here and fail clearly instead. !! Mirrors the 1D FVS guard in spatial_discretization.f90 (block-guarded !! check + parallel_fatal). Operands are guarded with nested ifs because !! Fortran .and. does not short-circuit (rho is a divisor in p). subroutine fvs_vacuum_guard_2d(q, gam, ctx) real(wp), intent(in) :: q(neq2d), gam character(len=*), intent(in) :: ctx real(wp) :: rho_i, p_i rho_i = q(1) if (rho_i <= 0.0_wp) then call parallel_fatal('compute_resid_fdm_fvs_2d: non-positive density in FVS precompute ('//ctx//')') else p_i = (q(4) - 0.5_wp * (q(2)**2 + q(3)**2) / rho_i) * (gam - 1.0_wp) if (p_i <= 0.0_wp) & call parallel_fatal('compute_resid_fdm_fvs_2d: non-positive pressure in FVS precompute ('//ctx//')') end if end subroutine fvs_vacuum_guard_2d !> FDM + flux-vector-splitting residual (Shu-style conservative finite !! difference). Per direction: split the physical flux at every node !! (interior + halo) into F^+/F^-, reconstruct F^+ with the LEFT-biased !! stencil and F^- with the RIGHT-biased (mirror) stencil, sum the two !! reconstructions at each face, then take the conservative difference. !! !! The stencil offsets are IDENTICAL to the uniform FDS sweeps in !! compute_resid_2d -- only the operand differs: FDS reconstructs the !! conserved state (left -> qL, right -> qR); FVS reconstructs the split !! fluxes (left -> F^+, right -> F^-). Matching offsets is what preserves !! the high-order accuracy (Task 6.3 convergence gate). !! !! No positivity limiter is applied: FVS reconstructs fluxes, not states, so !! the face-state Zhang-Shu guard used by the FDS path does not apply here !! (mirrors the 1D FVS path, which also leaves split fluxes unlimited). subroutine compute_resid_fdm_fvs_2d(state) type(solver_state_2d_t), intent(inout) :: state integer :: i, j, k, nx, ny, sw, so, h, col, row real(wp) :: fpL(neq2d), fmR(neq2d), gam nx = state % nx_local; ny = state % ny_local sw = state % stencil_width; so = state % stencil_start_offset h = state % halo_width gam = state % cfg % gam call ensure_resid_scratch_2d(state, nx, ny, sw) call ensure_fvs_scratch_2d(state, nx, ny, h) associate (fx => state % rs_fx, gy => state % rs_gy, stencil => state % rs_stencil, & fp => state % fp_all, fm => state % fm_all, ub => state % ub) call apply_halos_2d(state) state % resid(:, 1:nx, 1:ny) = 0.0_wp ! ---- x sweep: split the x-flux at every node read by the x stencils ---- ! Precompute over i = 1-h..nx+h (all columns the face stencils reach) and ! j = 1..ny (only interior rows are referenced in x). Corner halos are ! never touched, so non-physical corner ghosts cannot poison split_flux_2d. do j = 1, ny do i = 1 - h, nx + h call fvs_vacuum_guard_2d(ub(:, i, j), gam, 'x-sweep') call split_flux_2d(ub(:, i, j), 1, state % cfg % flux_scheme, gam, & & fp(:, i, j), fm(:, i, j)) end do end do do j = 1, ny do i = 1, nx + 1 ! Left-biased stencil on F^+ (same offsets as the FDS qL stencil). do k = 1, sw col = i + so + (k - 1) stencil(:, k) = fp(:, col, j) end do call state % reconstruct(stencil, fpL) ! Right-biased (mirror) stencil on F^- (same offsets as FDS qR). do k = 1, sw col = i - so - 1 - (k - 1) stencil(:, k) = fm(:, col, j) end do call state % reconstruct(stencil, fmR) fx(:, i) = fpL + fmR end do do i = 1, nx state % resid(:, i, j) = state % resid(:, i, j) & & - (fx(:, i + 1) - fx(:, i)) / state % dx end do end do ! ---- y sweep: split the y-flux at every node read by the y stencils ---- ! Precompute over i = 1..nx (interior columns) and j = 1-h..ny+h. do j = 1 - h, ny + h do i = 1, nx call fvs_vacuum_guard_2d(ub(:, i, j), gam, 'y-sweep') call split_flux_2d(ub(:, i, j), 2, state % cfg % flux_scheme, gam, & & fp(:, i, j), fm(:, i, j)) end do end do do i = 1, nx do j = 1, ny + 1 ! Left-biased stencil on F^+ along y. do k = 1, sw row = j + so + (k - 1) stencil(:, k) = fp(:, i, row) end do call state % reconstruct(stencil, fpL) ! Right-biased (mirror) stencil on F^- along y. do k = 1, sw row = j - so - 1 - (k - 1) stencil(:, k) = fm(:, i, row) end do call state % reconstruct(stencil, fmR) gy(:, j) = fpL + fmR end do do j = 1, ny state % resid(:, i, j) = state % resid(:, i, j) & & - (gy(:, j + 1) - gy(:, j)) / state % dy end do end do end associate end subroutine compute_resid_fdm_fvs_2d end module spatial_discretization_2d