!> @file solver_state.f90 !> @brief Solver instance state for the 1D Euler solver. !! !! Defines the @p solver_state_t derived type that bundles all per-simulation !! data: grid parameters, physical constants, time-stepping scalars, boundary !! conditions, ghost-state vectors, scheme dispatch pointers, and allocatable !! solution arrays. !! !! Replacing the former module-level global variables in @p euler_constants with !! an explicit derived type enables: !! - Multiple concurrent solver instances (e.g. Richardson extrapolation). !! - Self-contained unit tests with no shared global state. !! - Clear data ownership and lifetime through Fortran's allocatable semantics. !! !! The scheme dispatch pointers (@p reconstruct, @p smooth_reconstruct, !! @p flux_split, @p fds_solver) and their associated metadata !! (@p stencil_width, @p stencil_start_offset, @p coupling_radius, !! @p use_char_proj, @p use_fds) are owned per-instance so that two solver !! instances can run different schemes concurrently without interfering with !! each other. !! !! The integer parameter @p neq = 3 is exported from this module so that !! modules needing only the equation count can import it without dragging in !! the full derived type. !! !! Usage (driver): !! @code !! type(config_t) :: cfg !! type(solver_state_t) :: state !! call read_config('input.nml', cfg) !! call init_from_config(state, cfg) !! call allocate_work_arrays(state) !! call init_flux_scheme(state, 'lax_friedrichs') !! call init_recon_scheme(state, 'weno5') !! @endcode module solver_state use precision, only: wp use config, only: config_t use solver_interfaces, only: reconstructor_iface, flux_splitter_iface, & flux_splitter_both_iface, fds_iface use timer, only: timer_t implicit none private !> Number of conserved equations for 1D Euler: (rho, rho*u, E). integer, parameter, public :: neq = 3 public :: init_from_config, allocate_work_arrays, release_work_arrays !> Named wall-clock timers for the three hot regions of compute_resid(). !! !! Populated only when solver_state_t%do_timing is .true.; all timers are !! zero-initialised by default so they are safe to read even if unused. type, public :: perf_counters_t !> Time in compute_resid() (both FVS and FDS paths). type(timer_t) :: resid !> Time in the flux-split precompute loop (FVS path only). type(timer_t) :: fluxsplit !> Time in the face reconstruction loop (both paths). type(timer_t) :: faceloop end type perf_counters_t !> All per-simulation state for the 1D Euler solver. !! !! Immutable simulation parameters are stored in the embedded `cfg` field !! (a copy of the `config_t` struct); call `init_from_config` to populate it !! and the derived grid scalars `n_pt`, `dx`. Allocatable arrays must be !! explicitly allocated before use (see @p euler_1d and !! @p test_helpers). Procedure pointer components are initialised to null !! and must be bound via init_flux_scheme() and init_recon_scheme() before !! any residual call. type, public :: solver_state_t ! ------------------------------------------------------------------------- ! Immutable simulation parameters (copied from config_t at init time) ! ------------------------------------------------------------------------- !> Configuration object (copy of the parsed input). Access physics, !! grid bounds, BC types, scheme names, etc. via state%cfg%<field>. type(config_t) :: cfg ! ------------------------------------------------------------------------- ! Derived grid scalars (computed by init_from_config) ! ------------------------------------------------------------------------- !> Number of grid points = cfg%n_cell + 1. integer :: n_pt = 0 !> Uniform grid spacing Δx = (cfg%x_right - cfg%x_left) / cfg%n_cell [m]. real(wp) :: dx = 0.0_wp ! ------------------------------------------------------------------------- ! Mutable time-stepping scalars ! ------------------------------------------------------------------------- !> Current time step size Δt [s]. Initialised from cfg%dt; updated each !! step by CFL control or last-step clipping. real(wp) :: dt = 0.0_wp !> Global L2 residual norm, updated at the end of every time step. real(wp) :: resid_glob = 0.0_wp ! ------------------------------------------------------------------------- ! Derived BC flag ! ------------------------------------------------------------------------- !> True when periodic BCs are active (set by apply_bcs()). logical :: is_periodic = .false. ! ------------------------------------------------------------------------- ! Ghost conserved states (fixed size = neq) ! ------------------------------------------------------------------------- !> Left ghost conserved state Q_L = (rho, rho*u, E). real(wp) :: q_left(neq) = 0.0_wp !> Right ghost conserved state Q_R. real(wp) :: q_right(neq) = 0.0_wp ! ------------------------------------------------------------------------- ! Ghost split fluxes at boundaries (FVS path only) ! ------------------------------------------------------------------------- !> F^+ at the left boundary ghost cell. real(wp) :: fl_pos(neq) = 0.0_wp !> F^- at the left boundary ghost cell. real(wp) :: fl_neg(neq) = 0.0_wp !> F^+ at the right boundary ghost cell. real(wp) :: fr_pos(neq) = 0.0_wp !> F^- at the right boundary ghost cell. real(wp) :: fr_neg(neq) = 0.0_wp ! ------------------------------------------------------------------------- ! Scheme dispatch (per-instance; set by init_flux_scheme / init_recon_scheme) ! ------------------------------------------------------------------------- !> Procedure pointer to the active nonlinear reconstruction scheme. !! Set by init_recon_scheme(); null until then. procedure(reconstructor_iface), pointer, nopass :: reconstruct => null() !> Procedure pointer to the smooth-region reconstruction used by hybrid mode. !! For schemes without a dedicated smooth-region companion, this is bound to !! the same procedure as @p reconstruct. procedure(reconstructor_iface), pointer, nopass :: smooth_reconstruct => null() !> Procedure pointer to the active flux-splitting scheme (FVS path). !! Set by init_flux_scheme(); null when a FDS scheme is active. procedure(flux_splitter_iface), pointer, nopass :: flux_split => null() !> Fused flux-splitting procedure that computes F^+ and F^- in one call. !! Avoids extracting primitives twice per cell in the FVS precompute loop. !! Set by init_flux_scheme() alongside flux_split; null for FDS schemes. procedure(flux_splitter_both_iface), pointer, nopass :: split_both => null() !> Procedure pointer to the active FDS approximate Riemann solver. !! Set by init_flux_scheme(); null when a FVS scheme is active. procedure(fds_iface), pointer, nopass :: fds_solver => null() !> Active reconstruction stencil width. !! Examples: WENO5 = 5, WENO-CU6 = 6, WENO11 = 11. integer :: stencil_width = 5 !> Left-biased stencil start offset relative to the face index. !! For example, WENO5 uses -3 so the left-biased stencil starts at iface-3 !! and covers columns iface-3 : iface+1. integer :: stencil_start_offset = -3 !> Cell-to-cell coupling radius induced by the active reconstruction. !! Used by implicit steppers to size the banded Jacobian. integer :: coupling_radius = 3 !> When .true., spatial_discretization applies K^{-1}/K projection around !! the reconstruct() call. Set by init_recon_scheme() based on char_proj. logical :: use_char_proj = .false. !> True when a FDS scheme (AUSM+, HLL, HLLC, Roe) is active. !! When false the FVS path via flux_split is used. logical :: use_fds = .false. ! ------------------------------------------------------------------------- ! Allocatable solution arrays ! ------------------------------------------------------------------------- !> Conserved variables Q_i (neq × n_pt). real(wp), allocatable :: ub(:, :) !> Numerical flux at cell faces (neq × (n_pt+1)). real(wp), allocatable :: num_flux(:, :) !> Spatial residual R(Q) (neq × n_pt). real(wp), allocatable :: resid(:, :) !> Positive split flux F^+ (neq × n_pt). real(wp), allocatable :: fp(:, :) !> Negative split flux F^- (neq × n_pt). real(wp), allocatable :: fm(:, :) ! ------------------------------------------------------------------------- ! Profiling ! ------------------------------------------------------------------------- !> Accumulated wall-clock timers; populated only when cfg%do_timing is .true.. type(perf_counters_t) :: perf ! ------------------------------------------------------------------------- ! Scratch arrays for explicit Runge-Kutta steppers ! ------------------------------------------------------------------------- !> Scratch array 1 (neq × n_pt): stage save for SSPRK22, TVD-RK3, RK4, SSPRK54. !! Pre-allocated once to avoid per-step heap allocations in the time loop. real(wp), allocatable :: scratch1(:, :) !> Scratch array 2 (neq × n_pt): second stage save for RK4 (k_sum) and SSPRK54 (ub_stage2). real(wp), allocatable :: scratch2(:, :) ! ------------------------------------------------------------------------- ! BDF2 bootstrap state (owned per-instance to allow concurrent solvers) ! ------------------------------------------------------------------------- !> True after the first bdf2_step() call (bootstrap complete). logical :: bdf2_initialized = .false. !> Q^{n-1} storage for BDF2. Allocated lazily on first bdf2_step() call. real(wp), allocatable :: bdf2_ub_prev(:, :) contains !> Nullify procedure pointer components when the instance is finalised. !! Allocatable array components are automatically deallocated by Fortran. final :: destroy_solver_state end type solver_state_t contains ! --------------------------------------------------------------------------- !> Populate a solver_state_t from a config_t. !! !! Embeds a copy of @p cfg into state%cfg and computes the two derived grid !! scalars state%n_pt and state%dx. Also initialises state%dt from cfg%dt, !! applies problem-type-specific BC promotions (smooth_wave → periodic, !! woodward_colella → reflecting), and sets state%is_periodic. !! !! Call this once after read_config() and before init_flux_scheme() / !! init_recon_scheme() / apply_initial_condition(). !! !! @param[inout] state Solver state to populate. !! @param[in] cfg Parsed configuration. ! --------------------------------------------------------------------------- subroutine init_from_config(state, cfg) type(solver_state_t), intent(inout) :: state type(config_t), intent(in) :: cfg state % cfg = cfg state % n_pt = cfg % n_cell + 1 state % dx = (cfg % x_right - cfg % x_left) / real(cfg % n_cell, wp) state % dt = cfg % dt ! Problem-type-driven BC promotions (override user default 'dirichlet'). if (trim(cfg % problem_type) == 'smooth_wave' .or. & trim(cfg % problem_type) == 'linear_advection') then if (trim(state % cfg % bc_left) == 'dirichlet') state % cfg % bc_left = 'periodic' if (trim(state % cfg % bc_right) == 'dirichlet') state % cfg % bc_right = 'periodic' end if if (trim(cfg % problem_type) == 'woodward_colella') then if (trim(state % cfg % bc_left) == 'dirichlet') state % cfg % bc_left = 'reflecting' if (trim(state % cfg % bc_right) == 'dirichlet') state % cfg % bc_right = 'reflecting' end if state % is_periodic = (trim(state % cfg % bc_left) == 'periodic') end subroutine init_from_config ! --------------------------------------------------------------------------- !> Allocate all solver work arrays owned by a solver instance. !! !! This centralises array ownership so drivers and tests do not duplicate the !! same allocation block. Scratch arrays are preallocated here as part of the !! standard solver lifecycle. ! --------------------------------------------------------------------------- subroutine allocate_work_arrays(state) type(solver_state_t), intent(inout) :: state integer :: info if (allocated(state % ub)) & error stop 'solver_state: allocate_work_arrays called on already-allocated state' allocate (state % ub(neq, state % n_pt), stat=info) if (info /= 0) error stop 'solver_state: allocation of ub failed' allocate (state % num_flux(neq, state % n_pt + 1), stat=info) if (info /= 0) error stop 'solver_state: allocation of num_flux failed' allocate (state % fp(neq, state % n_pt), state % fm(neq, state % n_pt), & state % resid(neq, state % n_pt), stat=info) if (info /= 0) error stop 'solver_state: allocation of fp/fm/resid failed' allocate (state % scratch1(neq, state % n_pt), & state % scratch2(neq, state % n_pt), stat=info) if (info /= 0) error stop 'solver_state: allocation of scratch arrays failed' end subroutine allocate_work_arrays ! --------------------------------------------------------------------------- !> Release all allocatable arrays owned by a solver instance. ! --------------------------------------------------------------------------- subroutine release_work_arrays(state) type(solver_state_t), intent(inout) :: state integer :: info if (allocated(state % ub)) then deallocate (state % ub, stat=info) if (info /= 0) error stop 'solver_state: deallocation of ub failed' end if if (allocated(state % num_flux)) then deallocate (state % num_flux, stat=info) if (info /= 0) error stop 'solver_state: deallocation of num_flux failed' end if if (allocated(state % fp)) then deallocate (state % fp, stat=info) if (info /= 0) error stop 'solver_state: deallocation of fp failed' end if if (allocated(state % fm)) then deallocate (state % fm, stat=info) if (info /= 0) error stop 'solver_state: deallocation of fm failed' end if if (allocated(state % resid)) then deallocate (state % resid, stat=info) if (info /= 0) error stop 'solver_state: deallocation of resid failed' end if if (allocated(state % scratch1)) then deallocate (state % scratch1, stat=info) if (info /= 0) error stop 'solver_state: deallocation of scratch1 failed' end if if (allocated(state % scratch2)) then deallocate (state % scratch2, stat=info) if (info /= 0) error stop 'solver_state: deallocation of scratch2 failed' end if if (allocated(state % bdf2_ub_prev)) then deallocate (state % bdf2_ub_prev, stat=info) if (info /= 0) error stop 'solver_state: deallocation of bdf2_ub_prev failed' end if end subroutine release_work_arrays ! --------------------------------------------------------------------------- !> Nullify procedure pointer components of a solver_state_t instance. !! !! Called automatically by the compiler when a solver_state_t goes out of !! scope or is deallocated. Allocatable array components (ub, num_flux, !! resid, fp, fm, bdf2_ub_prev) are deallocated automatically by Fortran !! without any action here. !! !! @param[inout] self The solver state being finalised. ! --------------------------------------------------------------------------- subroutine destroy_solver_state(self) type(solver_state_t), intent(inout) :: self nullify (self % reconstruct) nullify (self % flux_split) nullify (self % split_both) nullify (self % fds_solver) end subroutine destroy_solver_state end module solver_state