!> @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 halo_width, !! @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 init_flux_scheme(state, 'lax_friedrichs') !! call init_recon_scheme(state, 'weno5') ! sets state%halo_width !! state%decomp = decompose(my_rank(), n_ranks(), & !! state%n_pt, state%halo_width, state%is_periodic) !! call allocate_work_arrays(state) ! requires state%decomp populated !! @endcode module solver_state use precision, only: wp use config, only: config_t use domain_decomposition, only: decomp_t use solver_interfaces, only: reconstructor_iface, flux_splitter_iface, & flux_splitter_both_iface, fds_iface use timer, only: timer_t use mpi_runtime, only: my_rank use mesh_1d, only: mesh_1d_t, build_mesh_global, build_mesh_local, & build_mesh_cellcentered_global, build_mesh_cellcentered_local use block_mod, only: block_t use option_registry, only: method_fdm, method_fvm 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) ! ------------------------------------------------------------------------- !> Per-rank interior count (= decomp%n_local at runtime). For FDM this is a !! node count (= cfg%n_cell+1 at -np 1); for FVM a cell count (= cfg%n_cell). integer :: n_pt = 0 !> Total global interior count across all ranks (= decomp%n_global). FDM: !! cfg%n_cell+1 (nodes); FVM: cfg%n_cell (cells). Used to normalise the !! global L2 residual norm in compute_resid_glob. integer :: n_pt_global = 0 !> Grid metric layer (node coordinates + dx/dξ Jacobian). Built by !! init_from_config (global) and allocate_work_arrays (local slice). type(mesh_1d_t) :: mesh !> Per-rank decomposition descriptor. Populated by solver_runtime at !! session init (Task 6 of the Phase B plan); read by halo_exchange, !! boundary_conditions, and any other rank-aware code path. At -np 1 !! this defaults to a single-rank, self-contained decomposition. type(decomp_t) :: decomp ! ------------------------------------------------------------------------- ! Block discretization descriptors ! ------------------------------------------------------------------------- !> Per-block discretization metadata (method + cell count). A single !! block spans the whole 1D domain today; `blocks(1) % method` selects the !! residual path via internal `select case` in spatial_discretization. !! Populated by init_from_config once cfg is available. type(block_t) :: blocks(1) ! ------------------------------------------------------------------------- ! 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. ! ------------------------------------------------------------------------- ! Boundary ghost state vectors and split fluxes were removed in Phase B ! Task 8. Boundary conditions now write directly into the halo cells of ! state%ub (and FVS split fluxes are computed over the halo range each ! residual evaluation), so no per-state scalar slots are needed. ! ------------------------------------------------------------------------- ! 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 !> Halo width: number of additional cells on each side of a rank's !! interior region that must be exchanged from neighbours to drive the !! reconstruction stencil. Sized by the chosen scheme's stencil metadata !! (determined by init_recon_scheme). integer :: halo_width = 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(:, :) !> Third RK stage-save scratch (ssprk54: holds u3; the L(u3) residual !! stays in resid, with its final-stage piece accumulated early). real(wp), allocatable :: scratch3(:, :) ! ------------------------------------------------------------------------- ! 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(:, :) ! ------------------------------------------------------------------------- ! Gather buffer for output writers (Phase D) ! ------------------------------------------------------------------------- !> Rank-0-owned gather buffer for write_solution_file / write_snapshot. !! Shape (neq, n_pt_global) on rank 0; (neq, 0) on every other rank. !! Allocated by allocate_work_arrays once decomp is known and reused on !! every writer call to avoid per-call alloc/dealloc churn. real(wp), allocatable :: gather_buf(:, :) 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 block logical :: mesh_ok character(len=256) :: mesh_msg call build_mesh_global(state % mesh, state % cfg % grid_type, state % cfg % grid_file, & state % cfg % n_cell, state % cfg % x_left, state % cfg % x_right, & mesh_ok, mesh_msg) if (.not. mesh_ok) error stop trim(mesh_msg) ! FVM on a non-uniform (file) grid: build the replicated GLOBAL ! cell-centered mesh from the cell faces that build_mesh_global has just ! loaded into mesh % x_global. allocate_work_arrays then slices this per ! rank. The uniform FVM path skips this (its local builder recomputes cell ! centers analytically, bit-for-bit with the original uniform formula), and ! the FDM path never touches the cell-centered arrays at all. if (trim(cfg % method) == method_fvm .and. .not. state % mesh % uniform) then call build_mesh_cellcentered_global(state % mesh, state % cfg % n_cell) end if end block state % n_pt = state % cfg % n_cell + 1 ! per-rank count; overwritten by allocate_work_arrays state % n_pt_global = state % cfg % n_cell + 1 ! global count; overwritten by allocate_work_arrays state % dt = cfg % dt ! Tag the (single) block with its discretization method and cell count. ! cfg has been validated upstream by read_config; the residual dispatcher ! reads blocks(1) % method to select the FDM or FVM path. ! block_t % method is sized to the canonical option_registry tokens ! (len=8); cfg % method is the wider raw config string (len=64), so copy ! the leading 8 chars explicitly — the validated value is 'fdm'/'fvm', well ! within 8 chars, and the slice keeps the assignment truncation-warning-free. ! Guard against a wider validated token silently losing characters in the ! len=8 slice: the assignment below would truncate, and the residual/mesh ! dispatch downstream would then mis-route on a corrupted method string. if (len_trim(cfg % method) > 8) & error stop 'init_from_config: method token exceeds block_t field width (8)' state % blocks(1) % method = cfg % method(1:8) state % blocks(1) % n_cell = cfg % n_cell ! 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' ! solver_runtime must instantiate state%decomp (via decompose()) before ! calling allocate_work_arrays. Enforce that invariant here. if (state % decomp % n_local == -1) then error stop 'solver_state: allocate_work_arrays called before solver_runtime instantiated state%decomp' end if ! Update per-rank and global point counts from the decomposition. state % n_pt = state % decomp % n_local state % n_pt_global = state % decomp % n_global ! Build this rank's local mesh in the layout the block's method needs: the ! nodal mesh for FDM, the cell-centered mesh for FVM. The work-array shapes ! below are halo-padded and shared by BOTH methods (n_local is the per-rank ! node count for FDM, cell count for FVM), so only the mesh build differs. select case (trim(state % blocks(1) % method)) case (method_fdm) call build_mesh_local(state % mesh, state % decomp) case (method_fvm) call build_mesh_cellcentered_local(state % mesh, state % decomp, & state % cfg % x_left, state % cfg % x_right) case default error stop 'solver_state: allocate_work_arrays unknown block method "'// & trim(state % blocks(1) % method)//'"' end select associate (n_local => state % decomp % n_local, h => state % decomp % halo_width) allocate (state % ub(neq, 1 - h:n_local + h), stat=info) if (info /= 0) error stop 'solver_state: allocate ub failed' state % ub = 0.0_wp allocate (state % num_flux(neq, 0:n_local + 1), stat=info) if (info /= 0) error stop 'solver_state: allocate num_flux failed' state % num_flux = 0.0_wp allocate (state % fp(neq, 1 - h:n_local + h), & state % fm(neq, 1 - h:n_local + h), & state % resid(neq, 1 - h:n_local + h), stat=info) if (info /= 0) error stop 'solver_state: allocate fp/fm/resid failed' state % fp = 0.0_wp state % fm = 0.0_wp state % resid = 0.0_wp allocate (state % scratch1(neq, 1 - h:n_local + h), & state % scratch2(neq, 1 - h:n_local + h), & state % scratch3(neq, 1 - h:n_local + h), stat=info) if (info /= 0) error stop 'solver_state: allocate scratch arrays failed' state % scratch1 = 0.0_wp state % scratch2 = 0.0_wp state % scratch3 = 0.0_wp end associate if (my_rank() == 0) then allocate (state % gather_buf(neq, state % decomp % n_global), stat=info) else allocate (state % gather_buf(neq, 0), stat=info) end if if (info /= 0) error stop 'solver_state: allocate gather_buf failed' state % gather_buf = 0.0_wp 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 % mesh % x_node)) then deallocate (state % mesh % x_node, stat=info) if (info /= 0) error stop 'solver_state: deallocation of mesh x_node failed' end if if (allocated(state % mesh % jac)) then deallocate (state % mesh % jac, stat=info) if (info /= 0) error stop 'solver_state: deallocation of mesh jac failed' end if if (allocated(state % mesh % x_cell)) then deallocate (state % mesh % x_cell, stat=info) if (info /= 0) error stop 'solver_state: deallocation of mesh x_cell failed' end if if (allocated(state % mesh % dx_cell)) then deallocate (state % mesh % dx_cell, stat=info) if (info /= 0) error stop 'solver_state: deallocation of mesh dx_cell failed' end if 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 % scratch3)) then deallocate (state % scratch3, stat=info) if (info /= 0) error stop 'solver_state: deallocation of scratch3 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 if (allocated(state % gather_buf)) then deallocate (state % gather_buf, stat=info) if (info /= 0) error stop 'solver_state: deallocation of gather_buf 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, scratch1, scratch2, scratch3, bdf2_ub_prev, gather_buf) 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