!> Parsing and validation of simulation parameters from `case.nml`. !! !! This module handles the ingestion of user-defined parameters using Fortran !! namelists. It is structured to read parameters in logical groups: !! - **`mesh_input`**: Path to the native hexahedral grid files. !! - **`time_input`**: Step counts, timestep size, and CFL control. !! - **`fluid_input`**: Constant/Cantera properties (density, viscosity). !! - **`solver_input`**: Linear solver tolerances and numerical schemes. !! - **`boundary_input`**: Mapping of BC types and values to mesh patches. !! - **`species_input`**: Multi-species transport and chemical mechanism selection. !! - **`output_input`**: Control over VTK and diagnostic file generation. !! - **`profiling_input`**: Control over profiling enable/disable and nested reporting. module mod_input use mod_precision, only : rk, zero, one, name_len, path_len, fatal_error, lowercase implicit none private integer, parameter, public :: max_patches = 64 !< Maximum number of boundary patches supported. integer, parameter, public :: max_species = 256 !< Maximum number of chemical species supported. !> Global container for all parsed simulation parameters. !! !! This structure acts as the source of truth for the entire solver. !! It is populated once at startup and accessed as a read-only object !! by the various solver modules. type, public :: case_params_t !> @name Mesh Configuration character(len=path_len) :: mesh_dir = "mesh_native" !< Directory containing `points.dat`, `cells.dat`, etc. !> @name Time Stepping integer :: nsteps = 0 !< Total number of timesteps to execute. real(rk) :: dt = zero !< Fixed timestep size [s]. real(rk) :: dt_old = zero !< Timestep size used in the previous step [s]. integer :: output_interval = 1 !< Frequency of VTK/PVD output (in steps). integer :: terminal_interval = 1 !< Frequency of lightweight terminal status updates (in steps). character(len=name_len) :: terminal_detail = "health" !< Terminal verbosity: brief, health, or full. logical :: use_dynamic_dt = .false.!< If true, `dt` scales to maintain target `max_cfl`. real(rk) :: max_cfl = 0.5_rk !< Target maximum CFL number for stability. real(rk) :: max_dt = huge(one) !< Absolute timestep ceiling for dynamic dt [s]; huge disables. real(rk) :: min_dt = 1.0e-12_rk !< Absolute timestep floor for automatic dt cuts [s]. real(rk) :: dt_growth_limit = 1.02_rk !< Per-step dynamic-dt growth factor ceiling. logical :: enable_step_rejection = .false. !< If true, reject and retry a timestep when health checks fail. integer :: max_step_retries = 6 !< Maximum rejected attempts for a single timestep. real(rk) :: step_reject_cut_factor = 0.5_rk !< Multiplicative dt cut applied after a rejected trial step. real(rk) :: step_reject_min_dt = 1.0e-12_rk !< Absolute lower bound for rejection retries [s]. real(rk) :: max_KE_growth_per_step = 10.0_rk !< Reject if KE grows by more than this factor; <=0 disables. real(rk) :: ke_reject_floor = 1.0e-12_rk !< Do not apply relative KE-growth rejection when KE before the step is below this floor. real(rk) :: max_cfl_overshoot_factor = 1.25_rk !< Reject if post-step CFL exceeds factor*max_cfl; <=0 disables. logical :: reject_on_pressure_max_iter = .true. !< Reject if the pressure solve reaches pressure_max_iter. logical :: stop_on_min_dt_failure = .true. !< Stop if a rejected step can no longer reduce dt. !> @name Fluid Properties real(rk) :: rho = one !< Constant active flow/projection density [kg/m^3] when `enable_variable_density=.false.`. logical :: enable_variable_density = .false. !< Enable experimental variable-density low-Mach coupling. logical :: enable_density_corrector_projection = .false. !< Post-chemistry pressure-only projection to re-enforce div(u)=S after density change. logical :: upwind_dilatation_density = .false. !< Use upwind density in the advective divergence source calculation. character(len=name_len) :: density_eos = "constant" !< Density source selector; `cantera` uses the selected YAML phase density. real(rk) :: nu = 1.0e-2_rk !< Constant kinematic viscosity [m^2/s]. logical :: enable_variable_nu = .false. !< Allow Cantera-updated flow viscosity/nu; default keeps validation Re fixed. integer :: transport_update_interval = 1 !< Cantera transport-property update interval for mu/D_k only [steps]. !> @name Linear Solver & Numerics integer :: pressure_max_iter = 300 !< Maximum Conjugate Gradient iterations for the Poisson solver. real(rk) :: pressure_tol = 1.0e-10_rk !< Relative residual tolerance for pressure convergence. real(rk) :: pressure_abs_tol = 1.0e-8_rk !< Absolute RMS residual tolerance; PCG exits when either criterion is met. real(rk) :: body_force(3) = zero !< Constant volumetric acceleration vector \((a_x, a_y, a_z)\) [m/s^2]. character(len=name_len) :: convection_scheme = "upwind" !< Legacy momentum advection selector: "upwind" or "central". character(len=name_len) :: momentum_convection_scheme = "" !< Optional momentum-specific scheme; blank falls back to convection_scheme. character(len=name_len) :: species_convection_scheme = "upwind" !< Species advection: upwind, central, bounded_central, bounded_linear, muscl, limited_linear. character(len=name_len) :: energy_convection_scheme = "upwind" !< Enthalpy advection: upwind, central, bounded_central, bounded_linear, muscl, limited_linear. character(len=name_len) :: scalar_limiter = "barth_jespersen" !< Local-bounds limiter for high-order scalar reconstruction. character(len=name_len) :: species_time_scheme = "euler" !< Species transport time scheme: euler or ab2. character(len=name_len) :: energy_time_scheme = "euler" !< Enthalpy transport time scheme: euler or ab2. !> @name Boundary Condition Mapping integer :: n_patches = 0 !< Total number of patches defined in the namelist. character(len=name_len) :: patch_name(max_patches) = "" !< Names identifying mesh patches. character(len=name_len) :: patch_type(max_patches) = "" !< Legacy BC type string. !> Field-specific BC Overrides character(len=name_len) :: patch_velocity_type(max_patches) = "" !< Override BC type for velocity. character(len=name_len) :: patch_pressure_type(max_patches) = "" !< Override BC type for pressure. character(len=name_len) :: patch_temperature_type(max_patches) = "" !< Override BC type for temperature/enthalpy. real(rk) :: patch_u(max_patches) = zero !< Specified x-velocity on patch [m/s]. real(rk) :: patch_v(max_patches) = zero !< Specified y-velocity on patch [m/s]. real(rk) :: patch_w(max_patches) = zero !< Specified z-velocity on patch [m/s]. !< Mass flux input [kg/m^2/s]. Signed for mass_flux; !! positive magnitude for inlet/outlet_mass_flux. real(rk) :: patch_mass_flux(max_patches) = zero real(rk) :: patch_p(max_patches) = zero !< Specified projection-pressure boundary value [Pa], not Cantera \(p_0\). real(rk) :: patch_dpdn(max_patches) = zero !< Specified pressure gradient on patch [Pa/m]. real(rk) :: patch_T(max_patches) = 300.0_rk !< Specified temperature on patch [K]. !> Species Boundary Conditions character(len=name_len) :: patch_species_type(max_patches) = "" !< Override BC type for mass fractions. real(rk) :: patch_Y(max_species, max_patches) = zero !< Legacy index-based \(Y_k\) mass fractions on patch. character(len=1024) :: patch_composition(max_patches) = "" !< Preferred named patch mixture, e.g. "CH4:1,O2:0.233,N2:0.767". !> @name Data Output character(len=path_len) :: output_dir = "output" !< Directory for result storage. logical :: write_vtu = .true. !< If true, generates Unstructured VTK files. logical :: write_diagnostics = .true. !< If true, writes global residuals to `diagnostics.csv`. character(len=name_len) :: vtu_format = "ascii" !< Output format: "ascii" or "binary". !> @name Restart Controls logical :: restart_from_file = .false. !< If true, initialize fields from restart_file after allocation. character(len=path_len) :: restart_file = "" !< Rank-independent restart file to read. logical :: write_restart = .false. !< If true, periodically write global restart files. integer :: restart_interval = 0 !< Restart write cadence [steps]; required when write_restart=.true. character(len=path_len) :: restart_output_dir = "restart" !< Restart directory, relative to output_dir unless absolute. integer :: restart_keep = 0 !< Number of newest restart files to keep; 0 keeps all. logical :: restart_use_file_dt = .false. !< If true, restore dt from restart; otherwise namelist dt wins. logical :: restart_nsteps_are_additional = .false. !< If true, add namelist nsteps to restart step. !> @name Multi-Species Transport logical :: enable_species = .false. !< Enable advection-diffusion of mass fractions. logical :: enable_reactions = .false. !< Reserved for chemical source terms; supported paths require non-reacting transport. integer :: nspecies = 0 !< Total number of transport species. character(len=name_len) :: species_name(max_species) = "" !< List of species names. real(rk) :: species_diffusivity(max_species) = 0.0_rk !< Constant species diffusivity \(D_k\) [m^2/s]. real(rk) :: initial_Y(max_species) = 0.0_rk !< Legacy index/name-matched global initial mass fractions. character(len=1024) :: initial_composition = "" !< Preferred named initial mixture, e.g. "CH4:0.028,O2:0.225,N2:0.747". !> @name Operator-split Cantera chemistry controls integer :: chemistry_update_interval = 1 !< Chemistry source update cadence [steps]. real(rk) :: chemistry_rtol = 1.0e-8_rk !< Relative tolerance for Cantera ReactorNet chemistry integration. real(rk) :: chemistry_atol = 1.0e-15_rk !< Absolute tolerance for Cantera ReactorNet chemistry integration. integer :: chemistry_max_steps = 10000 !< Maximum internal ReactorNet steps per cell chemistry advance. logical :: chemistry_energy_enabled = .true. !< If true, Cantera chemistry solves the reactor energy equation. real(rk) :: chemistry_temperature_cutoff = zero !< Skip chemistry below this temperature [K]; zero disables the cutoff. real(rk) :: chemistry_min_reactive_mass_fraction = zero !< Legacy skip if max(Y_k) is below this threshold; zero disables. integer :: chemistry_n_active_species = 0 !< Number of named species used for optional chemistry activity screening. character(len=name_len) :: chemistry_active_species_name(max_species) = "" !< Species names used for optional chemistry screening. real(rk) :: chemistry_active_species_threshold = zero !< Skip chemistry if all named active-species Y are below this threshold. real(rk) :: chemistry_max_dT_per_step = zero !< Optional chemistry dt/source limiter threshold [K]; zero disables. real(rk) :: chemistry_max_rel_rho_change_per_step = zero !< Optional chemistry density-change threshold; zero disables. real(rk) :: chemistry_max_dY_per_step = zero !< Optional chemistry species-change threshold; zero disables. real(rk) :: chemistry_dt_safety = 0.5_rk !< Safety factor for next-step chemistry dt cuts. real(rk) :: chemistry_dt_min_factor = 0.1_rk !< Minimum multiplicative dt cut in one chemistry-control event. logical :: chemistry_limit_source_update = .false. !< If true, under-relax/clip committed chemistry source update. real(rk) :: chemistry_source_relaxation = one !< Max fraction of Cantera chemistry update to commit when source limiting. logical :: enable_chemistry_load_balancing = .false. !< If true, distribute active chemistry cells globally round-robin. !> @name Chemistry and Reacting Flow Stability Upgrades logical :: enable_chemistry_subcycling = .false. !< If true, subcycle Cantera chemistry integration. integer :: max_chemistry_subcycles = 8 !< Upper limit on chemistry subcycles per timestep. real(rk) :: subcycling_dT_threshold = 10.0_rk !< Local temperature change threshold [K] to trigger subcycling. logical :: reject_on_nan = .true. !< If true, reject trial steps containing NaN/Inf values. real(rk) :: min_density_ratio_limit = 0.05_rk !< Minimum density ratio relative to background before step rejection. real(rk) :: max_temperature_limit = 3200.0_rk !< Maximum temperature limit [K] before step rejection. real(rk) :: min_temperature_limit = 200.0_rk !< Minimum temperature limit [K] before step rejection. real(rk) :: max_pressure_residual_limit = 1.0e-2_rk !< Maximum allowed PCG pressure residual before step rejection. !> @name Ignition-kernel initialization logical :: enable_ignition_kernel = .false. !< If true, overwrite an initial hot spherical kernel before the first step. character(len=name_len) :: ignition_kernel_shape = "sphere" !< First supported shape is sphere. real(rk) :: ignition_kernel_center(3) = zero !< Kernel center [m]. real(rk) :: ignition_kernel_radius = zero !< Kernel radius [m]. real(rk) :: ignition_kernel_T = zero !< Kernel temperature [K]. character(len=1024) :: ignition_kernel_composition = "" !< Optional named composition inside kernel. character(len=name_len) :: ignition_kernel_blend = "overwrite" !< First supported blend mode is overwrite. !> @name Internal Registry (Reacting discovery) integer :: namelist_nspecies = 0 !< Number of species specified in `case.nml`. character(len=name_len) :: namelist_species_name(max_species) = "" !< Names specified in `case.nml`. !> @name Cantera Bridge Integration logical :: enable_cantera_fluid = .false. !< Use Cantera for mixture-averaged viscosity; flow density remains `rho`. logical :: enable_cantera_species = .false. !< Use Cantera for mixture-averaged \(D_k\). character(len=path_len) :: cantera_mech_file = "gri30.yaml" !< Path to YAML/CTI mechanism file. character(len=name_len) :: cantera_phase_name = "" !< Optional Cantera phase name; blank uses first/default phase. integer :: cantera_reaction_count = 0 !< Number of reactions reported by the selected Cantera phase. real(rk) :: background_temp = 300.0_rk !< Fixed temperature for property evaluation [K]. real(rk) :: background_press = 101325.0_rk !< Uniform thermodynamic pressure \(p_0\) passed to Cantera [Pa]. !> @name Enthalpy Energy Equation Controls logical :: enable_energy = .false. !< Enable enthalpy/temperature field storage. logical :: enable_cantera_thermo = .false. !< Use Cantera for sensible h(T,Y,p0) and T(h,Y,p0). logical :: enable_species_enthalpy_diffusion = .false. !< Include optional species-enthalpy diffusion correction in the \(h\) equation. integer :: thermo_update_interval = 1 !< Reserved Cantera thermo interval [steps]; must remain 1. character(len=name_len) :: thermo_default_species = 'N2' !< Single-species Cantera thermo fallback when species transport is off. real(rk) :: initial_T = 300.0_rk !< Initial gas temperature [K]. real(rk) :: energy_reference_T = 298.15_rk !< Reference temperature for sensible enthalpy [K]. real(rk) :: energy_reference_h = zero !< Reference sensible enthalpy for constant-cp mode [J/kg]. real(rk) :: energy_cp = 1005.0_rk !< Constant heat capacity for non-Cantera thermo [J/kg/K]. real(rk) :: energy_lambda = 2.6e-2_rk !< Constant thermal conductivity for non-Cantera thermo [W/m/K]. !> @name Spectral radiation coupling controls logical :: enable_radiation = .false. !< Enable radiation-source updates on the radiation communicator. integer :: radiation_update_interval = 1 !< Recompute qrad every N flow steps; otherwise reuse last qrad. integer :: radiation_n_wavenumbers = 0 !< Number of spectral/wavenumber tasks decomposed over radiation MPI ranks. integer :: radiation_n_species = 0 !< Number of transported species used by the radiation model. character(len=name_len) :: radiation_species_name(max_species) = "" !< Radiation-relevant transported species names. character(len=name_len) :: radiation_source_model = "none" !< Radiation source model: none or spectral_test. real(rk) :: radiation_source_scale = 0.0_rk !< Test-model source scale [W/m^3]. integer :: radiation_n_scalars = 0 !< Reserved number of generic scalar arrays for radiation. character(len=name_len) :: radiation_scalar_name(max_species) = "" !< Reserved selected scalar names for future radiation models. character(len=name_len) :: radiation_pressure_source = "background" !< Pressure source for radiation state: background or system. logical :: radiation_debug = .false. !< Write extra gather/checksum radiation diagnostics. logical :: write_radiation_diagnostics = .true. !< Write radiation_setup.txt and radiation_diagnostics.csv. real(rk) :: radiation_absorption_coeff = 1.0_rk !< Constant grey gas absorption coefficient [1/m]. real(rk) :: radiation_emissivity = 1.0_rk !< Wall boundary surface emissivity [0-1]. character(len=name_len) :: radiation_dom_quadrature = "s2" !< DOM angular quadrature type: s2 or s4. !> @name Profiling Controls logical :: enable_profiling = .true. !< Enable wall-clock profiling. logical :: nested_profiling = .true. !< If true, print nested profiling tree. logical :: write_cantera_cache_stats = .false. !< If true, print final Cantera cache hit/miss statistics. logical :: write_chemistry_screening_stats = .false. !< If true, print final chemistry screening/skipped-cell statistics. logical :: variable_density_debug = .false. !< If true, print verbose experimental variable-density diagnostics. integer :: variable_density_debug_interval = 1 !< Print variable-density debug every N energy steps when enabled. end type case_params_t public :: read_case_params contains !> Orchestrates the reading of all namelist blocks from the configuration file. !! !! Performs a sequential read of all expected namelist groups and !! triggers a final validation pass to ensure physical consistency. !! !! @param filename Path to the `.nml` file (usually `case.nml`). !! @param params Container to be populated with parsed values. subroutine read_case_params(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params call read_mesh_input(filename, params) call read_time_input(filename, params) call read_fluid_input(filename, params) call read_solver_input(filename, params) call read_boundary_input(filename, params) call read_species_input(filename, params) call read_chemistry_input(filename, params) call read_energy_input(filename, params) call read_ignition_input(filename, params) call read_radiation_input(filename, params) call read_output_input(filename, params) call read_restart_input(filename, params) call read_profiling_input(filename, params) call validate_params(params) ! Density semantics: ! constant-density mode keeps transport%rho=params%rho and treats ! energy%rho_thermo as a diagnostic; variable-density density_eos="cantera" ! requires Cantera thermo and syncs transport%rho from energy%rho_thermo. ! The selected YAML phase defines the actual EOS. "ideal_gas" is parsed ! only to reserve the name and is rejected for active variable-density use. params%density_eos = trim(lowercase(params%density_eos)) if (len_trim(params%density_eos) == 0) params%density_eos = "constant" select case (trim(params%density_eos)) case ("constant", "cantera", "ideal_gas") continue case default call fatal_error("input", "density_eos must be one of: constant, cantera, ideal_gas") end select if (params%enable_variable_density) then select case (trim(params%density_eos)) case ("cantera") if (.not. params%enable_energy) then call fatal_error("input", "enable_variable_density with density_eos=cantera requires enable_energy=.true.") end if if (.not. params%enable_cantera_thermo) then call fatal_error("input", "enable_variable_density with density_eos=cantera requires enable_cantera_thermo=.true.") end if if (params%thermo_update_interval /= 1) then call fatal_error("input", "enable_variable_density requires thermo_update_interval=1") end if case ("ideal_gas") call fatal_error("input", "density_eos=ideal_gas is parsed but not implemented for variable-density mode yet") case default call fatal_error("input", "enable_variable_density currently requires density_eos=cantera") end select end if end subroutine read_case_params !> Reads the `&mesh_input` namelist block. subroutine read_mesh_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params character(len=path_len) :: mesh_dir integer :: unit_id, ios namelist /mesh_input/ mesh_dir mesh_dir = params%mesh_dir call open_namelist_file(filename, unit_id, ios) if (ios == 0) then read(unit_id, nml=mesh_input, iostat=ios) close(unit_id) end if if (ios > 0) call fatal_error('input', 'failed reading &mesh_input') params%mesh_dir = mesh_dir end subroutine read_mesh_input !> Reads the `&time_input` namelist block. subroutine read_time_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params integer :: nsteps, output_interval real(rk) :: dt, max_cfl, max_dt, min_dt, dt_growth_limit logical :: use_dynamic_dt logical :: enable_step_rejection, reject_on_pressure_max_iter, stop_on_min_dt_failure integer :: max_step_retries real(rk) :: step_reject_cut_factor, step_reject_min_dt real(rk) :: max_KE_growth_per_step, ke_reject_floor, max_cfl_overshoot_factor integer :: unit_id, ios namelist /time_input/ nsteps, dt, output_interval, use_dynamic_dt, max_cfl, & max_dt, min_dt, dt_growth_limit, & enable_step_rejection, max_step_retries, step_reject_cut_factor, & step_reject_min_dt, max_KE_growth_per_step, ke_reject_floor, max_cfl_overshoot_factor, & reject_on_pressure_max_iter, stop_on_min_dt_failure nsteps = params%nsteps dt = params%dt output_interval = params%output_interval use_dynamic_dt = params%use_dynamic_dt max_cfl = params%max_cfl max_dt = params%max_dt min_dt = params%min_dt dt_growth_limit = params%dt_growth_limit enable_step_rejection = params%enable_step_rejection max_step_retries = params%max_step_retries step_reject_cut_factor = params%step_reject_cut_factor step_reject_min_dt = params%step_reject_min_dt max_KE_growth_per_step = params%max_KE_growth_per_step ke_reject_floor = params%ke_reject_floor max_cfl_overshoot_factor = params%max_cfl_overshoot_factor reject_on_pressure_max_iter = params%reject_on_pressure_max_iter stop_on_min_dt_failure = params%stop_on_min_dt_failure call open_namelist_file(filename, unit_id, ios) if (ios == 0) then read(unit_id, nml=time_input, iostat=ios) close(unit_id) end if if (ios > 0) call fatal_error('input', 'failed reading &time_input') params%nsteps = nsteps params%dt = dt params%output_interval = output_interval params%use_dynamic_dt = use_dynamic_dt params%max_cfl = max_cfl params%max_dt = max_dt params%min_dt = min_dt params%dt_growth_limit = dt_growth_limit params%enable_step_rejection = enable_step_rejection params%max_step_retries = max_step_retries params%step_reject_cut_factor = step_reject_cut_factor params%step_reject_min_dt = step_reject_min_dt params%max_KE_growth_per_step = max_KE_growth_per_step params%ke_reject_floor = ke_reject_floor params%max_cfl_overshoot_factor = max_cfl_overshoot_factor params%reject_on_pressure_max_iter = reject_on_pressure_max_iter params%stop_on_min_dt_failure = stop_on_min_dt_failure end subroutine read_time_input !> Reads the `&fluid_input` namelist block. subroutine read_fluid_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params real(rk) :: rho, nu logical :: enable_cantera character(len=path_len) :: cantera_mech_file character(len=name_len) :: cantera_phase_name real(rk) :: background_temp real(rk) :: background_press integer :: transport_update_interval integer :: unit_id, ios logical :: enable_variable_density logical :: enable_density_corrector_projection logical :: upwind_dilatation_density character(len=name_len) :: density_eos = "constant" !< Density EOS selector: constant, cantera, ideal_gas. logical :: enable_variable_nu namelist /fluid_input/ rho, nu, enable_cantera, cantera_mech_file, cantera_phase_name, & background_temp, background_press, transport_update_interval, & enable_variable_density, enable_density_corrector_projection, & density_eos, enable_variable_nu, upwind_dilatation_density rho = params%rho nu = params%nu enable_cantera = params%enable_cantera_fluid cantera_mech_file = params%cantera_mech_file cantera_phase_name = params%cantera_phase_name background_temp = params%background_temp background_press = params%background_press transport_update_interval = params%transport_update_interval enable_variable_density = params%enable_variable_density enable_density_corrector_projection = params%enable_density_corrector_projection density_eos = params%density_eos enable_variable_nu = params%enable_variable_nu upwind_dilatation_density = params%upwind_dilatation_density call open_namelist_file(filename, unit_id, ios) if (ios == 0) then read(unit_id, nml=fluid_input, iostat=ios) close(unit_id) end if if (ios > 0) call fatal_error('input', 'failed reading &fluid_input') params%rho = rho params%nu = nu params%enable_cantera_fluid = enable_cantera params%enable_variable_density = enable_variable_density params%enable_density_corrector_projection = enable_density_corrector_projection params%density_eos = trim(density_eos) params%enable_variable_nu = enable_variable_nu params%upwind_dilatation_density = upwind_dilatation_density params%cantera_mech_file = cantera_mech_file params%cantera_phase_name = trim(cantera_phase_name) params%background_temp = background_temp params%background_press = background_press params%transport_update_interval = transport_update_interval end subroutine read_fluid_input !> Reads the `&solver_input` namelist block. subroutine read_solver_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params integer :: pressure_max_iter real(rk) :: pressure_tol real(rk) :: pressure_abs_tol real(rk) :: body_force_x, body_force_y, body_force_z character(len=name_len) :: convection_scheme character(len=name_len) :: momentum_convection_scheme character(len=name_len) :: species_convection_scheme character(len=name_len) :: energy_convection_scheme character(len=name_len) :: scalar_limiter character(len=name_len) :: species_time_scheme character(len=name_len) :: energy_time_scheme integer :: unit_id, ios namelist /solver_input/ pressure_max_iter, pressure_tol, pressure_abs_tol, & body_force_x, body_force_y, body_force_z, & convection_scheme, momentum_convection_scheme, & species_convection_scheme, energy_convection_scheme, & scalar_limiter, species_time_scheme, energy_time_scheme pressure_max_iter = params%pressure_max_iter pressure_tol = params%pressure_tol pressure_abs_tol = params%pressure_abs_tol body_force_x = params%body_force(1) body_force_y = params%body_force(2) body_force_z = params%body_force(3) convection_scheme = params%convection_scheme momentum_convection_scheme = params%momentum_convection_scheme species_convection_scheme = params%species_convection_scheme energy_convection_scheme = params%energy_convection_scheme scalar_limiter = params%scalar_limiter species_time_scheme = params%species_time_scheme energy_time_scheme = params%energy_time_scheme call open_namelist_file(filename, unit_id, ios) if (ios == 0) then read(unit_id, nml=solver_input, iostat=ios) close(unit_id) end if if (ios > 0) call fatal_error('input', 'failed reading &solver_input') params%pressure_max_iter = pressure_max_iter params%pressure_tol = pressure_tol params%pressure_abs_tol = pressure_abs_tol params%body_force = [body_force_x, body_force_y, body_force_z] params%convection_scheme = trim(lowercase(convection_scheme)) params%momentum_convection_scheme = trim(lowercase(momentum_convection_scheme)) if (len_trim(params%momentum_convection_scheme) == 0) then params%momentum_convection_scheme = params%convection_scheme end if ! The legacy field remains the value consumed by the existing projection solver. params%convection_scheme = params%momentum_convection_scheme params%species_convection_scheme = trim(lowercase(species_convection_scheme)) params%energy_convection_scheme = trim(lowercase(energy_convection_scheme)) params%scalar_limiter = trim(lowercase(scalar_limiter)) params%species_time_scheme = trim(lowercase(species_time_scheme)) params%energy_time_scheme = trim(lowercase(energy_time_scheme)) end subroutine read_solver_input !> Reads the `&boundary_input` namelist block. subroutine read_boundary_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params integer :: n_patches character(len=name_len) :: patch_name(max_patches) character(len=name_len) :: patch_type(max_patches) character(len=name_len) :: patch_velocity_type(max_patches) character(len=name_len) :: patch_pressure_type(max_patches) character(len=name_len) :: patch_temperature_type(max_patches) character(len=name_len) :: patch_species_type(max_patches) real(rk) :: patch_u(max_patches) real(rk) :: patch_v(max_patches) real(rk) :: patch_w(max_patches) real(rk) :: patch_mass_flux(max_patches) real(rk) :: patch_p(max_patches) real(rk) :: patch_dpdn(max_patches) real(rk) :: patch_T(max_patches) real(rk), save :: patch_Y(max_species, max_patches) = zero character(len=1024), save :: patch_composition(max_patches) = "" integer :: unit_id, ios, i namelist /boundary_input/ n_patches, patch_name, patch_type, & patch_velocity_type, patch_pressure_type, & patch_temperature_type, patch_species_type, & patch_u, patch_v, patch_w, patch_mass_flux, & patch_p, patch_dpdn, patch_T, patch_Y, patch_composition n_patches = params%n_patches patch_name = params%patch_name patch_type = params%patch_type patch_velocity_type = params%patch_velocity_type patch_pressure_type = params%patch_pressure_type patch_temperature_type = params%patch_temperature_type patch_species_type = params%patch_species_type patch_u = params%patch_u patch_v = params%patch_v patch_w = params%patch_w patch_mass_flux = params%patch_mass_flux patch_p = params%patch_p patch_dpdn = params%patch_dpdn patch_T = params%patch_T patch_Y = params%patch_Y patch_composition = params%patch_composition call open_namelist_file(filename, unit_id, ios) if (ios == 0) then read(unit_id, nml=boundary_input, iostat=ios) close(unit_id) end if if (ios > 0) call fatal_error('input', 'failed reading &boundary_input') params%n_patches = n_patches params%patch_name = patch_name params%patch_type = patch_type params%patch_velocity_type = patch_velocity_type params%patch_pressure_type = patch_pressure_type params%patch_temperature_type = patch_temperature_type params%patch_species_type = patch_species_type do i = 1, max_patches params%patch_type(i) = trim(lowercase(params%patch_type(i))) params%patch_velocity_type(i) = trim(lowercase(params%patch_velocity_type(i))) params%patch_pressure_type(i) = trim(lowercase(params%patch_pressure_type(i))) params%patch_temperature_type(i) = trim(lowercase(params%patch_temperature_type(i))) params%patch_species_type(i) = trim(lowercase(params%patch_species_type(i))) end do params%patch_u = patch_u params%patch_v = patch_v params%patch_w = patch_w params%patch_mass_flux = patch_mass_flux params%patch_p = patch_p params%patch_dpdn = patch_dpdn params%patch_T = patch_T params%patch_Y = patch_Y params%patch_composition = patch_composition end subroutine read_boundary_input !> Reads the `&profiling_input` block. !! !! This parser is intentionally explicit because profiling is optional and !! should be robust to block ordering in case.nml. subroutine read_profiling_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params character(len=512) :: line character(len=256) :: key, value integer :: unit_id, ios, eqpos, comment_pos, comma_pos logical :: in_block, found_block in_block = .false. found_block = .false. open(newunit=unit_id, file=trim(filename), status='old', action='read', iostat=ios) if (ios /= 0) return do read(unit_id, '(a)', iostat=ios) line if (ios /= 0) exit comment_pos = index(line, '!') if (comment_pos > 0) line = line(:comment_pos-1) line = adjustl(line) if (len_trim(line) == 0) cycle if (.not. in_block) then if (index(lowercase(line), '&profiling_input') == 1) then in_block = .true. found_block = .true. end if cycle end if if (index(line, '/') > 0) exit eqpos = index(line, '=') if (eqpos <= 0) cycle key = trim(adjustl(lowercase(line(:eqpos-1)))) value = trim(adjustl(lowercase(line(eqpos+1:)))) comma_pos = index(value, ',') if (comma_pos > 0) value = value(:comma_pos-1) value = trim(value) select case (trim(key)) case ('enable_profiling') call parse_logical_value(value, params%enable_profiling, 'enable_profiling') case ('nested_profiling') call parse_logical_value(value, params%nested_profiling, 'nested_profiling') case ('write_cantera_cache_stats') call parse_logical_value(value, params%write_cantera_cache_stats, 'write_cantera_cache_stats') case ('write_chemistry_screening_stats') call parse_logical_value(value, params%write_chemistry_screening_stats, 'write_chemistry_screening_stats') case ('variable_density_debug') call parse_logical_value(value, params%variable_density_debug, 'variable_density_debug') case ('variable_density_debug_interval') read(value, *, iostat=ios) params%variable_density_debug_interval if (ios /= 0) call fatal_error('input', 'invalid integer for variable_density_debug_interval: '//trim(value)) case default call fatal_error('input', 'unknown variable in &profiling_input: '//trim(key)) end select end do close(unit_id) if (found_block .and. in_block .and. ios /= 0) then call fatal_error('input', 'unterminated &profiling_input block') end if contains !> Parses a string representation of a logical value into a Fortran logical. !! !! Supports common formats like '.true.', 'true', 't', and their false counterparts. !! !! @param value_text The string to parse. !! @param output_value The resulting logical value. !! @param field_name The name of the field being parsed (used for error reporting). subroutine parse_logical_value(value_text, output_value, field_name) character(len=*), intent(in) :: value_text logical, intent(out) :: output_value character(len=*), intent(in) :: field_name select case (trim(value_text)) case ('.true.', 'true', 't', '.t.') output_value = .true. case ('.false.', 'false', 'f', '.f.') output_value = .false. case default call fatal_error('input', 'invalid logical for '//trim(field_name)//': '//trim(value_text)) end select end subroutine parse_logical_value end subroutine read_profiling_input !> Reads the optional `&chemistry_input` namelist block. !! !! Chemistry is operator-split from advection/diffusion and is solved by the !! Cantera C++ bridge as a cell-local constant-pressure reactor. The block !! is optional; defaults are conservative and do not enable chemistry by !! themselves. `enable_reactions` remains in `&species_input` for backward !! compatibility with existing cases. subroutine read_chemistry_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params integer :: chemistry_update_interval real(rk) :: chemistry_rtol, chemistry_atol integer :: chemistry_max_steps logical :: chemistry_energy_enabled real(rk) :: chemistry_temperature_cutoff real(rk) :: chemistry_min_reactive_mass_fraction integer :: chemistry_n_active_species character(len=name_len) :: chemistry_active_species_name(max_species) real(rk) :: chemistry_active_species_threshold real(rk) :: chemistry_max_dT_per_step, chemistry_max_rel_rho_change_per_step real(rk) :: chemistry_max_dY_per_step, chemistry_dt_safety, chemistry_dt_min_factor logical :: chemistry_limit_source_update real(rk) :: chemistry_source_relaxation logical :: enable_chemistry_subcycling integer :: max_chemistry_subcycles real(rk) :: subcycling_dT_threshold logical :: reject_on_nan real(rk) :: min_density_ratio_limit real(rk) :: max_temperature_limit real(rk) :: min_temperature_limit real(rk) :: max_pressure_residual_limit logical :: enable_chemistry_load_balancing integer :: unit_id, ios, k namelist /chemistry_input/ chemistry_update_interval, chemistry_rtol, chemistry_atol, & chemistry_max_steps, chemistry_energy_enabled, & chemistry_temperature_cutoff, chemistry_min_reactive_mass_fraction, & chemistry_n_active_species, chemistry_active_species_name, & chemistry_active_species_threshold, & chemistry_max_dT_per_step, chemistry_max_rel_rho_change_per_step, & chemistry_max_dY_per_step, chemistry_dt_safety, chemistry_dt_min_factor, & chemistry_limit_source_update, chemistry_source_relaxation, & enable_chemistry_subcycling, max_chemistry_subcycles, subcycling_dT_threshold, & reject_on_nan, min_density_ratio_limit, max_temperature_limit, & min_temperature_limit, max_pressure_residual_limit, & enable_chemistry_load_balancing chemistry_update_interval = params%chemistry_update_interval chemistry_rtol = params%chemistry_rtol chemistry_atol = params%chemistry_atol chemistry_max_steps = params%chemistry_max_steps chemistry_energy_enabled = params%chemistry_energy_enabled chemistry_temperature_cutoff = params%chemistry_temperature_cutoff chemistry_min_reactive_mass_fraction = params%chemistry_min_reactive_mass_fraction chemistry_n_active_species = params%chemistry_n_active_species chemistry_active_species_name = params%chemistry_active_species_name chemistry_active_species_threshold = params%chemistry_active_species_threshold chemistry_max_dT_per_step = params%chemistry_max_dT_per_step chemistry_max_rel_rho_change_per_step = params%chemistry_max_rel_rho_change_per_step chemistry_max_dY_per_step = params%chemistry_max_dY_per_step chemistry_dt_safety = params%chemistry_dt_safety chemistry_dt_min_factor = params%chemistry_dt_min_factor chemistry_limit_source_update = params%chemistry_limit_source_update chemistry_source_relaxation = params%chemistry_source_relaxation enable_chemistry_subcycling = params%enable_chemistry_subcycling max_chemistry_subcycles = params%max_chemistry_subcycles subcycling_dT_threshold = params%subcycling_dT_threshold reject_on_nan = params%reject_on_nan min_density_ratio_limit = params%min_density_ratio_limit max_temperature_limit = params%max_temperature_limit min_temperature_limit = params%min_temperature_limit max_pressure_residual_limit = params%max_pressure_residual_limit enable_chemistry_load_balancing = params%enable_chemistry_load_balancing call open_namelist_file(filename, unit_id, ios) if (ios == 0) then read(unit_id, nml=chemistry_input, iostat=ios) close(unit_id) end if if (ios /= 0 .and. ios /= -1) then call fatal_error('input', 'failed reading &chemistry_input. Check for unknown variables or typos.') end if if (ios == 0) then params%chemistry_update_interval = chemistry_update_interval params%chemistry_rtol = chemistry_rtol params%chemistry_atol = chemistry_atol params%chemistry_max_steps = chemistry_max_steps params%chemistry_energy_enabled = chemistry_energy_enabled params%chemistry_temperature_cutoff = chemistry_temperature_cutoff params%chemistry_min_reactive_mass_fraction = chemistry_min_reactive_mass_fraction params%chemistry_n_active_species = chemistry_n_active_species params%chemistry_active_species_name = chemistry_active_species_name params%chemistry_active_species_threshold = chemistry_active_species_threshold params%chemistry_max_dT_per_step = chemistry_max_dT_per_step params%chemistry_max_rel_rho_change_per_step = chemistry_max_rel_rho_change_per_step params%chemistry_max_dY_per_step = chemistry_max_dY_per_step params%chemistry_dt_safety = chemistry_dt_safety params%chemistry_dt_min_factor = chemistry_dt_min_factor params%chemistry_limit_source_update = chemistry_limit_source_update params%chemistry_source_relaxation = chemistry_source_relaxation params%enable_chemistry_subcycling = enable_chemistry_subcycling params%max_chemistry_subcycles = max_chemistry_subcycles params%subcycling_dT_threshold = subcycling_dT_threshold params%reject_on_nan = reject_on_nan params%min_density_ratio_limit = min_density_ratio_limit params%max_temperature_limit = max_temperature_limit params%min_temperature_limit = min_temperature_limit params%max_pressure_residual_limit = max_pressure_residual_limit params%enable_chemistry_load_balancing = enable_chemistry_load_balancing if (params%chemistry_n_active_species <= 0) then do k = 1, max_species if (len_trim(params%chemistry_active_species_name(k)) > 0) then params%chemistry_n_active_species = k end if end do end if end if end subroutine read_chemistry_input !> Reads the optional `&radiation_input` namelist block. !! !! Radiation is intentionally configured independently from flow MPI. Flow !! remains spatially decomposed while radiation decomposes spectral tasks !! (`radiation_n_wavenumbers`) over the radiation communicator and evaluates !! source terms on the full replicated mesh. subroutine read_radiation_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params logical :: enable_radiation integer :: radiation_update_interval, radiation_n_wavenumbers, radiation_n_species integer :: radiation_n_scalars character(len=name_len), save :: radiation_species_name(max_species) = "" character(len=name_len), save :: radiation_scalar_name(max_species) = "" character(len=name_len) :: radiation_source_model, radiation_pressure_source, radiation_dom_quadrature real(rk) :: radiation_source_scale logical :: radiation_debug, write_radiation_diagnostics real(rk) :: radiation_absorption_coeff, radiation_emissivity integer :: unit_id, ios, i namelist /radiation_input/ enable_radiation, radiation_update_interval, & radiation_n_wavenumbers, radiation_n_species, & radiation_species_name, radiation_source_model, & radiation_source_scale, radiation_n_scalars, & radiation_scalar_name, radiation_pressure_source, & radiation_debug, write_radiation_diagnostics, & radiation_absorption_coeff, radiation_emissivity, & radiation_dom_quadrature enable_radiation = params%enable_radiation radiation_update_interval = params%radiation_update_interval radiation_n_wavenumbers = params%radiation_n_wavenumbers radiation_n_species = params%radiation_n_species radiation_species_name = params%radiation_species_name radiation_scalar_name = params%radiation_scalar_name radiation_source_model = params%radiation_source_model radiation_source_scale = params%radiation_source_scale radiation_n_scalars = params%radiation_n_scalars radiation_pressure_source = params%radiation_pressure_source radiation_debug = params%radiation_debug write_radiation_diagnostics = params%write_radiation_diagnostics radiation_absorption_coeff = params%radiation_absorption_coeff radiation_emissivity = params%radiation_emissivity radiation_dom_quadrature = params%radiation_dom_quadrature call open_namelist_file(filename, unit_id, ios) if (ios == 0) then read(unit_id, nml=radiation_input, iostat=ios) close(unit_id) end if if (ios /= 0 .and. ios /= -1) then call fatal_error('input', 'failed reading &radiation_input. Check for unknown variables or typos.') end if if (ios == 0) then params%enable_radiation = enable_radiation params%radiation_update_interval = radiation_update_interval params%radiation_n_wavenumbers = radiation_n_wavenumbers params%radiation_n_species = radiation_n_species do i = 1, max_species params%radiation_species_name(i) = trim(radiation_species_name(i)) end do params%radiation_source_model = trim(lowercase(radiation_source_model)) params%radiation_source_scale = radiation_source_scale params%radiation_n_scalars = radiation_n_scalars do i = 1, max_species params%radiation_scalar_name(i) = trim(radiation_scalar_name(i)) end do params%radiation_pressure_source = trim(lowercase(radiation_pressure_source)) params%radiation_debug = radiation_debug params%write_radiation_diagnostics = write_radiation_diagnostics params%radiation_absorption_coeff = radiation_absorption_coeff params%radiation_emissivity = radiation_emissivity params%radiation_dom_quadrature = trim(lowercase(radiation_dom_quadrature)) end if end subroutine read_radiation_input !> Reads the `&output_input` namelist block. !> Reads the `&energy_input` namelist block. !! !! Reads controls for sensible-enthalpy transport. Cantera thermo, when !! enabled, preserves transported `h` and recovers `T(h,Y,p0)` every energy !! step; `thermo_update_interval` is reserved and must remain 1. subroutine read_energy_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params logical :: enable_energy logical :: enable_cantera_thermo logical :: enable_species_enthalpy_diffusion integer :: thermo_update_interval character(len=name_len) :: thermo_default_species real(rk) :: initial_T real(rk) :: energy_reference_T real(rk) :: energy_reference_h real(rk) :: energy_cp real(rk) :: energy_lambda integer :: unit_id, ios namelist /energy_input/ enable_energy, enable_cantera_thermo, enable_species_enthalpy_diffusion, & thermo_default_species, thermo_update_interval, & initial_T, energy_reference_T, energy_reference_h, & energy_cp, energy_lambda enable_energy = params%enable_energy enable_cantera_thermo = params%enable_cantera_thermo enable_species_enthalpy_diffusion = params%enable_species_enthalpy_diffusion thermo_update_interval = params%thermo_update_interval thermo_default_species = params%thermo_default_species initial_T = params%initial_T energy_reference_T = params%energy_reference_T energy_reference_h = params%energy_reference_h energy_cp = params%energy_cp energy_lambda = params%energy_lambda call open_namelist_file(filename, unit_id, ios) if (ios == 0) then read(unit_id, nml=energy_input, iostat=ios) close(unit_id) end if if (ios /= 0 .and. ios /= -1) then call fatal_error('input', 'failed reading &energy_input. Check for unknown variables or typos.') end if if (ios == 0) then params%enable_energy = enable_energy params%enable_cantera_thermo = enable_cantera_thermo params%enable_species_enthalpy_diffusion = enable_species_enthalpy_diffusion params%thermo_update_interval = thermo_update_interval params%thermo_default_species = trim(thermo_default_species) params%initial_T = initial_T params%energy_reference_T = energy_reference_T params%energy_reference_h = energy_reference_h params%energy_cp = energy_cp params%energy_lambda = energy_lambda end if end subroutine read_energy_input !> Reads the optional `&ignition_input` namelist block. !! !! The ignition kernel is an initialization aid, not a time-dependent source. !! It overwrites temperature and optionally composition inside a spherical !! region before the first timestep, then the caller refreshes thermodynamics. subroutine read_ignition_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params logical :: enable_ignition_kernel character(len=name_len) :: ignition_kernel_shape real(rk) :: ignition_kernel_center(3) real(rk) :: ignition_kernel_radius real(rk) :: ignition_kernel_T character(len=1024) :: ignition_kernel_composition character(len=name_len) :: ignition_kernel_blend integer :: unit_id, ios namelist /ignition_input/ enable_ignition_kernel, ignition_kernel_shape, & ignition_kernel_center, ignition_kernel_radius, & ignition_kernel_T, ignition_kernel_composition, & ignition_kernel_blend enable_ignition_kernel = params%enable_ignition_kernel ignition_kernel_shape = params%ignition_kernel_shape ignition_kernel_center = params%ignition_kernel_center ignition_kernel_radius = params%ignition_kernel_radius ignition_kernel_T = params%ignition_kernel_T ignition_kernel_composition = params%ignition_kernel_composition ignition_kernel_blend = params%ignition_kernel_blend call open_namelist_file(filename, unit_id, ios) if (ios == 0) then read(unit_id, nml=ignition_input, iostat=ios) close(unit_id) end if if (ios /= 0 .and. ios /= -1) then call fatal_error('input', 'failed reading &ignition_input. Check for unknown variables or typos.') end if if (ios == 0) then params%enable_ignition_kernel = enable_ignition_kernel params%ignition_kernel_shape = trim(lowercase(ignition_kernel_shape)) params%ignition_kernel_center = ignition_kernel_center params%ignition_kernel_radius = ignition_kernel_radius params%ignition_kernel_T = ignition_kernel_T params%ignition_kernel_composition = trim(ignition_kernel_composition) params%ignition_kernel_blend = trim(lowercase(ignition_kernel_blend)) end if end subroutine read_ignition_input !> Reads the optional `&restart_input` block. subroutine read_restart_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params logical :: restart_from_file character(len=path_len) :: restart_file logical :: write_restart integer :: restart_interval character(len=path_len) :: restart_output_dir integer :: restart_keep logical :: restart_use_file_dt, restart_nsteps_are_additional integer :: unit_id, ios namelist /restart_input/ restart_from_file, restart_file, write_restart, & restart_interval, restart_output_dir, restart_keep, & restart_use_file_dt, restart_nsteps_are_additional restart_from_file = params%restart_from_file restart_file = params%restart_file write_restart = params%write_restart restart_interval = params%restart_interval restart_output_dir = params%restart_output_dir restart_keep = params%restart_keep restart_use_file_dt = params%restart_use_file_dt restart_nsteps_are_additional = params%restart_nsteps_are_additional call open_namelist_file(filename, unit_id, ios) if (ios == 0) then read(unit_id, nml=restart_input, iostat=ios) close(unit_id) end if if (ios /= 0 .and. ios /= -1) then call fatal_error('input', 'failed reading &restart_input. Check for unknown variables or typos.') end if if (ios == 0) then params%restart_from_file = restart_from_file params%restart_file = trim(restart_file) params%write_restart = write_restart params%restart_interval = restart_interval params%restart_output_dir = trim(restart_output_dir) params%restart_keep = restart_keep params%restart_use_file_dt = restart_use_file_dt params%restart_nsteps_are_additional = restart_nsteps_are_additional end if end subroutine read_restart_input subroutine read_output_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params character(len=path_len) :: output_dir character(len=name_len) :: terminal_detail logical :: write_vtu, write_diagnostics integer :: terminal_interval character(len=name_len) :: vtu_format integer :: unit_id, ios namelist /output_input/ output_dir, write_vtu, write_diagnostics, terminal_interval, terminal_detail, vtu_format output_dir = params%output_dir write_vtu = params%write_vtu write_diagnostics = params%write_diagnostics terminal_interval = params%terminal_interval terminal_detail = params%terminal_detail vtu_format = params%vtu_format call open_namelist_file(filename, unit_id, ios) if (ios == 0) then read(unit_id, nml=output_input, iostat=ios) close(unit_id) end if if (ios > 0) call fatal_error('input', 'failed reading &output_input') params%output_dir = output_dir params%write_vtu = write_vtu params%write_diagnostics = write_diagnostics params%terminal_interval = terminal_interval params%terminal_detail = trim(lowercase(terminal_detail)) params%vtu_format = trim(lowercase(vtu_format)) end subroutine read_output_input !> Reads the `&species_input` namelist block. !! !! Implements strict error checking to catch typos in chemical !! species names or advection settings. subroutine read_species_input(filename, params) character(len=*), intent(in) :: filename type(case_params_t), intent(inout) :: params logical :: enable_species, enable_reactions, enable_cantera integer :: nspecies character(len=name_len), save :: species_name(max_species) = "" real(rk), save :: species_diffusivity(max_species) = zero real(rk), save :: initial_Y(max_species) = zero character(len=1024) :: initial_composition integer :: unit_id, ios namelist /species_input/ enable_species, enable_reactions, enable_cantera, & nspecies, species_name, species_diffusivity, initial_Y, & initial_composition enable_species = params%enable_species enable_reactions = params%enable_reactions enable_cantera = params%enable_cantera_species nspecies = params%nspecies species_name = params%species_name species_diffusivity = params%species_diffusivity initial_Y = params%initial_Y initial_composition = params%initial_composition call open_namelist_file(filename, unit_id, ios) if (ios == 0) then read(unit_id, nml=species_input, iostat=ios) close(unit_id) end if if (ios /= 0 .and. ios /= -1) then call fatal_error('input', 'failed reading &species_input. Check for unknown variables or typos.') end if if (ios == 0) then params%enable_species = enable_species params%enable_reactions = enable_reactions params%enable_cantera_species = enable_cantera params%nspecies = nspecies params%species_name = species_name params%species_diffusivity = species_diffusivity params%initial_Y = initial_Y params%initial_composition = trim(initial_composition) ! Save original namelist values params%namelist_nspecies = nspecies params%namelist_species_name = species_name end if end subroutine read_species_input !> Helper routine to safely open a namelist file for reading. subroutine open_namelist_file(filename, unit_id, ios) character(len=*), intent(in) :: filename integer, intent(out) :: unit_id integer, intent(out) :: ios open(newunit=unit_id, file=trim(filename), status='old', action='read', iostat=ios) if (ios /= 0) unit_id = -1 end subroutine open_namelist_file !> Validates all parsed parameters against physical and algorithmic limits. !! !! Triggers `fatal_error` for non-physical values (e.g., negative density, !! zero timestep) or unsupported configurations. subroutine validate_params(params) type(case_params_t), intent(in) :: params if (params%nsteps < 0) call fatal_error('input', 'nsteps must be non-negative') if (params%dt <= zero) call fatal_error('input', 'dt must be positive') if (params%max_dt <= zero) call fatal_error('input', 'max_dt must be positive') if (params%min_dt <= zero) call fatal_error('input', 'min_dt must be positive') if (params%min_dt > params%max_dt) call fatal_error('input', 'min_dt must be <= max_dt') if (params%dt_growth_limit < one) call fatal_error('input', 'dt_growth_limit must be >= 1') if (params%max_step_retries < 0) call fatal_error('input', 'max_step_retries must be non-negative') if (params%step_reject_cut_factor <= zero .or. params%step_reject_cut_factor >= one) then call fatal_error('input', 'step_reject_cut_factor must be in (0,1)') end if if (params%step_reject_min_dt <= zero) call fatal_error('input', 'step_reject_min_dt must be positive') if (params%max_KE_growth_per_step < zero) call fatal_error('input', 'max_KE_growth_per_step must be non-negative') if (params%ke_reject_floor < zero) call fatal_error('input', 'ke_reject_floor must be non-negative') if (params%max_cfl_overshoot_factor < zero) call fatal_error('input', 'max_cfl_overshoot_factor must be non-negative') if (params%output_interval <= 0) call fatal_error('input', 'output_interval must be positive') if (params%terminal_interval <= 0) call fatal_error('input', 'terminal_interval must be positive') select case (trim(params%terminal_detail)) case ('brief', 'health', 'full') continue case default call fatal_error('input', 'terminal_detail must be one of: brief, health, full') end select if (params%rho <= zero) call fatal_error('input', 'rho must be positive') if (params%nu < zero) call fatal_error('input', 'nu must be non-negative') if (params%enable_variable_nu .and. .not. params%enable_cantera_fluid) then call fatal_error('input', 'enable_variable_nu=.true. requires fluid_input enable_cantera=.true.; otherwise use constant nu') end if if (params%transport_update_interval <= 0) call fatal_error('input', 'transport_update_interval must be positive') if (params%pressure_max_iter <= 0) call fatal_error('input', 'pressure_max_iter must be positive') call validate_advection_scheme_local(params%momentum_convection_scheme, 'momentum_convection_scheme', .false.) call validate_advection_scheme_local(params%convection_scheme, 'convection_scheme', .false.) call validate_advection_scheme_local(params%species_convection_scheme, 'species_convection_scheme', .true.) call validate_advection_scheme_local(params%energy_convection_scheme, 'energy_convection_scheme', .true.) call validate_limiter_name_local(params%scalar_limiter) call validate_time_scheme_local(params%species_time_scheme, 'species_time_scheme') call validate_time_scheme_local(params%energy_time_scheme, 'energy_time_scheme') if (params%initial_T <= zero) call fatal_error('input', 'initial_T must be positive') if (params%energy_reference_T <= zero) call fatal_error('input', 'energy_reference_T must be positive') if (params%energy_cp <= zero) call fatal_error('input', 'energy_cp must be positive') if (params%energy_lambda < zero) call fatal_error('input', 'energy_lambda must be non-negative') if (params%thermo_update_interval <= 0) call fatal_error('input', 'thermo_update_interval must be positive') if (params%enable_cantera_thermo .and. params%thermo_update_interval /= 1) then call fatal_error('input', 'thermo_update_interval values other than 1 are reserved; Cantera thermo is currently updated every energy step') end if if (params%enable_cantera_thermo .and. len_trim(params%thermo_default_species) == 0) & call fatal_error('input', 'thermo_default_species must not be empty when enable_cantera_thermo is true') if (params%enable_species_enthalpy_diffusion) then if (.not. params%enable_energy) then call fatal_error('input', 'enable_species_enthalpy_diffusion requires enable_energy=.true.') end if if (.not. params%enable_species) then call fatal_error('input', 'enable_species_enthalpy_diffusion requires transported species') end if if (params%nspecies <= 0 .and. .not. params%enable_reactions .and. & .not. (params%enable_species .and. params%enable_cantera_species)) then call fatal_error('input', 'enable_species_enthalpy_diffusion requires nspecies > 0 unless Cantera auto-discovers mechanism species') end if if (.not. params%enable_cantera_thermo) then call fatal_error('input', 'enable_species_enthalpy_diffusion currently requires enable_cantera_thermo=.true. for species sensible enthalpies') end if end if if (params%n_patches < 0 .or. params%n_patches > max_patches) then call fatal_error('input', 'n_patches is outside supported range') end if if (len_trim(params%mesh_dir) == 0) then call fatal_error('input', 'mesh_dir cannot be empty') end if if (len_trim(params%output_dir) == 0) then call fatal_error('input', 'output_dir cannot be empty') end if if (params%vtu_format /= "ascii" .and. params%vtu_format /= "binary") then call fatal_error('input', 'vtu_format must be either "ascii" or "binary"') end if if (params%nspecies < 0 .or. params%nspecies > max_species) then call fatal_error('input', 'nspecies is outside supported range') end if if (params%enable_reactions) then if (.not. params%enable_species) call fatal_error('input', 'enable_reactions requires enable_species=.true.') if (.not. params%enable_energy) call fatal_error('input', 'enable_reactions requires enable_energy=.true.') if (.not. params%enable_cantera_thermo) call fatal_error('input', 'enable_reactions requires enable_cantera_thermo=.true.') if (params%chemistry_update_interval <= 0) call fatal_error('input', 'chemistry_update_interval must be positive') if (params%chemistry_rtol <= zero) call fatal_error('input', 'chemistry_rtol must be positive') if (params%chemistry_atol <= zero) call fatal_error('input', 'chemistry_atol must be positive') if (params%chemistry_max_steps <= 0) call fatal_error('input', 'chemistry_max_steps must be positive') if (params%chemistry_temperature_cutoff < zero) call fatal_error('input', 'chemistry_temperature_cutoff must be non-negative') if (params%chemistry_min_reactive_mass_fraction < zero) call fatal_error('input', 'chemistry_min_reactive_mass_fraction must be non-negative') if (params%chemistry_active_species_threshold < zero) call fatal_error('input', 'chemistry_active_species_threshold must be non-negative') if (params%chemistry_max_dT_per_step < zero) call fatal_error('input', 'chemistry_max_dT_per_step must be non-negative') if (params%chemistry_max_rel_rho_change_per_step < zero) call fatal_error('input', 'chemistry_max_rel_rho_change_per_step must be non-negative') if (params%chemistry_max_dY_per_step < zero) call fatal_error('input', 'chemistry_max_dY_per_step must be non-negative') if (params%chemistry_dt_safety <= zero .or. params%chemistry_dt_safety > one) call fatal_error('input', 'chemistry_dt_safety must be in (0,1]') if (params%chemistry_dt_min_factor <= zero .or. params%chemistry_dt_min_factor > one) call fatal_error('input', 'chemistry_dt_min_factor must be in (0,1]') if (params%chemistry_source_relaxation <= zero .or. params%chemistry_source_relaxation > one) call fatal_error('input', 'chemistry_source_relaxation must be in (0,1]') if (params%chemistry_n_active_species < 0 .or. params%chemistry_n_active_species > max_species) then call fatal_error('input', 'chemistry_n_active_species must be between 0 and max_species') end if end if if (params%variable_density_debug_interval <= 0) then call fatal_error('input', 'variable_density_debug_interval must be positive') end if if (params%enable_ignition_kernel) then if (.not. params%enable_energy) call fatal_error('input', 'enable_ignition_kernel requires enable_energy=.true.') select case (trim(lowercase(params%ignition_kernel_shape))) case ('sphere') continue case default call fatal_error('input', 'ignition_kernel_shape currently supports only sphere') end select select case (trim(lowercase(params%ignition_kernel_blend))) case ('overwrite') continue case default call fatal_error('input', 'ignition_kernel_blend currently supports only overwrite') end select if (params%ignition_kernel_radius <= zero) call fatal_error('input', 'ignition_kernel_radius must be positive') if (params%ignition_kernel_T <= zero) call fatal_error('input', 'ignition_kernel_T must be positive') if (len_trim(params%ignition_kernel_composition) > 0 .and. .not. params%enable_species) then call fatal_error('input', 'ignition_kernel_composition requires enable_species=.true.') end if end if if (params%restart_from_file .and. len_trim(params%restart_file) == 0) then call fatal_error('input', 'restart_from_file=.true. requires restart_file') end if if (params%write_restart) then if (params%restart_interval <= 0) call fatal_error('input', 'write_restart=.true. requires restart_interval > 0') if (len_trim(params%restart_output_dir) == 0) call fatal_error('input', 'restart_output_dir cannot be empty') if (params%restart_keep < 0) call fatal_error('input', 'restart_keep must be non-negative') end if if (params%enable_radiation) then if (.not. params%enable_energy) call fatal_error('input', 'enable_radiation requires enable_energy=.true.') if (params%radiation_update_interval <= 0) call fatal_error('input', 'radiation_update_interval must be positive') if (params%radiation_n_wavenumbers <= 0) call fatal_error('input', 'radiation_n_wavenumbers must be positive when radiation is enabled') if (params%radiation_n_species < 0 .or. params%radiation_n_species > max_species) then call fatal_error('input', 'radiation_n_species is outside supported range') end if if (params%radiation_n_species > 0 .and. .not. params%enable_species) then call fatal_error('input', 'radiation species selection requires enable_species=.true.') end if select case (trim(lowercase(params%radiation_source_model))) case ('none', 'spectral_test', 'external', 'p1', 'dom') continue case default call fatal_error('input', 'radiation_source_model must be one of: none, spectral_test, external, p1, dom') end select if (trim(lowercase(params%radiation_source_model)) == 'dom') then if (trim(lowercase(params%radiation_dom_quadrature)) /= 's2' .and. & trim(lowercase(params%radiation_dom_quadrature)) /= 's4' .and. & trim(lowercase(params%radiation_dom_quadrature)) /= 's6' .and. & trim(lowercase(params%radiation_dom_quadrature)) /= 's8') then call fatal_error('input', 'radiation_dom_quadrature must be s2, s4, s6, or s8') end if end if if (params%radiation_n_scalars < 0 .or. params%radiation_n_scalars > max_species) then call fatal_error('input', 'radiation_n_scalars is outside supported range') end if if (trim(lowercase(params%radiation_pressure_source)) /= 'background' .and. & trim(lowercase(params%radiation_pressure_source)) /= 'system') then call fatal_error('input', 'radiation_pressure_source must be one of: background, system') end if end if call validate_boundary_arrays(params) end subroutine validate_params subroutine validate_advection_scheme_local(scheme_name, field_name, allow_scalar_high_order) character(len=*), intent(in) :: scheme_name character(len=*), intent(in) :: field_name logical, intent(in) :: allow_scalar_high_order character(len=len(scheme_name)) :: scheme scheme = trim(lowercase(scheme_name)) if (len_trim(scheme) == 0) scheme = 'upwind' select case (trim(scheme)) case ('upwind', 'first_order_upwind', 'first-order-upwind', & 'central', 'central_difference', 'central-difference') return case ('bounded_central', 'bounded-central', 'bounded_linear', 'bounded-linear', & 'deferred_central', 'deferred-correction', & 'muscl', 'limited_linear', 'limited-linear', 'linear_upwind', 'linear-upwind', & 'higher_order_bounded', 'higher-order-bounded', & 'quick', 'hlpa', 'quick_limited', 'quick-limited') if (allow_scalar_high_order) return case default continue end select if (allow_scalar_high_order) then call fatal_error('input', trim(field_name)//' must be one of: upwind, central, bounded_central, bounded_linear, muscl, limited_linear, quick, hlpa') else call fatal_error('input', trim(field_name)//' must be one of: upwind, central') end if end subroutine validate_advection_scheme_local subroutine validate_limiter_name_local(limiter_name) character(len=*), intent(in) :: limiter_name character(len=len(limiter_name)) :: limiter limiter = trim(lowercase(limiter_name)) if (len_trim(limiter) == 0) limiter = 'barth_jespersen' select case (trim(limiter)) case ('none', 'minmod', 'vanleer', 'van_leer', 'mc', 'monotonized_central', & 'superbee', 'barth_jespersen', 'barth-jespersen', 'venkatakrishnan', 'hlpa') return case default call fatal_error('input', 'scalar_limiter must be one of: none, minmod, vanleer, mc, superbee, barth_jespersen, venkatakrishnan, hlpa') end select end subroutine validate_limiter_name_local subroutine validate_time_scheme_local(scheme_name, field_name) character(len=*), intent(in) :: scheme_name character(len=*), intent(in) :: field_name character(len=len(scheme_name)) :: scheme scheme = trim(lowercase(scheme_name)) if (len_trim(scheme) == 0) scheme = 'euler' select case (trim(scheme)) case ('euler', 'forward_euler', 'forward-euler', 'ab2', 'adams_bashforth2', 'adams-bashforth2') return case default call fatal_error('input', trim(field_name)//' must be one of: euler, ab2') end select end subroutine validate_time_scheme_local !> Ensures all boundary patch names and types are non-empty. subroutine validate_boundary_arrays(params) type(case_params_t), intent(in) :: params integer :: i do i = 1, params%n_patches if (len_trim(params%patch_name(i)) == 0) then call fatal_error('input', 'patch_name entry cannot be empty') end if if (len_trim(params%patch_type(i)) == 0) then call fatal_error('input', 'patch_type entry cannot be empty') end if if (params%enable_energy .and. params%patch_T(i) <= zero) then call fatal_error('input', 'patch_T entry must be positive when energy is enabled') end if end do end subroutine validate_boundary_arrays end module mod_input