!> Solver execution lifecycle and I/O orchestration. !! !! This module centralises the solver lifecycle behind the thin !! `app/euler_1d.f90` program: CLI parsing, config loading, solver-state !! initialisation, adaptive-dt time loop, checkpoint/restart, live snapshots, !! iteration logging, and performance reporting. !! !! Typical call sequence: !! ```fortran !! call resolve_cli_namelist(ctx % nml_file) !! call read_config(ctx % nml_file, cfg) !! call initialize_runtime(ctx, cfg, ctx % nml_file) !! call run_solver(ctx) !! call finalize_runtime(ctx) !! call teardown_runtime(ctx) !! ``` module solver_runtime use run_step_engine, only: run_loop_ctx_t, run_bounded_steps, progress_callback_i, tick_seconds use precision, only: wp, safe_vel use config, only: config_t use option_registry, only: problem_sod, problem_lax, problem_woodward_colella, & problem_shu_osher, problem_acoustic_pulse, & problem_from_file, problem_udf, method_fdm, method_fvm use solver_state, only: solver_state_t, allocate_work_arrays, release_work_arrays, & init_from_config, neq use euler_physics, only: init_flux_scheme, compute_max_wave_speed use parallel_reductions, only: par_max_real, par_lor use solution_gather, only: gather_solution_to_root use reconstruction, only: init_recon_scheme use mpi_runtime, only: my_rank, n_ranks, parallel_fatal use domain_decomposition, only: decompose, validate_decomp use time_integration, only: stepper_iface, resolve_time_scheme use initial_conditions, only: apply_initial_condition use timer, only: timer_t, timer_start, timer_stop, timer_elapsed_s, & timer_elapsed_running_s, timer_reset, timer_report use checkpoint, only: write_checkpoint, read_checkpoint use logger, only: log_init, log_finalize, log_info, log_warn use run_banner, only: log_run_banner use tecplot_writer, only: write_tecplot_ascii, write_tecplot_binary use output_format_list, only: parse_format_list, derived_filename, max_formats implicit none private public :: solver_run_context_t, resolve_cli_namelist, resolve_cli, consume_cli_arg public :: initialize_runtime, run_solver, run_solver_steps public :: finalize_runtime, report_performance_summary, teardown_runtime, copy_current_solution public :: write_solution_file ! Re-exported from run_step_engine (Q2c) so existing consumers of the ! historical solver_runtime names are unaffected. public :: progress_callback_i public :: tick_seconds !> Aggregate context bundling all state needed for a single solver run. !! !! One instance is owned by the driver and passed by reference through the !! entire lifecycle. Fields are populated in order by `resolve_cli_namelist`, !! `read_config`, and `initialize_runtime`. !! The authoritative runtime configuration is `state % cfg`; all config reads !! go through that field after `initialize_runtime` returns. !! `t`, `t_comp`, `iter`, and `run_complete` are inherited from !! `run_loop_ctx_t` (Q2c) under their historical names. type, extends(run_loop_ctx_t), public :: solver_run_context_t type(solver_state_t) :: state !< Solution arrays and scheme handles. type(timer_t) :: t_total !< Wall-clock timer for the entire run. type(timer_t) :: t_io !< Wall-clock timer for output I/O. type(timer_t) :: t_iter !< Wall-clock timer for a single step. integer :: print_freq = 50 !< Iteration interval between log lines. character(len=512) :: nml_file = 'input.nml' !< Path to the namelist file. procedure(stepper_iface), pointer, nopass :: stepper => null() !< Bound time integrator. logical :: timers_started = .false. !< True after starting the total timer. logical :: performance_reported = .false. !< True after printing timings. contains procedure :: raw_dt => ctx1d_raw_dt procedure :: set_dt => ctx1d_set_dt procedure :: step => ctx1d_step procedure :: write_checkpoint_now => ctx1d_checkpoint procedure :: log_iteration_line => ctx1d_log_line procedure :: begin_work => ctx1d_begin_work procedure :: write_snapshot_now => ctx1d_snapshot procedure :: finish_work => ctx1d_finish_work procedure :: progress_residual => ctx1d_progress_residual end type solver_run_context_t contains !> Read the namelist filename from the first command-line argument. !! Falls back to `'input.nml'` when no argument is supplied. !! @param nml_file On return: path to the namelist file. subroutine resolve_cli_namelist(nml_file) character(len=*), intent(out) :: nml_file integer :: arg_len, arg_stat call get_command_argument(1, nml_file, arg_len, arg_stat) if (arg_stat /= 0 .or. arg_len == 0) nml_file = 'input.nml' end subroutine resolve_cli_namelist !> Parse the worker's CLI: --cortex-socket=URI plus a positional !! namelist path. Either may be absent. !! @param nml_file On return: namelist path (defaults to 'input.nml'). !! @param cortex_uri On return: URI string if --cortex-socket was given, !! or empty. !! @param attached On return: .true. iff --cortex-socket was given. subroutine resolve_cli(nml_file, cortex_uri, attached) character(len=*), intent(out) :: nml_file, cortex_uri logical, intent(out) :: attached integer :: i, arg_len, arg_stat character(len=512) :: arg nml_file = '' cortex_uri = '' attached = .false. do i = 1, command_argument_count() call get_command_argument(i, arg, arg_len, arg_stat) if (arg_stat /= 0) cycle call consume_cli_arg(arg, arg_len, nml_file, cortex_uri, attached) end do if (.not. attached .and. len_trim(nml_file) == 0) then nml_file = 'input.nml' end if end subroutine resolve_cli !> Dispatch one parsed CLI token onto the worker-CLI state. !! !! `get_command_argument` returns `arg_len` as the argument's TRUE length, !! which can exceed the fixed `arg` buffer for an over-long argument; the !! length is therefore clamped to `len(arg)` before any substring is taken !! so `arg(17:arg_len)` / `arg(1:arg_len)` can never index past the buffer !! (an out-of-bounds read otherwise). An over-long value is truncated to the !! buffer rather than overflowing it. Factored out of `resolve_cli` so the !! clamp is exercised directly by the runtime-check (bounds) gate. !! !! @param arg The argument buffer as filled by `get_command_argument`. !! @param arg_len The argument's reported TRUE length (may exceed !! `len(arg)`); the routine clamps it internally to !! `len(arg)` before slicing — callers pass it as-is. !! @param nml_file Updated in place with a positional namelist path. !! @param cortex_uri Updated in place with a `--cortex-socket=` value. !! @param attached Set `.true.` iff a `--cortex-socket=` token was seen. pure subroutine consume_cli_arg(arg, arg_len, nml_file, cortex_uri, attached) character(len=*), intent(in) :: arg integer, intent(in) :: arg_len character(len=*), intent(inout) :: nml_file, cortex_uri logical, intent(inout) :: attached integer :: eff_len eff_len = min(arg_len, len(arg)) if (eff_len < 1) return if (eff_len >= 16) then if (arg(1:16) == "--cortex-socket=") then cortex_uri = arg(17:eff_len) attached = .true. return end if end if if (arg(1:1) /= "-") nml_file = arg(1:eff_len) end subroutine consume_cli_arg !> Populate `ctx` and bring the solver to a ready-to-run state. !! !! Steps performed in order: copy config into context; initialise logger; !! initialise solver state from config; bind flux and recon scheme pointers !! (sets state%halo_width); instantiate the domain decomposition (decompose + !! validate_decomp, uses halo_width); allocate work arrays (uses decomp); !! bind the time-integration scheme pointer; apply initial conditions; !! optionally load a restart checkpoint; log the effective configuration. !! !! @param ctx Context to initialise. !! @param cfg Configuration loaded via `read_config`. !! @param nml_file Optional namelist path stored in `ctx % nml_file`. subroutine initialize_runtime(ctx, cfg, nml_file, is_ok, message) type(solver_run_context_t), intent(inout) :: ctx type(config_t), intent(in) :: cfg character(len=*), intent(in), optional :: nml_file logical, intent(out), optional :: is_ok character(len=*), intent(out), optional :: message logical :: ok character(len=256) :: err integer :: n_global_decomp if (present(is_ok)) is_ok = .true. if (present(message)) message = '' if (present(nml_file)) ctx % nml_file = nml_file call log_init(cfg % verbosity, cfg % log_file) call log_run_banner() call log_info('config: loaded "'//trim(ctx % nml_file)//'"') call init_from_config(ctx % state, cfg) ! init_recon_scheme must run before the decompose() call below so that ! state%halo_width is correctly resolved from the chosen scheme before ! decompose() uses it. call init_flux_scheme(ctx % state, ctx % state % cfg % flux_scheme) call init_recon_scheme(ctx % state, ctx % state % cfg % recon_scheme, & ctx % state % cfg % char_proj, ctx % state % cfg % limiter) ! Instantiate the per-rank decomposition from the real MPI identity ! (Phase B Task 6). Must run after init_recon_scheme (sets halo_width) ! and before allocate_work_arrays (consumes decomp%n_local). ! ! The decomposition's global count is the one array-related quantity that ! differs by method: FDM is node-based (n_cell + 1 nodes), FVM is cell-based ! (n_cell cells). Everything downstream (n_local, the halo-padded work-array ! shapes, n_pt/n_pt_global) flows from this count. select case (trim(ctx % state % blocks(1) % method)) case (method_fdm) n_global_decomp = ctx % state % cfg % n_cell + 1 case (method_fvm) n_global_decomp = ctx % state % cfg % n_cell case default call parallel_fatal('initialize_runtime: unknown block method "'// & trim(ctx % state % blocks(1) % method)//'"') n_global_decomp = ctx % state % cfg % n_cell + 1 ! unreachable; silences -Wuninitialized end select ctx % state % decomp = decompose( & my_rank=my_rank(), & n_ranks=n_ranks(), & n_global=n_global_decomp, & halo_width=ctx % state % halo_width, & is_periodic=ctx % state % is_periodic) call validate_decomp(ctx % state % decomp) call allocate_work_arrays(ctx % state) call resolve_time_scheme(ctx % stepper, ctx % state % cfg % time_scheme) call apply_initial_condition(ctx % state, ctx % state % cfg) ctx % print_freq = ctx % state % cfg % print_freq ctx % t = ctx % state % cfg % time_start ctx % t_comp = 0.0_wp ctx % iter = 0 call timer_reset(ctx % t_total) call timer_reset(ctx % t_io) call timer_reset(ctx % t_iter) ctx % timers_started = .false. ctx % run_complete = .false. ctx % performance_reported = .false. if (len_trim(ctx % state % cfg % restart_file) /= 0) then ! read_checkpoint agrees `ok` across ranks via par_lor, so all ranks take ! the same branch here — no rank can be orphaned at a later collective. call read_checkpoint(ctx % state, ctx % state % cfg % restart_file, ctx % t, ctx % iter, ok, err, & t_comp=ctx % t_comp) if (.not. ok) then if (present(is_ok)) is_ok = .false. if (present(message)) message = trim(err) if (.not. present(is_ok) .and. .not. present(message)) error stop trim(err) return end if end if call log_effective_config(ctx) end subroutine initialize_runtime !> Execute the main time-marching loop until `time_stop` is reached. !! !! Each iteration: optionally recomputes `dt` from the CFL condition; !! clips `dt` at the final time; calls `step`; advances time using !! Kahan compensated summation to minimise floating-point drift; !! writes checkpoints and live snapshots at the configured intervals. !! !! `initialize_runtime` must have been called before `run_solver`. !! On return from `initialize_runtime`, `ctx % t` and `ctx % iter` are !! already set correctly: `ctx % t` equals `time_start` for a fresh run or !! the checkpointed time for a restart; `ctx % iter` is 0 or the restored !! iteration count respectively. !! @param ctx Solver context; `ctx % t` and `ctx % iter` are updated in place. subroutine run_solver(ctx, on_progress, every_steps, every_seconds) type(solver_run_context_t), intent(inout) :: ctx procedure(progress_callback_i), optional :: on_progress integer, intent(in), optional :: every_steps real(wp), intent(in), optional :: every_seconds integer :: steps_taken logical :: finished call run_solver_steps(ctx, huge(steps_taken), steps_taken, finished, & on_progress=on_progress, every_steps=every_steps, & every_seconds=every_seconds) end subroutine run_solver !> Advance the main time loop by at most `max_steps` iterations. !! !! This polling-friendly entry point preserves checkpointing and live-snapshot !! behaviour while allowing adapters to inspect the state between calls. !! Polling is the canonical live-update path for the library ABI; the !! snapshot-file writes remain as a compatibility side channel. !! Lightweight iteration timing (`iter_s`, `elapsed_s`) is always tracked so !! log lines stay meaningful even when `do_timing` is `.false.`; the heavier !! end-of-run summary and fine-grained region timers remain opt-in. !! The loop skeleton itself lives in `run_step_engine` (Q2c); everything 1D !! arrives there through the `ctx1d_*` type-bound hooks below. subroutine run_solver_steps(ctx, max_steps, steps_taken, finished, is_ok, message, & on_progress, every_steps, every_seconds) type(solver_run_context_t), intent(inout) :: ctx integer, intent(in) :: max_steps integer, intent(out) :: steps_taken logical, intent(out) :: finished logical, intent(out), optional :: is_ok character(len=*), intent(out), optional :: message procedure(progress_callback_i), optional :: on_progress integer, intent(in), optional :: every_steps real(wp), intent(in), optional :: every_seconds ! Default status: assume success unless an error is detected below. if (present(is_ok)) is_ok = .true. if (present(message)) message = '' if (ctx % run_complete) then steps_taken = 0 finished = .true. return end if if (.not. associated(ctx % stepper)) then steps_taken = 0 finished = .true. if (present(is_ok)) is_ok = .false. if (present(message)) message = 'solver_runtime: time integrator not initialised' return end if call run_bounded_steps(ctx, max_steps, steps_taken, finished, & ctx % state % cfg % time_stop, ctx % state % cfg % max_iter, & ctx % state % cfg % checkpoint_freq, ctx % print_freq, & ctx % state % cfg % snapshot_freq, max_iter_top_check=.false., & on_progress=on_progress, every_steps=every_steps, every_seconds=every_seconds) end subroutine run_solver_steps !> 1D `raw_dt` hook: recompute `state % dt` from the CFL condition when !! CFL-driven, then return the stored dt (Q2c; body lifted verbatim from !! the pre-engine `run_solver_steps` loop). function ctx1d_raw_dt(ctx) result(dt) class(solver_run_context_t), intent(inout) :: ctx real(wp) :: dt if (ctx % state % cfg % cfl > 0.0_wp) then block real(wp) :: max_ws, dx_min max_ws = par_max_real(compute_max_wave_speed(ctx % state)) if (max_ws <= 0.0_wp) error stop 'solver_runtime: zero max wave speed; CFL dt undefined' ! Method-aware minimum cell width for the CFL stability constraint. ! FVM: minval over interior cells on this rank. On a UNIFORM grid every ! dx_cell is identical, so the rank-local min already equals the global ! min and no collective is needed. On a NON-UNIFORM (file) grid the ! smallest cell may live on another rank, so reduce to the global min ! across ranks (expressed as -max(-x) to reuse par_max_real). FDM: the ! pre-computed global h_min is unchanged. select case (trim(ctx % state % blocks(1) % method)) case (method_fvm) dx_min = minval(ctx % state % mesh % dx_cell(1:ctx % state % decomp % n_local)) if (.not. ctx % state % mesh % uniform) dx_min = -par_max_real(-dx_min) case default ! FDM: use the global minimum node spacing dx_min = ctx % state % mesh % h_min end select ctx % state % dt = ctx % state % cfg % cfl * dx_min / max_ws end block end if dt = ctx % state % dt end function ctx1d_raw_dt !> 1D `set_dt` hook: store the engine's clipped dt. subroutine ctx1d_set_dt(ctx, dt) class(solver_run_context_t), intent(inout) :: ctx real(wp), intent(in) :: dt ctx % state % dt = dt end subroutine ctx1d_set_dt !> 1D `step` hook: advance one step, wrapped in the per-iteration timer. subroutine ctx1d_step(ctx) class(solver_run_context_t), intent(inout) :: ctx call timer_reset(ctx % t_iter) call timer_start(ctx % t_iter) call ctx % stepper(ctx % state) call timer_stop(ctx % t_iter) end subroutine ctx1d_step !> 1D `write_checkpoint_now` hook: write the checkpoint and own the !! failure policy. subroutine ctx1d_checkpoint(ctx) class(solver_run_context_t), intent(inout) :: ctx logical :: ckpt_ok character(len=256) :: ckpt_msg ! par_lor inside write_checkpoint makes ckpt_ok identical on every ! rank, so parallel_fatal fires uniformly (no orphaned ranks). call write_checkpoint(ctx % state, ctx % state % cfg % checkpoint_file, & ctx % t, ctx % iter, ckpt_ok, ckpt_msg, t_comp=ctx % t_comp) if (.not. ckpt_ok) call parallel_fatal(trim(ckpt_msg)) end subroutine ctx1d_checkpoint !> 1D `log_iteration_line` hook: one iteration summary line. subroutine ctx1d_log_line(ctx) class(solver_run_context_t), intent(inout) :: ctx call log_iteration(ctx) end subroutine ctx1d_log_line !> 1D `begin_work` hook: start total wall-time accounting only when real !! work is about to begin, so a zero-step polling call does not consume !! elapsed time. subroutine ctx1d_begin_work(ctx) class(solver_run_context_t), intent(inout) :: ctx if (.not. ctx % timers_started) then call timer_start(ctx % t_total) ctx % timers_started = .true. end if end subroutine ctx1d_begin_work !> 1D `write_snapshot_now` hook: live snapshot side channel. subroutine ctx1d_snapshot(ctx) class(solver_run_context_t), intent(inout) :: ctx call write_snapshot(ctx) end subroutine ctx1d_snapshot !> 1D `finish_work` hook: stop the total timer once the run finishes. subroutine ctx1d_finish_work(ctx) class(solver_run_context_t), intent(inout) :: ctx if (ctx % timers_started) call timer_stop(ctx % t_total) end subroutine ctx1d_finish_work !> 1D `progress_residual` hook: state residual norm for progress ticks. function ctx1d_progress_residual(ctx) result(r) class(solver_run_context_t), intent(inout) :: ctx real(wp) :: r r = ctx % state % resid_glob end function ctx1d_progress_residual !> Write the final solution to the output file and print the performance table. !! !! Output format: four columns (x, ρ, u, p) in `ES20.12` notation, one row !! per interior cell. Must be called after `run_solver`. !! @param ctx Solver context (solution written from `ctx % state % ub`). subroutine finalize_runtime(ctx) type(solver_run_context_t), intent(inout) :: ctx call write_solution_file(ctx, trim(ctx % state % cfg % output_file)) call report_performance_summary(ctx) end subroutine finalize_runtime !> Print the end-of-run performance table once after a completed run. !! !! This wrapper keeps the policy in one place for the legacy runtime driver, !! the session API, and the C ABI: only completed runs emit the heavier !! summary, and only when `do_timing` explicitly enables it. !! @param ctx Solver context (read-only except for the once-only guard flag). subroutine report_performance_summary(ctx) type(solver_run_context_t), intent(inout) :: ctx if (ctx % run_complete .and. ctx % state % cfg % do_timing .and. .not. ctx % performance_reported) then call log_performance_summary(ctx) ctx % performance_reported = .true. end if end subroutine report_performance_summary !> Release all allocations held by `ctx` and finalise the logger. !! Must be the last call in the lifecycle; `ctx` is invalid afterwards. !! @param ctx Solver context to tear down. subroutine teardown_runtime(ctx) type(solver_run_context_t), intent(inout) :: ctx call release_work_arrays(ctx % state) call log_finalize() end subroutine teardown_runtime !> Write the current solution to a text file. !! @param ctx Solver context (solution written from `ctx % state % ub`). !! @param filename Output path to overwrite. subroutine write_solution_file(ctx, filename, is_ok, message) type(solver_run_context_t), intent(inout) :: ctx character(len=*), intent(in) :: filename logical, intent(out), optional :: is_ok character(len=*), intent(out), optional :: message real(wp), allocatable :: local_prim(:, :) integer :: ipt, io_unit, info, n_local real(wp) :: x, rho, u_vel, p logical :: local_ok, failed if (present(is_ok)) is_ok = .true. if (present(message)) message = '' if (ctx % state % cfg % do_timing) call timer_start(ctx % t_io) n_local = ctx % state % n_pt allocate (local_prim(neq, n_local), stat=info) if (info /= 0) error stop 'solver_runtime: write_solution_file local_prim allocation failed' do ipt = 1, n_local rho = ctx % state % ub(1, ipt) ! safe_vel emits 0 velocity (no division) for a vacuum cell (rho <= 0) so ! output cannot raise a SIGFPE; no-op for rho > 0. u_vel = safe_vel(ctx % state % ub(2, ipt), rho) p = (ctx % state % ub(3, ipt) - 0.5_wp * rho * u_vel**2) * & (ctx % state % cfg % gam - 1.0_wp) local_prim(1, ipt) = rho local_prim(2, ipt) = u_vel local_prim(3, ipt) = p end do call gather_solution_to_root(local_prim, n_local, ctx % state % decomp, & ctx % state % gather_buf) deallocate (local_prim, stat=info) local_ok = .true. if (my_rank() == 0) then block character(len=8) :: tokens(max_formats) character(len=8) :: tp_names(4) integer :: n_fmt, kf, ng logical :: fmt_ok, w_ok character(len=256) :: fmt_msg, w_msg character(len=:), allocatable :: fname real(wp), allocatable :: tp_data(:, :) real(wp) :: dx_out !< FVM uniform cell width; unused for FDM. ng = ctx % state % decomp % n_global ! Cell width for the FVM cell-centre output coordinate. ! For FDM, dx_out is never read (the select case below routes to ! x_global), so computing it unconditionally is harmless. dx_out = (ctx % state % cfg % x_right - ctx % state % cfg % x_left) / real(ng, wp) call parse_format_list(ctx % state % cfg % output_format, & [character(len=8) :: 'dat', 'tec', 'plt'], & tokens, n_fmt, fmt_ok, fmt_msg) if (.not. fmt_ok) then call log_warn('write_solution_file: '//trim(fmt_msg)) local_ok = .false. else if (any(tokens(1:n_fmt) == 'tec') .or. any(tokens(1:n_fmt) == 'plt')) then tp_names = [character(len=8) :: 'X', 'rho', 'u', 'p'] allocate (tp_data(4, ng), stat=info) if (info /= 0) error stop 'solver_runtime: write_solution_file tecplot buffer allocation failed' do ipt = 1, ng ! Method-aware x-coordinate: FVM writes cell centres; FDM writes ! nodal coordinates. FDM path is byte-identical to before. select case (trim(ctx % state % blocks(1) % method)) case (method_fvm) ! Non-uniform FVM: write the true cell centres (replicated ! global array); uniform FVM keeps the analytic formula (the ! global array is not built there), bit-for-bit as before. if (allocated(ctx % state % mesh % x_cell_global)) then tp_data(1, ipt) = ctx % state % mesh % x_cell_global(ipt) else tp_data(1, ipt) = ctx % state % cfg % x_left + (real(ipt, wp) - 0.5_wp) * dx_out end if case default tp_data(1, ipt) = ctx % state % mesh % x_global(ipt) end select tp_data(2, ipt) = ctx % state % gather_buf(1, ipt) tp_data(3, ipt) = ctx % state % gather_buf(2, ipt) tp_data(4, ipt) = ctx % state % gather_buf(3, ipt) end do end if do kf = 1, n_fmt fname = derived_filename(trim(filename), tokens(kf), n_fmt) select case (trim(tokens(kf))) case ('tec') call write_tecplot_ascii(fname, 'euler_1d', tp_names, ng, 1, tp_data, w_ok, w_msg) if (.not. w_ok) then call log_warn('write_solution_file: '//trim(w_msg)) local_ok = .false. end if case ('plt') call write_tecplot_binary(fname, 'euler_1d', tp_names, ng, 1, tp_data, w_ok, w_msg) if (.not. w_ok) then call log_warn('write_solution_file: '//trim(w_msg)) local_ok = .false. end if case default ! 'dat' — exact legacy columnar path open (newunit=io_unit, file=fname, status='replace', action='write', iostat=info) if (info /= 0) then local_ok = .false. else do ipt = 1, ng ! Method-aware x-coordinate (same dispatch as Tecplot path). ! FDM path is byte-identical to the pre-Task-5.2 behaviour. select case (trim(ctx % state % blocks(1) % method)) case (method_fvm) if (allocated(ctx % state % mesh % x_cell_global)) then x = ctx % state % mesh % x_cell_global(ipt) else x = ctx % state % cfg % x_left + (real(ipt, wp) - 0.5_wp) * dx_out end if case default x = ctx % state % mesh % x_global(ipt) end select write (io_unit, '(4ES20.12)', iostat=info) x, & ctx % state % gather_buf(1, ipt), ctx % state % gather_buf(2, ipt), & ctx % state % gather_buf(3, ipt) if (info /= 0) then local_ok = .false.; exit end if end do close (io_unit, iostat=info) if (info /= 0) local_ok = .false. end if end select if (.not. local_ok) exit end do if (allocated(tp_data)) deallocate (tp_data, stat=info) end if end block end if ! Collective failure broadcast: every rank must agree on the outcome. failed = par_lor(.not. local_ok) if (failed) then if (present(is_ok)) is_ok = .false. if (present(message)) message = 'solver_runtime: write_solution_file: I/O failed on rank 0' if (.not. present(is_ok) .and. .not. present(message)) & error stop 'solver_runtime: write_solution_file: I/O failed on rank 0' end if if (ctx % state % cfg % do_timing) call timer_stop(ctx % t_io) end subroutine write_solution_file !> Copy the current solution into primitive-variable arrays. !! @param ctx Solver context (read-only). !! @param x Cell-centre coordinates. !! @param rho Density field. !! @param u Velocity field. !! @param p Pressure field. subroutine copy_current_solution(ctx, x, rho, u, p, is_ok, message) type(solver_run_context_t), intent(in) :: ctx real(wp), intent(out) :: x(:), rho(:), u(:), p(:) logical, intent(out), optional :: is_ok character(len=*), intent(out), optional :: message integer :: ipt real(wp) :: u_vel ! The session/C-ABI boundary requires exact-size caller-owned buffers. if (size(x) /= ctx % state % n_pt .or. size(rho) /= ctx % state % n_pt .or. & size(u) /= ctx % state % n_pt .or. size(p) /= ctx % state % n_pt) then if (present(is_ok)) is_ok = .false. if (present(message)) message = 'solver_runtime: copy_current_solution: array size mismatch' return end if if (present(is_ok)) is_ok = .true. if (present(message)) message = '' do ipt = 1, ctx % state % n_pt x(ipt) = ctx % state % mesh % x_node(ipt) rho(ipt) = ctx % state % ub(1, ipt) ! safe_vel emits 0 velocity (no division) for a vacuum cell (rho <= 0) so ! this cannot raise a SIGFPE; no-op for rho > 0. u_vel = safe_vel(ctx % state % ub(2, ipt), rho(ipt)) u(ipt) = u_vel p(ipt) = (ctx % state % ub(3, ipt) - 0.5_wp * rho(ipt) * u_vel**2) * (ctx % state % cfg % gam - 1.0_wp) end do end subroutine copy_current_solution !> Log one iteration summary line (time, residual, per-iteration and elapsed wall time). !! The lightweight `iter_s` and `elapsed_s` fields are always populated. !! Called every `ctx % print_freq` iterations and once unconditionally after the loop. !! @param ctx Solver context (read-only). subroutine log_iteration(ctx) class(solver_run_context_t), intent(in) :: ctx character(len=256) :: msg write (msg, '(A,I6,A,ES12.5,A,ES12.5,A,ES10.3,A,ES10.3)') & 'iter=', ctx % iter, ' t=', ctx % t, ' residual=', ctx % state % resid_glob, & ' iter_s=', timer_elapsed_s(ctx % t_iter), & ' elapsed_s=', timer_elapsed_running_s(ctx % t_total) call log_info(trim(msg)) end subroutine log_iteration !> Print the wall-clock performance table after the run completes. !! Reports total time, I/O time, and the three hot-path timers !! (`compute_resid`, flux precompute, face loop) plus overall throughput. !! Only called when `ctx % state % cfg % do_timing` is `.true.`. !! @param ctx Solver context (read-only). subroutine log_performance_summary(ctx) type(solver_run_context_t), intent(in) :: ctx real(wp) :: elapsed, mcell_steps_per_s character(len=256) :: msg elapsed = real(timer_elapsed_s(ctx % t_total), wp) if (elapsed > 0.0_wp) then mcell_steps_per_s = real(ctx % state % n_pt_global - 1, wp) * real(ctx % iter, wp) / elapsed / 1.0e6_wp else mcell_steps_per_s = 0.0_wp end if call log_info('') call log_info('=== Performance Summary ===') write (msg, '(A,I0,A)') 'grid : ', ctx % state % n_pt_global - 1, ' cells' call log_info(trim(msg)) write (msg, '(A,I0,A,ES10.3,A)') 'steps : ', ctx % iter, ' (t_final =', ctx % t, ' s)' call log_info(trim(msg)) call log_info('-----------------------------------------------') call timer_report(ctx % t_total, 'Total (wall)') call timer_report(ctx % t_io, 'Output I/O') call timer_report(ctx % state % perf % resid, 'compute_resid') call timer_report(ctx % state % perf % fluxsplit, 'flux precompute') call timer_report(ctx % state % perf % faceloop, 'face loop') call log_info('-----------------------------------------------') write (msg, '(A,ES9.3,A)') 'Throughput ', mcell_steps_per_s, ' M cell-steps/s' call log_info(trim(msg)) end subroutine log_performance_summary !> Write a live snapshot of the current solution to `snapshot_file`. !! The file is overwritten on each call. A header comment records the !! current iteration and time; subsequent rows are (x, ρ, u, p) in !! `'(4ES20.12)'` format. Open/write failures are logged as warnings !! rather than fatal errors so that a snapshot glitch does not abort !! the run. !! This file path is retained for compatibility with legacy file-watching !! tooling; new integrations should prefer the polling session/API surface. !! @param ctx Solver context (gather_buf is updated; all other fields are read-only). subroutine write_snapshot(ctx) class(solver_run_context_t), intent(inout) :: ctx real(wp), allocatable :: local_prim(:, :) integer :: u, ipt, io_stat, n_local, ng real(wp) :: x, rho, u_vel, p, dx_snap logical :: local_ok, failed n_local = ctx % state % n_pt ng = ctx % state % decomp % n_global ! Cell width for FVM cell-centre snapshot coordinates. For FDM dx_snap is ! never read (the method dispatch below routes to x_global), so computing ! it unconditionally is harmless — mirrors write_solution_file. dx_snap = (ctx % state % cfg % x_right - ctx % state % cfg % x_left) / real(ng, wp) allocate (local_prim(neq, n_local), stat=io_stat) if (io_stat /= 0) error stop 'solver_runtime: write_snapshot local_prim allocation failed' do ipt = 1, n_local rho = ctx % state % ub(1, ipt) ! safe_vel emits 0 velocity (no division) for a vacuum cell (rho <= 0) so ! output cannot raise a SIGFPE; no-op for rho > 0. u_vel = safe_vel(ctx % state % ub(2, ipt), rho) p = (ctx % state % ub(3, ipt) - 0.5_wp * rho * u_vel**2) * & (ctx % state % cfg % gam - 1.0_wp) local_prim(1, ipt) = rho local_prim(2, ipt) = u_vel local_prim(3, ipt) = p end do call gather_solution_to_root(local_prim, n_local, ctx % state % decomp, & ctx % state % gather_buf) deallocate (local_prim, stat=io_stat) local_ok = .true. if (my_rank() == 0) then open (newunit=u, file=trim(ctx % state % cfg % snapshot_file), & status='replace', action='write', iostat=io_stat) if (io_stat /= 0) then local_ok = .false. else write (u, '(A,I0,A,ES20.12)', iostat=io_stat) '# iter=', ctx % iter, ' t=', ctx % t if (io_stat /= 0) local_ok = .false. if (local_ok) then do ipt = 1, ng ! Method-aware x-coordinate: FVM writes cell centres; FDM writes ! nodal coordinates. FDM path is byte-identical to before. select case (trim(ctx % state % blocks(1) % method)) case (method_fvm) if (allocated(ctx % state % mesh % x_cell_global)) then x = ctx % state % mesh % x_cell_global(ipt) else x = ctx % state % cfg % x_left + (real(ipt, wp) - 0.5_wp) * dx_snap end if case default x = ctx % state % mesh % x_global(ipt) end select write (u, '(4ES20.12)', iostat=io_stat) & x, & ctx % state % gather_buf(1, ipt), & ctx % state % gather_buf(2, ipt), & ctx % state % gather_buf(3, ipt) if (io_stat /= 0) then local_ok = .false. exit end if end do end if close (u, iostat=io_stat) if (io_stat /= 0) local_ok = .false. end if end if ! Every rank participates in the collective; only rank 0 logs the ! warning. Snapshot failures are non-fatal by contract (see header). failed = par_lor(.not. local_ok) if (failed .and. my_rank() == 0) then call log_warn('snapshot: cannot write "'// & trim(ctx % state % cfg % snapshot_file)//'"') end if end subroutine write_snapshot !> Log the effective (post-promotion) configuration at startup. !! Reads from the authoritative copy `ctx % state % cfg`, which may differ !! from the original namelist values when BC promotions have been applied !! (e.g. smooth_wave → periodic, woodward_colella → reflecting). !! @param ctx Solver context (read-only). subroutine log_effective_config(ctx) type(solver_run_context_t), intent(in) :: ctx character(len=256) :: msg associate (cfg => ctx % state % cfg) call log_info('--- Effective Configuration ---') if (trim(cfg % grid_type) == 'file') then write (msg, '(A,A)') 'grid : type=file file=', trim(cfg % grid_file) else write (msg, '(A)') 'grid : type=uniform' end if call log_info(trim(msg)) write (msg, '(A,I0,A,ES11.4,A,ES11.4,A,ES11.4)') & 'grid : n_cell=', cfg % n_cell, ' x=[', cfg % x_left, ', ', cfg % x_right, ']' call log_info(trim(msg)) write (msg, '(A,ES11.4)') ' dx_min=', ctx % state % mesh % h_min call log_info(trim(msg)) write (msg, '(A,A,A,ES11.4,A,ES11.4)') & 'time : scheme=', trim(cfg % time_scheme), ' dt=', cfg % dt, ' cfl=', cfg % cfl call log_info(trim(msg)) write (msg, '(A,ES11.4,A,ES11.4,A)') & ' t=[', cfg % time_start, ', ', cfg % time_stop, ']' call log_info(trim(msg)) write (msg, '(A,A,A,A,A,L1)') & 'schemes : recon=', trim(cfg % recon_scheme), ' flux=', trim(cfg % flux_scheme), & ' char_proj=', ctx % state % use_char_proj call log_info(trim(msg)) write (msg, '(A,A,A,L1)') ' limiter=', trim(cfg % limiter), & ' positivity=', cfg % use_positivity_limiter call log_info(trim(msg)) if (trim(cfg % time_scheme) == 'bdf2') then write (msg, '(A,L1)') ' lapack_solver=', cfg % lapack_solver call log_info(trim(msg)) end if write (msg, '(A,A,A,A,A,A)') & 'ic/bc : problem=', trim(cfg % problem_type), ' bc_left=', trim(cfg % bc_left), & ' bc_right=', trim(cfg % bc_right) call log_info(trim(msg)) select case (trim(cfg % problem_type)) case (problem_sod, problem_lax, problem_woodward_colella, problem_shu_osher, problem_acoustic_pulse) write (msg, '(A,ES11.4,A,ES11.4,A,ES11.4)') & ' rho_L=', cfg % rho_left, ' u_L=', cfg % u_left, ' p_L=', cfg % p_left call log_info(trim(msg)) write (msg, '(A,ES11.4,A,ES11.4,A,ES11.4)') & ' rho_R=', cfg % rho_right, ' u_R=', cfg % u_right, ' p_R=', cfg % p_right call log_info(trim(msg)) case (problem_from_file) call log_info(' ic_file='//trim(cfg % ic_file)) case (problem_udf) call log_info(' ic_udf_src='//trim(cfg % ic_udf_src)) end select if (trim(cfg % bc_left) == 'nonreflecting' .or. trim(cfg % bc_right) == 'nonreflecting') then write (msg, '(A,ES11.4,A,ES11.4,A,ES11.4)') & ' sigma_nrbc=', cfg % sigma_nrbc, ' p_ref_L=', cfg % p_ref_left, & ' p_ref_R=', cfg % p_ref_right call log_info(trim(msg)) end if if (cfg % checkpoint_freq > 0 .or. len_trim(cfg % restart_file) > 0) then write (msg, '(A,I0,A,A)') & 'checkpoint: every ', cfg % checkpoint_freq, ' restart=', trim(cfg % restart_file) call log_info(trim(msg)) end if write (msg, '(A,F8.5)') 'physics : gamma=', cfg % gam call log_info(trim(msg)) call log_info('output : file='//trim(cfg % output_file)) call log_info('-------------------------------') end associate end subroutine log_effective_config end module solver_runtime