!> @file spatial_discretization.f90 !> @brief Spatial residual computation via flux splitting and reconstruction. module spatial_discretization use precision, only: wp use option_registry, only: hybrid_sensor_jameson, hybrid_sensor_density_gradient, & hybrid_sensor_weno_beta, method_fdm, method_fvm implicit none private public :: compute_resid contains ! --------------------------------------------------------------------------- !> Compute the spatial residual R(Q), dispatching on the block's method. !! !! Thin seam over the per-method residual kernels. The discretization method !! is data on the solver state (`state % blocks(1) % method`); this routine !! selects the matching kernel via an internal `select case` rather than a !! state-typed procedure pointer (which would create a circular module !! dependency). The public signature is unchanged, so all existing callers !! keep working. The 'fdm' arm runs the finite-difference kernel; the 'fvm' !! arm runs the cell-centered finite-volume kernel. !! !! @param[inout] state Solver instance state. ! --------------------------------------------------------------------------- subroutine compute_resid(state) use solver_state, only: solver_state_t type(solver_state_t), intent(inout) :: state select case (trim(state % blocks(1) % method)) case (method_fdm) call compute_resid_fdm(state) case (method_fvm) call compute_resid_fvm_1d(state) case default error stop 'spatial_discretization: unknown block method "'// & trim(state % blocks(1) % method)//'"' end select end subroutine compute_resid ! --------------------------------------------------------------------------- !> Compute the spatial residual R(Q) = -1/dx * (F_{i+1/2} - F_{i-1/2}). !! !! Two code paths are supported, selected by state%use_fds: !! !! **FVS path** (Lax-Friedrichs, Steger-Warming, van Leer): !! 1. Computes F^+(Q_j) and F^-(Q_j) at all grid points AND halo cells. !! 2. Averages adjacent states to get K and K^{-1} at the interface. !! 3. Assembles the stencil for both F^+ (left-biased) and F^- (right-biased) !! directly from the halo'd fp/fm arrays. !! 4. Calls the active reconstruction scheme for each split flux. !! 5. Sums the contributions: F_hat = F_hat^+ + F_hat^-. !! !! **FDS path** (AUSM+, HLL, HLLC, Roe): !! 1. Reconstructs Q_L and Q_R at each face from the halo'd ub array. !! 2. Calls the active FDS solver: F_hat = fds_solver(Q_L, Q_R, gam). !! !! Boundary values are pulled from the halo cells of state%ub (and state%fp, !! state%fm for FVS); apply_bcs writes those at edge ranks and halo_exchange !! fills the inner-rank halos with neighbour interior data. !! !! See course notes, Sec. "Finite Difference WENO Scheme", Eq. (28)-(30). !! !! @param[inout] state Solver instance state. ! --------------------------------------------------------------------------- subroutine compute_resid_fdm(state) use solver_state, only: solver_state_t, neq use euler_physics, only: compute_eigensystem use reconstruction, only: max_stencil_width use boundary_conditions, only: apply_bcs use halo_exchange, only: exchange_halos use positivity_limiter, only: limit_positivity use mpi_runtime, only: parallel_fatal use timer, only: timer_start, timer_stop type(solver_state_t), intent(inout) :: state integer :: i, iface, h, n_local real(wp) :: flux_pos(neq), flux_neg(neq) real(wp) :: r_mat(neq, neq), r_inv(neq, neq), q_avg(neq) real(wp) :: f_stencil(neq, max_stencil_width) real(wp) :: q_stencil(neq, max_stencil_width) real(wp) :: q_face_L(neq), q_face_R(neq) integer :: stencil_width real(wp) :: sensor_val logical :: is_smooth stencil_width = state % stencil_width h = state % decomp % halo_width n_local = state % n_pt if (state % cfg % do_timing) call timer_start(state % perf % resid) ! Update halo cells of state%ub from the current solution and BC types. ! Phase B: apply_bcs writes into halo cells of state%ub directly; Task 9 ! prepends a call to halo_exchange to fill inner-rank halos with neighbour ! interior data before apply_bcs runs on the edge ranks. call exchange_halos(state, state % decomp) call apply_bcs(state) if (state % use_fds) then ! ----------------------------------------------------------------------- ! FDS path (AUSM+, HLL, HLLC, Roe): reconstruct Q_L and Q_R at each face. ! ----------------------------------------------------------------------- if (state % cfg % do_timing) call timer_start(state % perf % faceloop) face_loop_ausm: do iface = 1, n_local + 1 ! Left-biased stencil of Q for Q_L (same index pattern as F^+ stencil). ! The halo'd ub array is passed with its declared bounds; fill_stencil ! receives field(neq, 1-h:n_local+h), so indices like iface + offset ! that fall into the halo range map to the correct halo cell. call fill_stencil(state % ub, h, & & iface + state % stencil_start_offset, 1, & & q_stencil(:, 1:stencil_width)) ! Hybrid sensor: classify face as smooth or non-smooth. if (state % cfg % use_hybrid_recon) then sensor_val = eval_face_sensor(state, iface, q_stencil(:, 1:stencil_width)) is_smooth = (sensor_val <= state % cfg % hybrid_sensor_threshold) else is_smooth = .false. end if if (is_smooth) then call state % smooth_reconstruct(q_stencil(:, 1:stencil_width), q_face_L) else call state % reconstruct(q_stencil(:, 1:stencil_width), q_face_L) end if if (state % cfg % use_positivity_limiter .and. iface > 1) & call limit_positivity(q_face_L, state % ub(:, iface - 1), state % cfg % gam) ! Right-biased (reversed) stencil of Q for Q_R. call fill_stencil(state % ub, h, & & iface - state % stencil_start_offset - 1, -1, & & q_stencil(:, 1:stencil_width)) if (is_smooth) then call state % smooth_reconstruct(q_stencil(:, 1:stencil_width), q_face_R) else call state % reconstruct(q_stencil(:, 1:stencil_width), q_face_R) end if if (state % cfg % use_positivity_limiter .and. iface <= n_local) & call limit_positivity(q_face_R, state % ub(:, iface), state % cfg % gam) ! Guard: every FDS solver (AUSM+, HLL, HLLC, Roe) computes sqrt(gamma p/rho) ! on the reconstructed faces. A reconstruction overshoot into negative ! density/pressure (e.g. WENO at a strong discontinuity with the ! positivity limiter disabled) would otherwise propagate a silent NaN ! through num_flux. Detect it here and abort cleanly, matching the FVS ! precompute guard below. Placed AFTER limit_positivity so it only fires ! on a still-non-physical face. block real(wp) :: rho_f, p_f rho_f = q_face_L(1) if (rho_f <= 0.0_wp) & call parallel_fatal('compute_resid: non-positive density in FDS reconstructed face (left)') p_f = (q_face_L(3) - 0.5_wp * q_face_L(2)**2 / rho_f) * (state % cfg % gam - 1.0_wp) if (p_f <= 0.0_wp) & call parallel_fatal('compute_resid: non-positive pressure in FDS reconstructed face (left)') rho_f = q_face_R(1) if (rho_f <= 0.0_wp) & call parallel_fatal('compute_resid: non-positive density in FDS reconstructed face (right)') p_f = (q_face_R(3) - 0.5_wp * q_face_R(2)**2 / rho_f) * (state % cfg % gam - 1.0_wp) if (p_f <= 0.0_wp) & call parallel_fatal('compute_resid: non-positive pressure in FDS reconstructed face (right)') end block call state % fds_solver(q_face_L, q_face_R, state % num_flux(:, iface), state % cfg % gam) end do face_loop_ausm if (state % cfg % do_timing) call timer_stop(state % perf % faceloop) else ! ----------------------------------------------------------------------- ! FVS path: split fluxes at cells (interior + halo), reconstruct, sum at ! faces. The split-flux precompute extends to halo cells so that ! fill_stencil reads of fp/fm near the boundary find valid data. ! ----------------------------------------------------------------------- if (state % cfg % do_timing) call timer_start(state % perf % fluxsplit) do i = 1 - h, n_local + h ! Guard: FVS split routines call sqrt(γp/ρ); detect a vacuum or ! non-physical pressure before that to produce a clear error rather ! than propagating NaN through the residual. block real(wp) :: rho_i, p_i rho_i = state % ub(1, i) if (rho_i <= 0.0_wp) & call parallel_fatal('compute_resid: non-positive density in FVS precompute') p_i = (state % ub(3, i) - 0.5_wp * state % ub(2, i)**2 / rho_i) & * (state % cfg % gam - 1.0_wp) if (p_i <= 0.0_wp) & call parallel_fatal('compute_resid: non-positive pressure in FVS precompute') end block call state % split_both(state % ub(:, i), state % fp(:, i), state % fm(:, i), & state % cfg % gam) end do if (state % cfg % do_timing) call timer_stop(state % perf % fluxsplit) ! Loop over all cell faces to compute the numerical flux. if (state % cfg % do_timing) call timer_start(state % perf % faceloop) face_loop: do iface = 1, n_local + 1 ! Compute averaged state at the interface for the eigensystem. ! With halo cells filled, the same Q-average formula works at the ! boundary faces as in the interior — halo cell at index 0 holds the ! left ghost; halo cell at index n_local+1 holds the right ghost. ! See course notes, Sec. "Generalization to Systems", Step 1. q_avg = 0.5_wp * (state % ub(:, iface - 1) + state % ub(:, iface)) ! Step 2: Compute eigenvector matrices at this interface. call compute_eigensystem(q_avg, r_mat, r_inv, state % cfg % gam) ! Hybrid sensor: classify face as smooth or non-smooth using Q (ub), ! not the split-flux stencil, so the same sensor logic applies to ! both the FVS and FDS paths. No q_stencil is available here, so ! weno_beta falls back to density_gradient automatically inside ! eval_face_sensor when the stencil argument has fewer than 5 columns. if (state % cfg % use_hybrid_recon) then sensor_val = eval_face_sensor(state, iface, state % ub(:, 1:0)) is_smooth = (sensor_val <= state % cfg % hybrid_sensor_threshold) else is_smooth = .false. end if ! --- Positive flux F^+: left-biased stencil --- call fill_stencil(state % fp, h, & & iface + state % stencil_start_offset, 1, & & f_stencil(:, 1:stencil_width)) if (is_smooth) then call reconstruct_with_proj(state % smooth_reconstruct, f_stencil(:, 1:stencil_width), & & r_mat, r_inv, state % use_char_proj, flux_pos) else call reconstruct_with_proj(state % reconstruct, f_stencil(:, 1:stencil_width), & & r_mat, r_inv, state % use_char_proj, flux_pos) end if ! --- Negative flux F^-: right-biased (mirror) stencil --- call fill_stencil(state % fm, h, & & iface - state % stencil_start_offset - 1, -1, & & f_stencil(:, 1:stencil_width)) if (is_smooth) then call reconstruct_with_proj(state % smooth_reconstruct, f_stencil(:, 1:stencil_width), & & r_mat, r_inv, state % use_char_proj, flux_neg) else call reconstruct_with_proj(state % reconstruct, f_stencil(:, 1:stencil_width), & & r_mat, r_inv, state % use_char_proj, flux_neg) end if ! Total numerical flux at this face. state % num_flux(:, iface) = flux_pos + flux_neg end do face_loop if (state % cfg % do_timing) call timer_stop(state % perf % faceloop) end if ! Compute the residual: R_i = -1/J_i * (F_{i+1/2} - F_{i-1/2}). ! On uniform grids J_i == dx_uniform exactly, so this is bit-for-bit identical. ! See course notes, Eq. (28). do i = 1, n_local state % resid(:, i) = -(state % num_flux(:, i + 1) - state % num_flux(:, i)) / state % mesh % jac(i) end do if (state % cfg % do_timing) call timer_stop(state % perf % resid) end subroutine compute_resid_fdm ! --------------------------------------------------------------------------- !> Compute the cell-centered finite-volume residual !! R_i = -1/dx_i * (F_{i+1/2} - F_{i-1/2}) on a block of cell averages. !! !! Structurally this is the **FDS path** of `compute_resid_fdm` applied to !! cell averages: reconstruct Q_L and Q_R at each face from the halo'd ub !! array, optionally hybrid-select between the smooth and shock-capturing !! reconstructions, positivity-limit each face, guard against non-physical !! reconstructed states, and call the active FDS / Riemann solver. Two !! deliberate differences from the FDM kernel: !! !! 1. The conservative divide uses the **cell width** `mesh % dx_cell(i)` !! rather than the nodal Jacobian `mesh % jac(i)`. On a uniform grid the !! two coincide, so this reduces to the same operator there. !! 2. FVM is **unconditionally FDS**: a finite-volume block always uses a !! Riemann (FDS) flux, so there is no FVS branch here (config validation !! guarantees an FDS `flux_scheme` for method='fvm'). !! !! The face loop is a faithful copy of the FDM FDS branch (lines in !! `compute_resid_fdm` guarded by `state % use_fds`), reusing the same private !! helpers (`fill_stencil`, `eval_face_sensor`) and the same module imports; !! `compute_resid_fdm` is intentionally left untouched to preserve its !! byte-for-byte FDM behavior, so the loop body is duplicated rather than !! factored into a shared helper. See task-3.1 report for the refactor note. !! !! @param[inout] state Solver instance state. ! --------------------------------------------------------------------------- subroutine compute_resid_fvm_1d(state) use solver_state, only: solver_state_t, neq use reconstruction, only: max_stencil_width use boundary_conditions, only: apply_bcs use halo_exchange, only: exchange_halos use positivity_limiter, only: limit_positivity use mpi_runtime, only: parallel_fatal use timer, only: timer_start, timer_stop type(solver_state_t), intent(inout) :: state integer :: i, iface, h, n_local real(wp) :: q_stencil(neq, max_stencil_width) real(wp) :: q_face_L(neq), q_face_R(neq) integer :: stencil_width real(wp) :: sensor_val logical :: is_smooth stencil_width = state % stencil_width h = state % decomp % halo_width n_local = state % n_pt if (state % cfg % do_timing) call timer_start(state % perf % resid) ! Fill ghost cells: halo_exchange supplies inner-rank halos from neighbour ! interior data, then apply_bcs overwrites the edge-rank physical ghosts ! (dispatching to the FVM cell-centered BC kernel by block method). call exchange_halos(state, state % decomp) call apply_bcs(state) ! ----------------------------------------------------------------------- ! FDS face loop: reconstruct Q_L and Q_R at each face from cell averages. ! ----------------------------------------------------------------------- if (state % cfg % do_timing) call timer_start(state % perf % faceloop) face_loop: do iface = 1, n_local + 1 ! Left-biased stencil of Q for Q_L. call fill_stencil(state % ub, h, & & iface + state % stencil_start_offset, 1, & & q_stencil(:, 1:stencil_width)) ! Hybrid sensor: classify face as smooth or non-smooth. if (state % cfg % use_hybrid_recon) then sensor_val = eval_face_sensor(state, iface, q_stencil(:, 1:stencil_width)) is_smooth = (sensor_val <= state % cfg % hybrid_sensor_threshold) else is_smooth = .false. end if if (is_smooth) then call state % smooth_reconstruct(q_stencil(:, 1:stencil_width), q_face_L) else call state % reconstruct(q_stencil(:, 1:stencil_width), q_face_L) end if if (state % cfg % use_positivity_limiter .and. iface > 1) & call limit_positivity(q_face_L, state % ub(:, iface - 1), state % cfg % gam) ! Right-biased (reversed) stencil of Q for Q_R. call fill_stencil(state % ub, h, & & iface - state % stencil_start_offset - 1, -1, & & q_stencil(:, 1:stencil_width)) if (is_smooth) then call state % smooth_reconstruct(q_stencil(:, 1:stencil_width), q_face_R) else call state % reconstruct(q_stencil(:, 1:stencil_width), q_face_R) end if if (state % cfg % use_positivity_limiter .and. iface <= n_local) & call limit_positivity(q_face_R, state % ub(:, iface), state % cfg % gam) ! Guard: the FDS solver computes sqrt(gamma p/rho) on the reconstructed ! faces; a reconstruction overshoot into negative density/pressure would ! otherwise propagate a silent NaN through num_flux. Placed AFTER ! limit_positivity so it only fires on a still-non-physical face. block real(wp) :: rho_f, p_f rho_f = q_face_L(1) if (rho_f <= 0.0_wp) & call parallel_fatal('compute_resid_fvm_1d: non-positive density in FDS reconstructed face (left)') p_f = (q_face_L(3) - 0.5_wp * q_face_L(2)**2 / rho_f) * (state % cfg % gam - 1.0_wp) if (p_f <= 0.0_wp) & call parallel_fatal('compute_resid_fvm_1d: non-positive pressure in FDS reconstructed face (left)') rho_f = q_face_R(1) if (rho_f <= 0.0_wp) & call parallel_fatal('compute_resid_fvm_1d: non-positive density in FDS reconstructed face (right)') p_f = (q_face_R(3) - 0.5_wp * q_face_R(2)**2 / rho_f) * (state % cfg % gam - 1.0_wp) if (p_f <= 0.0_wp) & call parallel_fatal('compute_resid_fvm_1d: non-positive pressure in FDS reconstructed face (right)') end block call state % fds_solver(q_face_L, q_face_R, state % num_flux(:, iface), state % cfg % gam) end do face_loop if (state % cfg % do_timing) call timer_stop(state % perf % faceloop) ! Conservative finite-volume update: divide by the cell width dx_cell(i) ! (NOT the nodal Jacobian). On a uniform grid dx_cell(i) == dx == jac(i). do i = 1, n_local state % resid(:, i) = -(state % num_flux(:, i + 1) - state % num_flux(:, i)) / state % mesh % dx_cell(i) end do if (state % cfg % do_timing) call timer_stop(state % perf % resid) end subroutine compute_resid_fvm_1d ! --------------------------------------------------------------------------- !> Fill a reconstruction stencil from a halo'd field. !! !! Iterates over stencil positions i = 1 .. size(stencil, 2) and maps each !! to grid index pt_idx = pt_start + (i-1)*pt_step. All indices are !! expected to fall within the halo'd range [1-h, n_local+h]; this is the !! caller's responsibility (guaranteed by the choice of stencil_start_offset !! relative to the scheme's coupling radius). !! !! Use pt_step = +1 for a left-biased stencil (F^+, Q_L reconstruction). !! Use pt_step = -1 for a right-biased (mirror) stencil (F^-, Q_R). !! !! @param field Source halo'd data array, declared with bounds !! (neq, 1-h:n_local+h). !! @param h Halo width (= state%decomp%halo_width). !! @param n_local Interior cell count for this rank (= state%n_pt). !! @param pt_start Grid index for stencil position i = 1. !! @param pt_step Index stride per position: +1 (left-biased) or -1 (right-biased). !! @param stencil Output stencil (neq x stencil_width), overwritten. ! --------------------------------------------------------------------------- pure subroutine fill_stencil(field, h, pt_start, pt_step, stencil) integer, intent(in) :: h !> Halo'd source array. The dummy's column lower bound is explicitly !! set to 1-h so that stencil indices like `iface + stencil_start_offset` !! that fall into the halo range map to the correct cell. real(wp), intent(in) :: field(:, 1 - h:) integer, intent(in) :: pt_start, pt_step real(wp), intent(out) :: stencil(:, :) integer :: i, pt_idx do i = 1, size(stencil, 2) pt_idx = pt_start + (i - 1) * pt_step stencil(:, i) = field(:, pt_idx) end do end subroutine fill_stencil ! --------------------------------------------------------------------------- !> Reconstruct a face flux or state with optional characteristic projection. !! !! If proj is .true.: !! 1. Transform each stencil column to characteristic space: K^{-1} * f. !! 2. Reconstruct the face value in characteristic space. !! 3. Map back to physical space: K * flux_out. !! If proj is .false.: !! Call recon directly on the physical-space stencil. !! !! See course notes, Sec. "Generalization to Systems", Steps 3-4. !! !! @param recon Active reconstruction procedure pointer. !! @param stencil Input flux or state stencil (neq x stencil_width). !! @param r_mat Right eigenvector matrix K (neq x neq). !! @param r_inv Inverse eigenvector matrix K^{-1} (neq x neq). !! @param proj Enable characteristic-space projection. !! @param flux_out Reconstructed face value (neq), overwritten. ! --------------------------------------------------------------------------- subroutine reconstruct_with_proj(recon, stencil, r_mat, r_inv, proj, flux_out) use solver_interfaces, only: reconstructor_iface procedure(reconstructor_iface) :: recon real(wp), intent(in) :: stencil(:, :), r_mat(:, :), r_inv(:, :) logical, intent(in) :: proj real(wp), intent(out) :: flux_out(:) real(wp) :: proj_stencil(size(stencil, 1), size(stencil, 2)) real(wp) :: tmp1, tmp2, tmp3 integer :: k if (proj) then ! Hand-unrolled 3x3 matrix-vector products to avoid runtime library dispatch ! (matmul on assumed-shape arrays generates a _gfortran_matmul_r8 call even at ! -O3 because the rank is not known inside the subroutine; with BLAS linked the ! compiler may additionally dispatch to dgemm, which has non-trivial startup cost ! for a 3x3 kernel). The unrolled form exposes 9 FMAs per column for the ! compiler's vectorizer. do k = 1, size(stencil, 2) proj_stencil(1, k) = r_inv(1, 1) * stencil(1, k) + r_inv(1, 2) * stencil(2, k) & + r_inv(1, 3) * stencil(3, k) proj_stencil(2, k) = r_inv(2, 1) * stencil(1, k) + r_inv(2, 2) * stencil(2, k) & + r_inv(2, 3) * stencil(3, k) proj_stencil(3, k) = r_inv(3, 1) * stencil(1, k) + r_inv(3, 2) * stencil(2, k) & + r_inv(3, 3) * stencil(3, k) end do call recon(proj_stencil, flux_out) ! Back-project: flux_out = r_mat * flux_out (3x3 matmul, hand-unrolled) tmp1 = r_mat(1, 1) * flux_out(1) + r_mat(1, 2) * flux_out(2) + r_mat(1, 3) * flux_out(3) tmp2 = r_mat(2, 1) * flux_out(1) + r_mat(2, 2) * flux_out(2) + r_mat(2, 3) * flux_out(3) tmp3 = r_mat(3, 1) * flux_out(1) + r_mat(3, 2) * flux_out(2) + r_mat(3, 3) * flux_out(3) flux_out(1) = tmp1 flux_out(2) = tmp2 flux_out(3) = tmp3 else call recon(stencil, flux_out) end if end subroutine reconstruct_with_proj ! --------------------------------------------------------------------------- !> Evaluate the hybrid shock sensor for face iface. !! !! Dispatches to the sensor type set in state%cfg%hybrid_sensor: !! 'jameson' — Jameson-Schmidt-Turkel pressure second-derivative !! sensor (JST 1981, AIAA Paper 81-1259). !! 'density_gradient' — Normalised density jump across the face. !! 'weno_beta' — WENO5 smoothness-indicator ratio on density !! (Jiang & Shu 1996, Eq. 24-26); falls back to !! density_gradient for non-5-point schemes or when !! q_stencil has < 5 columns. !! !! @param state Solver state (solution arrays and config). !! @param iface Face index in [1, n_pt+1]. !! @param q_stencil Assembled Q stencil (neq x ≥5) for weno_beta sensor; !! pass a zero-column slice for the FVS path. !! @return Non-negative sensor value; larger = more non-smooth. ! --------------------------------------------------------------------------- real(wp) function eval_face_sensor(state, iface, q_stencil) result(s) use solver_state, only: solver_state_t, neq use weno_family, only: weno5_smoothness, eps_weno type(solver_state_t), intent(in) :: state integer, intent(in) :: iface real(wp), intent(in) :: q_stencil(:, :) real(wp) :: q_m2(neq), q_m1(neq), q_p1(neq), q_p2(neq) real(wp) :: p_m2, p_m1, p_p1, p_p2, gm1 real(wp) :: rho_stencil(1, 5), beta0(1), beta1(1), beta2(1) q_m1 = get_q_at_face(state, iface - 1) q_p1 = get_q_at_face(state, iface) select case (trim(state % cfg % hybrid_sensor)) case (hybrid_sensor_jameson) q_m2 = get_q_at_face(state, iface - 2) q_p2 = get_q_at_face(state, iface + 1) gm1 = state % cfg % gam - 1.0_wp p_m2 = gm1 * (q_m2(3) - 0.5_wp * q_m2(2)**2 / q_m2(1)) p_m1 = gm1 * (q_m1(3) - 0.5_wp * q_m1(2)**2 / q_m1(1)) p_p1 = gm1 * (q_p1(3) - 0.5_wp * q_p1(2)**2 / q_p1(1)) p_p2 = gm1 * (q_p2(3) - 0.5_wp * q_p2(2)**2 / q_p2(1)) s = (abs(p_m2 - 2.0_wp * p_m1 + p_p1) + abs(p_m1 - 2.0_wp * p_p1 + p_p2)) & / (abs(p_m2) + 2.0_wp * abs(p_m1) + 2.0_wp * abs(p_p1) + abs(p_p2) & + tiny(1.0_wp)) case (hybrid_sensor_density_gradient) s = abs(q_p1(1) - q_m1(1)) & / (0.5_wp * (abs(q_p1(1)) + abs(q_m1(1))) + tiny(1.0_wp)) case (hybrid_sensor_weno_beta) if (state % stencil_width == 5 .and. size(q_stencil, 2) >= 5) then ! WENO5 smoothness-indicator ratio on density (equation 1). ! Large value indicates a discontinuity in the reconstruction stencil. rho_stencil(1, :) = q_stencil(1, 1:5) call weno5_smoothness(rho_stencil, beta0, beta1, beta2) s = max(beta0(1), beta1(1), beta2(1)) & / (min(beta0(1), beta1(1), beta2(1)) + eps_weno) else ! No Q stencil (FVS path) — fall back to density_gradient. s = abs(q_p1(1) - q_m1(1)) & / (0.5_wp * (abs(q_p1(1)) + abs(q_m1(1))) + tiny(1.0_wp)) end if case default s = 0.0_wp ! unreachable: config validation guards this path end select end function eval_face_sensor ! --------------------------------------------------------------------------- !> Return Q at grid cell idx, reading directly from the halo'd state%ub. !! !! Indices in [1-h, n_local+h] read the appropriate halo or interior cell. !! For periodic single-rank runs apply_bcs has already written the wrap into !! the halo cells; for multi-rank runs halo_exchange handles it. Indices !! that fall outside that range (e.g. iface + 1 == n_local + 2 when !! h == 1) are clamped to the nearest halo cell as a defensive fallback. !! !! @param state Solver state. !! @param idx Cell index (may lie outside [1, n_pt]). !! @return Q(neq) for that cell, or the appropriate ghost halo cell. ! --------------------------------------------------------------------------- function get_q_at_face(state, idx) result(q) use solver_state, only: solver_state_t, neq type(solver_state_t), intent(in) :: state integer, intent(in) :: idx real(wp) :: q(neq) integer :: lo, hi, idx_use lo = lbound(state % ub, 2) hi = ubound(state % ub, 2) idx_use = max(lo, min(hi, idx)) q = state % ub(:, idx_use) end function get_q_at_face end module spatial_discretization