!> Main entry point for the LowMachReact-Hex finite-volume solver. !! !! This program orchestrates the entire simulation lifecycle: !! 1. **Initialization**: Starts MPI, parses the `case.nml` namelist, and reads the mesh. !! 2. **Domain Decomposition**: Sets up the replicated mesh MPI ranks for flow and radiation. !! 3. **Field Setup**: Allocates and initializes velocity, pressure, species, and energy fields. !! 4. **Time Integration Loop**: Executes projection, non-reacting species transport, !! and sensible-enthalpy transport. Constant-density mode is the stable !! baseline; guarded variable-density mode syncs active density from !! Cantera thermo density and projects to \( \nabla \cdot u=S \). !! 5. **Diagnostics & Output**: Computes global observables and writes VTU/CSV files. !! 6. **Finalization**: Safely releases all allocated memory and shuts down MPI. program lowmach_react_hex use mpi_f08 use mod_precision, only : rk, zero, one, tiny_safe, output_unit, fatal_error use mod_input, only : case_params_t, read_case_params, max_species use mod_mesh, only : mesh_t, mesh_finalize use mod_mesh_io, only : read_native_mesh use mod_mpi_flow, only : flow_mpi_t, mpi_flow_startup, mpi_flow_shutdown, & flow_mpi_initialize, flow_mpi_finalize, & flow_allgather_owned_scalar, flow_allgather_owned_matrix use mod_mpi_radiation, only : radiation_mpi_t, radiation_mpi_initialize, & radiation_mpi_finalize, radiation_configure_spectral use mod_bc, only : bc_set_t, build_bc_set, finalize_bc_set use mod_species_composition, only : parse_named_composition, composition_string_is_set use mod_flow_fields, only : flow_fields_t, initialize_fields, finalize_fields use mod_flow_projection, only : solver_stats_t, advance_projection_step, compute_flow_diagnostics, compute_and_update_cfl, & update_lowmach_divergence_source_from_density, correct_projection_to_lowmach_source use mod_output, only : & prepare_output, & write_diagnostics_header, & write_diagnostics_row, & write_vtu_unstructured, & write_mesh_summary, & write_pvd_collection, & write_variable_density_diagnostics, & write_boundary_mass_balance_diagnostics, & write_species_energy_conservation_diagnostics, & load_existing_pvd use mod_species, only : species_fields_t, initialize_species, finalize_species, & advance_species_transport use mod_energy, only : & energy_fields_t, & initialize_energy, & finalize_energy, & refresh_energy_from_temperature, & advance_energy_transport, & advance_cantera_chemistry, & write_cantera_cache_stats, write_chemistry_screening_stats, & write_energy_diagnostics_header, & write_energy_diagnostics_row, & write_chemistry_diagnostics_header use mod_cantera, only : transport_properties_t, initialize_transport, & finalize_transport, update_transport_properties, sync_active_density_from_thermo use mod_restart, only : read_restart_file, write_restart_file use mod_radiation, only : radiation_context_t, radiation_initialize, radiation_update_source, radiation_finalize use mod_profiling, only : profiler_start, profiler_stop, profiler_report, profiler_configure, & profiler_register_standard_timers implicit none type :: step_backup_t type(flow_fields_t) :: fields type(species_fields_t) :: species type(energy_fields_t) :: energy type(transport_properties_t) :: transport type(solver_stats_t) :: stats real(rk) :: time = zero real(rk) :: dt = zero real(rk) :: chemistry_accum_dt = zero logical :: valid = .false. end type step_backup_t type(case_params_t) :: params !< Parsed simulation parameters. type(mesh_t) :: mesh !< Global computational mesh. type(flow_mpi_t) :: flow_mpi !< MPI context for flow solver. type(radiation_mpi_t) :: rad_mpi !< MPI context for radiation solver. type(radiation_context_t) :: rad_context !< Radiation model/data-exchange context. type(bc_set_t) :: bc !< Boundary condition definitions. type(flow_fields_t) :: fields !< Hydrodynamic fields (U, P). type(solver_stats_t) :: stats !< Runtime diagnostics and solver performance metrics. type(species_fields_t) :: species !< Species mass fractions. type(energy_fields_t) :: energy !< Enthalpy/temperature/radiation-source fields. type(transport_properties_t) :: transport !< Physical properties (rho, mu, Dk). character(len=256) :: case_file integer :: step, start_step, ierr, world_rank logical :: do_output_step, do_terminal_step, do_status_step, do_cfl_step logical :: do_terminal_full, do_terminal_health logical :: is_chem_step, is_rad_step real(rk) :: time, t0, chemistry_accum_dt real(rk) :: console_max_divu, console_sproj_max real(rk) :: console_lm_res_projection, console_lm_res_current real(rk) :: console_mass_balance_defect real(rk) :: terminal_projection_metric, terminal_ke, terminal_tmin, terminal_tmax real(rk) :: chem_max_dT, chem_max_dY, chem_max_rel_drho, chem_min_alpha real(rk) :: chem_raw_max_dT, chem_raw_max_dY, chem_raw_max_rel_drho type(step_backup_t) :: step_backup integer :: step_retry_count logical :: step_accepted, reject_step, block_dynamic_dt character(len=256) :: reject_reason real(rk) :: step_ke_before, step_ke_after, step_cfl_after, step_umax_after real(rk) :: step_trial_dt, step_cut_dt ! 1. Initialize MPI and read runtime configuration. call mpi_flow_startup() call MPI_Comm_rank(MPI_COMM_WORLD, world_rank, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI_Comm_rank failed') t0 = real(mpi_wtime(), rk) step = 0 start_step = 0 time = zero chemistry_accum_dt = zero call get_case_filename(case_file) call read_case_params(trim(case_file), params) call profiler_configure(params%enable_profiling, params%nested_profiling) call profiler_register_standard_timers() if (world_rank == 0) then write(output_unit,'(a,l1,a,l1,a,l1)') 'profiling config: enable=', & params%enable_profiling, ' nested=', params%nested_profiling, & ' cantera_cache_stats=', params%write_cantera_cache_stats end if call profiler_start('Total_Simulation') ! 2. Setup simulation environment and read input data. call read_native_mesh(trim(params%mesh_dir), mesh) ! 3. Initialize parallel contexts. call initialize_transport(mesh, params, transport) call flow_mpi_initialize(mesh, flow_mpi, MPI_COMM_WORLD, max(4, params%nspecies)) call radiation_mpi_initialize(rad_mpi, MPI_COMM_WORLD) if (params%enable_radiation) then call radiation_configure_spectral(rad_mpi, params%radiation_n_wavenumbers) end if ! 4. Prepare physical fields and BCs. call build_bc_set(mesh, params, bc) call initialize_fields(mesh, fields) if (params%enable_species) then call initialize_species(mesh, params, species) end if if (params%enable_energy) then if (params%enable_species) then call initialize_energy(mesh, params, energy, species%Y) else call initialize_energy(mesh, params, energy) end if if (params%enable_ignition_kernel .and. .not. params%restart_from_file) then call apply_ignition_kernel(mesh, flow_mpi, params, species, energy) if (params%enable_species) then call refresh_energy_from_temperature(mesh, params, energy, species%Y) else call refresh_energy_from_temperature(mesh, params, energy) end if end if ! No-op for the default constant-density solver. In guarded ! density_eos="cantera" variable-density mode, energy%rho_thermo becomes ! the active flow density after thermo sync. call sync_active_density_from_thermo(mesh, flow_mpi, params, transport, energy%rho_thermo) ! Initial variable-density state: rho_old must equal active rho at t=0. ! The sync routine stores the pre-sync placeholder density in rho_old; ! that is useful during timestepping but not for the initial state. if (params%enable_variable_density) then if (allocated(transport%rho_old)) transport%rho_old = transport%rho end if end if ! Optional rank-independent restart. The file stores global arrays in mesh ! order, so this works even when the new run uses a different MPI rank count. call read_restart_file(mesh, flow_mpi, params, fields, species, energy, transport, start_step, time, chemistry_accum_dt) ! If ignition kernel is requested on a restarted run, overwrite the fields ! immediately after loading and re-refresh enthalpy. if (params%enable_ignition_kernel .and. params%restart_from_file) then call apply_ignition_kernel(mesh, flow_mpi, params, species, energy) if (params%enable_species) then call refresh_energy_from_temperature(mesh, params, energy, species%Y) else call refresh_energy_from_temperature(mesh, params, energy) end if end if ! Align active flow density with thermodynamic density immediately after restart ! to prevent catastrophic startup history-source transients. if (params%enable_variable_density) then call sync_active_density_from_thermo(mesh, flow_mpi, params, transport, energy%rho_thermo) if (allocated(transport%rho_old)) transport%rho_old = transport%rho end if params%dt_old = params%dt if (flow_mpi%rank == 0 .and. params%enable_species) then print *, "species: ", species%nspecies end if ! 5. Setup output files. call prepare_output(params, flow_mpi) call execute_command_line('mkdir -p ' // trim(params%output_dir) // '/diagnostics') call execute_command_line('mkdir -p ' // trim(params%output_dir) // '/VTK') call write_mesh_summary(params, flow_mpi, mesh) call write_diagnostics_header(params, flow_mpi) call write_energy_diagnostics_header(params, flow_mpi) call write_chemistry_diagnostics_header(params, flow_mpi) call load_existing_pvd(params, flow_mpi, start_step) call radiation_initialize(mesh, rad_mpi, params, rad_context) call write_run_config_summary(mesh, flow_mpi, rad_mpi, params, trim(case_file)) ! 6. Initial diagnostics and visualization snapshot. step = start_step call compute_flow_diagnostics(mesh, flow_mpi, bc, params, fields, stats) call write_diagnostics_row(params, flow_mpi, step, time, stats) call write_boundary_mass_balance_diagnostics(mesh, flow_mpi, params, fields, transport, step, time) if (params%enable_variable_density .and. step > 0 .and. time > 0.0_rk) then call write_variable_density_diagnostics(mesh, flow_mpi, params, fields, energy, transport, step, time) end if call write_energy_diagnostics_row(mesh, flow_mpi, params, energy, step, time, write_file=.true.) call write_vtu_unstructured(params, flow_mpi, mesh, fields, species, energy, transport, step, time) if (flow_mpi%rank == 0) then write(output_unit,'(a)') 'LowMachReact-Hex hexahedral low-Mach FV solver' write(output_unit,'(a,a)') 'case: ', trim(case_file) write(output_unit,'(a,i0)') 'cells: ', mesh%ncells write(output_unit,'(a,i0)') 'flow/chemistry MPI ranks: ', flow_mpi%nprocs write(output_unit,'(a,i0)') 'radiation MPI ranks: ', rad_mpi%nprocs if (params%enable_radiation) then write(output_unit,'(a,i0,a,i0,a,i0)') 'Radiation spectral tasks: n_wavenumbers=', & params%radiation_n_wavenumbers, ' local(first)=', rad_mpi%first_wavenumber, & ' update_interval=', params%radiation_update_interval write(output_unit,'(a,a)') 'Radiation source model: ', trim(params%radiation_source_model) write(output_unit,'(a)') 'Radiation coupling: spectral communicator computes qrad on full replicated mesh and reuses qrad between updates' end if if (params%enable_variable_density) then write(output_unit,'(a,a)') 'Flow density mode: EXPERIMENTAL variable-density low-Mach, density_eos = ', & trim(params%density_eos) write(output_unit,'(a)') 'Flow density source: transport%rho <- energy%rho_thermo after thermo sync' write(output_unit,'(a)') 'Low-Mach constraint: div(u) = S with conservative density-history plus advective-density source' if (params%variable_density_debug) then write(output_unit,'(a)') 'Variable-density debug diagnostics: enabled' else write(output_unit,'(a)') 'Variable-density debug diagnostics: disabled' end if else write(output_unit,'(a,es12.5,a)') 'Flow density mode: constant rho = ', params%rho, ' kg/m^3' end if if (params%enable_variable_nu) then write(output_unit,'(a)') 'Flow viscosity mode: variable Cantera mu/nu enabled' else write(output_unit,'(a,es12.5,a)') 'Flow viscosity mode: constant nu = ', params%nu, ' m^2/s' end if if (params%enable_cantera_fluid .or. params%enable_cantera_species .or. params%enable_cantera_thermo) then if (len_trim(params%cantera_phase_name) > 0) then write(output_unit,'(a,a)') 'Cantera phase selector: ', trim(params%cantera_phase_name) else write(output_unit,'(a)') 'Cantera phase selector: ' end if end if if (params%enable_energy .and. params%enable_cantera_thermo) then if (params%enable_variable_density .and. trim(params%density_eos) == 'cantera') then write(output_unit,'(a)') 'Cantera rho_thermo: ACTIVE EOS density for experimental low-Mach coupling' else write(output_unit,'(a)') 'Cantera rho_thermo: diagnostic only, not used by projection' end if end if if (params%enable_cantera_fluid .or. params%enable_cantera_species) then write(output_unit,'(a,i0,a)') 'Cantera transport update interval: ', & params%transport_update_interval, ' step(s) [mu/D_k only]' if (params%enable_energy) then write(output_unit,'(a)') 'Cantera transport temperature source: energy%T' else write(output_unit,'(a,es12.5,a)') 'Cantera transport temperature source: background_temp = ', & params%background_temp, ' K' end if end if if (params%enable_energy .and. params%enable_species_enthalpy_diffusion) then if (params%enable_variable_density) then write(output_unit,'(a)') 'Species enthalpy diffusion: enabled with density-weighted variable-density branch' else write(output_unit,'(a)') 'Species enthalpy diffusion: enabled [-div(sum_k h_k J_k)/rho]' end if end if if (params%enable_energy .and. params%enable_cantera_thermo) then write(output_unit,'(a)') 'Cantera thermo update interval: every energy step' write(output_unit,'(a,es12.5,a)') 'Cantera thermodynamic pressure p0: ', & params%background_press, ' Pa' end if if (params%enable_reactions) then write(output_unit,'(a,i0,a,es10.3,a,es10.3,a,i0)') 'Cantera chemistry: enabled, update_interval=', & params%chemistry_update_interval, ' rtol=', params%chemistry_rtol, & ' atol=', params%chemistry_atol, ' max_steps=', params%chemistry_max_steps end if end if ! 7. Main Time-Stepping Loop do step = start_step + 1, params%nsteps step_retry_count = 0 step_accepted = .false. do while (.not. step_accepted) if (params%enable_step_rejection) then call save_step_backup(fields, species, energy, transport, stats, time, params%dt, chemistry_accum_dt, step_backup) end if call compute_step_health(mesh, flow_mpi, params, fields, transport, step_ke_before, step_cfl_after, step_umax_after) chem_max_dT = zero chem_max_dY = zero chem_max_rel_drho = zero chem_min_alpha = 1.0_rk chem_raw_max_dT = zero chem_raw_max_dY = zero chem_raw_max_rel_drho = zero do_output_step = (mod(step, params%output_interval) == 0) .or. step == params%nsteps do_terminal_step = (mod(step, params%terminal_interval) == 0) .or. step == params%nsteps do_terminal_full = trim(params%terminal_detail) == 'full' do_terminal_health = trim(params%terminal_detail) == 'health' do_status_step = do_output_step .or. (do_terminal_step .and. do_terminal_full) do_cfl_step = params%use_dynamic_dt .or. do_output_step .or. & (do_terminal_step .and. (do_terminal_health .or. do_terminal_full)) ! A. Adaptive time-stepping / CFL diagnostics. ! Dynamic dt requires CFL every step. Health/full terminal modes also ! report CFL, while brief mode avoids this diagnostic reduction unless ! output or dynamic timestep control needs it. if (do_cfl_step) then call profiler_start('CFL_Update') if (step_retry_count > 0) then block_dynamic_dt = params%use_dynamic_dt params%use_dynamic_dt = .false. call compute_and_update_cfl(mesh, flow_mpi, params, fields, stats) params%use_dynamic_dt = block_dynamic_dt else call compute_and_update_cfl(mesh, flow_mpi, params, fields, stats) end if call profiler_stop('CFL_Update') end if ! B. Dynamic transport-property update via Cantera. ! transport_update_interval gates Cantera mu/D_k refresh; enable_variable_nu controls whether mu affects flow. ! In variable-density mode, active flow density is refreshed after the energy thermo sync. if (mod(step-1, params%transport_update_interval) == 0 .or. step == 1) then call profiler_start('Transport_Update') if (params%enable_species) then if (params%enable_energy) then call update_transport_properties(mesh, flow_mpi, params, species%Y, transport, T_state=energy%T) else call update_transport_properties(mesh, flow_mpi, params, species%Y, transport) end if else if (params%enable_energy) then call update_transport_properties(mesh, flow_mpi, params, transport=transport, T_state=energy%T) else call update_transport_properties(mesh, flow_mpi, params, transport=transport) end if end if call profiler_stop('Transport_Update') end if ! C. Advance Momentum & Pressure (Projection Method). call profiler_start('Projection_Step') if (params%enable_variable_density .and. allocated(fields%projection_rho)) then fields%projection_rho = transport%rho end if call advance_projection_step(mesh, flow_mpi, bc, params, transport, fields, stats) call profiler_stop('Projection_Step') ! D. Advance Species Transport. if (params%enable_species) then call profiler_start('Species_Transport') call advance_species_transport(mesh, flow_mpi, bc, params, fields, species, transport) call profiler_stop('Species_Transport') end if ! Unified Physics Synchronization Step for load-balancing if (params%enable_chemistry_load_balancing) then is_chem_step = params%enable_reactions .and. & (mod(step - start_step, params%chemistry_update_interval) == 0 .or. step == params%nsteps) is_rad_step = params%enable_energy .and. params%enable_radiation .and. & (mod(step - 1, params%radiation_update_interval) == 0 .or. step == 1) if (is_chem_step .or. is_rad_step) then call profiler_start('Unified_Physics_Sync') call flow_allgather_owned_scalar(flow_mpi, energy%T, energy%T) call flow_allgather_owned_matrix(flow_mpi, species%Y, species%Y) call profiler_stop('Unified_Physics_Sync') end if end if ! E. Optional cell-local Cantera chemistry after species transport. ! When chemistry_update_interval > 1, integrate over the accumulated ! skipped flow time instead of only one CFD dt. if (params%enable_reactions) then chemistry_accum_dt = chemistry_accum_dt + params%dt if (mod(step - start_step, params%chemistry_update_interval) == 0 .or. step == params%nsteps) then call profiler_start('Chemistry_Update') call advance_cantera_chemistry(mesh, flow_mpi, params, species%Y, energy, step, & chemistry_dt=chemistry_accum_dt, & max_dT_chem=chem_max_dT, max_dY_chem=chem_max_dY, & max_rel_drho_chem=chem_max_rel_drho, min_source_alpha=chem_min_alpha, & raw_max_dT_chem=chem_raw_max_dT, raw_max_dY_chem=chem_raw_max_dY, & raw_max_rel_drho_chem=chem_raw_max_rel_drho) chemistry_accum_dt = zero call profiler_stop('Chemistry_Update') end if end if ! F. Optional spectral radiation update before energy transport. On ! non-update steps this call preserves the existing qrad field so the ! last radiation source is reused. if (params%enable_energy .and. params%enable_radiation) then if (params%enable_species) then call radiation_update_source(mesh, bc, flow_mpi, rad_mpi, rad_context, params, energy, species%Y, step, time, params%dt) else call radiation_update_source(mesh, bc, flow_mpi, rad_mpi, rad_context, params, energy, step=step, time=time, dt=params%dt) end if end if ! G. Advance passive/reacting sensible-enthalpy transport. if (params%enable_energy) then call profiler_start('Energy_Transport') if (params%enable_variable_density) then call sync_active_density_from_thermo(mesh, flow_mpi, params, transport, energy%rho_thermo) end if if (params%enable_species) then call advance_energy_transport(mesh, flow_mpi, bc, params, fields, energy, species%Y, transport, step) else call advance_energy_transport(mesh, flow_mpi, bc, params, fields, energy, transport=transport, step=step) end if ! No-op for enable_variable_density=.false.; guarded variable-density ! mode refreshes transport%rho after the thermo sync. call sync_active_density_from_thermo(mesh, flow_mpi, params, transport, energy%rho_thermo) if (params%enable_species) then call update_lowmach_divergence_source_from_density(mesh, flow_mpi, bc, params, transport, fields, species%Y, energy%T) else call update_lowmach_divergence_source_from_density(mesh, flow_mpi, bc, params, transport, fields, T_state=energy%T) end if ! Post-density corrector projection (no-op unless enable_density_corrector_projection=.true.). if (params%enable_density_corrector_projection) then call profiler_start('Density_Corrector_Projection') call correct_projection_to_lowmach_source(mesh, flow_mpi, bc, params, transport, fields, stats) call profiler_stop('Density_Corrector_Projection') end if call profiler_stop('Energy_Transport') end if call compute_step_health(mesh, flow_mpi, params, fields, transport, step_ke_after, step_cfl_after, step_umax_after) stats%cfl = step_cfl_after stats%kinetic_energy = step_ke_after stats%max_velocity = step_umax_after call evaluate_step_rejection(params, stats, chem_max_dT, chem_max_dY, chem_max_rel_drho, & step_ke_before, step_ke_after, step_cfl_after, fields, energy, transport, & reject_step, reject_reason) if (reject_step) then step_trial_dt = params%dt step_retry_count = step_retry_count + 1 call restore_step_backup(step_backup, fields, species, energy, transport, stats, time, params%dt, chemistry_accum_dt) fields%rhs_old_valid = .false. if (params%enable_species) species%rhs_old_valid = .false. if (params%enable_energy) energy%rhs_old_valid = .false. step_cut_dt = max(max(params%step_reject_min_dt, params%min_dt), step_trial_dt * params%step_reject_cut_factor) params%dt = min(step_cut_dt, params%max_dt) if (flow_mpi%rank == 0) then write(output_unit,'(a,i0,a,i0,a,es12.5,a,es12.5,a,a)') & 'step rejection: step=', step, ' retry=', step_retry_count, & ' old_dt=', step_trial_dt, ' new_dt=', params%dt, ' reason=', trim(reject_reason) end if if (step_retry_count > params%max_step_retries) then call fatal_error('main', 'maximum step retries exceeded at step '//trim(int_to_string(step))//': '//trim(reject_reason)) end if if (params%stop_on_min_dt_failure .and. params%dt <= params%step_reject_min_dt * (1.0_rk + 10.0_rk*tiny_safe)) then call fatal_error('main', 'step rejected at minimum dt at step '//trim(int_to_string(step))//': '//trim(reject_reason)) end if cycle end if step_accepted = .true. end do time = time + params%dt params%dt_old = params%dt if (params%enable_reactions) then call apply_chemistry_timestep_control(flow_mpi, params, chem_max_dT, chem_max_dY, & chem_max_rel_drho, chem_min_alpha) end if call write_restart_file(mesh, flow_mpi, params, fields, species, energy, transport, step, time, chemistry_accum_dt) ! E. Periodic terminal status, diagnostics, and file output. if (do_status_step) then stats%wall_time = real(mpi_wtime(), rk) - t0 call profiler_start('Flow_Diagnostics') call compute_global_observables(mesh, flow_mpi, fields, species, transport, stats) call profiler_stop('Flow_Diagnostics') if (do_output_step) then call profiler_start('Diagnostics_Write_Flow') call write_diagnostics_row(params, flow_mpi, step, time, stats) call write_boundary_mass_balance_diagnostics(mesh, flow_mpi, params, fields, transport, step, time) if (params%enable_variable_density .and. step > 0 .and. time > 0.0_rk) then call write_variable_density_diagnostics(mesh, flow_mpi, params, fields, energy, transport, step, time) end if call profiler_stop('Diagnostics_Write_Flow') call profiler_start('Output_Write_VTU') call write_vtu_unstructured(params, flow_mpi, mesh, fields, species, energy, transport, step, time) call profiler_stop('Output_Write_VTU') end if if (params%enable_energy .and. do_output_step) then call profiler_start('Diagnostics_Write_Energy') call write_energy_diagnostics_row(mesh, flow_mpi, params, energy, step, time, & write_file=do_output_step) call profiler_stop('Diagnostics_Write_Energy') end if end if if (do_terminal_step) then stats%wall_time = real(mpi_wtime(), rk) - t0 select case (trim(params%terminal_detail)) case ('brief') if (flow_mpi%rank == 0) then write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,i0,2x,a,a,2x,a,a,2x,a,a)') & 'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, & 'piter=', stats%pressure_iterations, 'out=', yn(do_output_step), & 'rst=', yn(is_restart_write_step(params, step)), & 'rad=', yn(is_radiation_update_step(params, step)) end if case ('health') call compute_terminal_health(mesh, flow_mpi, params, fields, transport, energy, & terminal_projection_metric, terminal_ke, terminal_tmin, terminal_tmax) if (flow_mpi%rank == 0) then if (params%enable_variable_density) then if (params%enable_energy) then write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es10.3,2x,a,i0,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,a,2x,a,a,2x,a,a)') & 'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, & 'cfl=', stats%cfl, 'piter=', stats%pressure_iterations, & 'lm_res_proj=', terminal_projection_metric, 'Tmin=', terminal_tmin, & 'Tmax=', terminal_tmax, 'KE=', terminal_ke, & 'out=', yn(do_output_step), 'rst=', yn(is_restart_write_step(params, step)), & 'rad=', yn(is_radiation_update_step(params, step)) else write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es10.3,2x,a,i0,2x,a,es12.5,2x,a,es12.5,2x,a,a,2x,a,a,2x,a,a)') & 'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, & 'cfl=', stats%cfl, 'piter=', stats%pressure_iterations, & 'lm_res_proj=', terminal_projection_metric, 'KE=', terminal_ke, & 'out=', yn(do_output_step), 'rst=', yn(is_restart_write_step(params, step)), & 'rad=', yn(is_radiation_update_step(params, step)) end if else if (params%enable_energy) then write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es10.3,2x,a,i0,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,a,2x,a,a,2x,a,a)') & 'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, & 'cfl=', stats%cfl, 'piter=', stats%pressure_iterations, & 'max_div=', terminal_projection_metric, 'Tmin=', terminal_tmin, & 'Tmax=', terminal_tmax, 'KE=', terminal_ke, & 'out=', yn(do_output_step), 'rst=', yn(is_restart_write_step(params, step)), & 'rad=', yn(is_radiation_update_step(params, step)) else write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es10.3,2x,a,i0,2x,a,es12.5,2x,a,es12.5,2x,a,a,2x,a,a,2x,a,a)') & 'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, & 'cfl=', stats%cfl, 'piter=', stats%pressure_iterations, & 'max_div=', terminal_projection_metric, 'KE=', terminal_ke, & 'out=', yn(do_output_step), 'rst=', yn(is_restart_write_step(params, step)), & 'rad=', yn(is_radiation_update_step(params, step)) end if end if end if case ('full') if (params%enable_variable_density) then call compute_variable_density_console_diagnostics(mesh, flow_mpi, params, fields, transport, time, & console_max_divu, console_sproj_max, & console_lm_res_projection, console_lm_res_current, & console_mass_balance_defect) end if if (flow_mpi%rank == 0) then if (params%enable_variable_density) then write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,i0,2x,a,es12.5)') & 'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, 'max_divu=', console_max_divu, & 'lm_res_proj=', console_lm_res_projection, 'lm_res_cur=', console_lm_res_current, & 'Sproj_max=', console_sproj_max, 'mass_def=', console_mass_balance_defect, & 'piter=', stats%pressure_iterations, '|U|max=', stats%max_velocity else write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,i0,2x,a,es12.5)') & 'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, 'max_div=', stats%max_divergence, & 'cfl=', stats%cfl, 'ke=', stats%kinetic_energy, 'piter=', stats%pressure_iterations, '|U|max=', stats%max_velocity end if end if end select end if end do ! 8. Clean up and finalize simulation. call finalize_bc_set(bc) call finalize_fields(fields) if (params%enable_species) then call finalize_species(species) end if if (params%enable_energy) then call finalize_energy(energy) end if call finalize_transport(transport) call mesh_finalize(mesh) call profiler_stop('Total_Simulation') call profiler_report(flow_mpi%comm, flow_mpi%rank, flow_mpi%nprocs) if (params%write_cantera_cache_stats) then call write_chemistry_screening_stats(flow_mpi, params) call write_cantera_cache_stats(flow_mpi) end if call radiation_finalize(rad_context) call radiation_mpi_finalize(rad_mpi) call flow_mpi_finalize(flow_mpi) call mpi_flow_shutdown() contains subroutine save_step_backup(fields, species, energy, transport, stats, time, dt, chemistry_accum_dt, backup) type(flow_fields_t), intent(in) :: fields type(species_fields_t), intent(in) :: species type(energy_fields_t), intent(in) :: energy type(transport_properties_t), intent(in) :: transport type(solver_stats_t), intent(in) :: stats real(rk), intent(in) :: time, dt, chemistry_accum_dt type(step_backup_t), intent(inout) :: backup backup%fields = fields backup%species = species backup%energy = energy backup%transport = transport backup%stats = stats backup%time = time backup%dt = dt backup%chemistry_accum_dt = chemistry_accum_dt backup%valid = .true. end subroutine save_step_backup subroutine restore_step_backup(backup, fields, species, energy, transport, stats, time, dt, chemistry_accum_dt) type(step_backup_t), intent(in) :: backup type(flow_fields_t), intent(inout) :: fields type(species_fields_t), intent(inout) :: species type(energy_fields_t), intent(inout) :: energy type(transport_properties_t), intent(inout) :: transport type(solver_stats_t), intent(inout) :: stats real(rk), intent(inout) :: time, dt, chemistry_accum_dt if (.not. backup%valid) call fatal_error('main', 'attempted to restore an invalid step backup') fields = backup%fields species = backup%species energy = backup%energy transport = backup%transport stats = backup%stats time = backup%time dt = backup%dt chemistry_accum_dt = backup%chemistry_accum_dt end subroutine restore_step_backup subroutine compute_step_health(mesh, flow, params, fields, transport, kinetic_energy, cfl_value, max_velocity) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(case_params_t), intent(in) :: params type(flow_fields_t), intent(in) :: fields type(transport_properties_t), intent(in) :: transport real(rk), intent(out) :: kinetic_energy real(rk), intent(out) :: cfl_value real(rk), intent(out) :: max_velocity integer :: c, lf, f, ierr real(rk) :: rho_cell, vmag, local_ke, global_ke real(rk) :: local_umax, global_umax real(rk) :: local_cfl_rate, global_cfl_rate, cell_outward_flux local_ke = zero local_umax = zero local_cfl_rate = zero do c = flow%first_cell, flow%last_cell rho_cell = params%rho if (params%enable_variable_density .and. allocated(transport%rho)) then if (c <= size(transport%rho)) rho_cell = max(transport%rho(c), tiny_safe) end if vmag = sqrt(dot_product(fields%u(:, c), fields%u(:, c))) local_umax = max(local_umax, vmag) local_ke = local_ke + 0.5_rk * rho_cell * mesh%cells(c)%volume * vmag * vmag if (allocated(fields%face_flux)) then cell_outward_flux = zero do lf = 1, mesh%ncell_faces(c) f = mesh%cell_faces(lf, c) cell_outward_flux = cell_outward_flux + abs(fields%face_flux(f)) end do if (mesh%cells(c)%volume > tiny_safe) then local_cfl_rate = max(local_cfl_rate, 0.5_rk * cell_outward_flux / mesh%cells(c)%volume) end if end if end do call MPI_Allreduce(local_ke, global_ke, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing step kinetic energy') call MPI_Allreduce(local_umax, global_umax, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing step max velocity') call MPI_Allreduce(local_cfl_rate, global_cfl_rate, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing step CFL rate') kinetic_energy = global_ke max_velocity = global_umax cfl_value = global_cfl_rate * params%dt end subroutine compute_step_health !> Evaluate multi-guard physical and numerical safety gates for step rejection and retry. !! !! Evaluates multiple safety constraints over the trial timestep. If any limits are !! violated (e.g. NaN/Inf presence, chemistry increments, density vacuum floors, !! temperature excursions, or PCG linear solver convergence bounds), the trial step !! is rejected, triggering a temporal step reduction and retry. !! !! @param params Input case parameters and user namelist settings. !! @param stats Cumulative solver and PCG convergence metrics. !! @param max_dT_chem Maximum temperature increment during the chemistry advance [K]. !! @param max_dY_chem Maximum species mass fraction increment during the chemistry advance. !! @param max_rel_drho_chem Maximum relative density change during the chemistry advance. !! @param ke_before Kinetic energy of the domain before the trial timestep [J]. !! @param ke_after Trial domain kinetic energy after the trial timestep [J]. !! @param cfl_after Peak domain CFL value computed at the trial state. !! @param fields Struct containing current flow fields (u, p). !! @param energy Struct containing current thermodynamic/energy fields (T, h). !! @param transport Struct containing current transport/thermo properties (rho, mu). !! @param reject Output logical indicating if the trial step must be rejected. !! @param reason Output character string describing the first tripped safety gate. subroutine evaluate_step_rejection(params, stats, max_dT_chem, max_dY_chem, max_rel_drho_chem, & ke_before, ke_after, cfl_after, fields, energy, transport, & reject, reason) use, intrinsic :: ieee_arithmetic, only: ieee_is_nan type(case_params_t), intent(in) :: params type(solver_stats_t), intent(in) :: stats real(rk), intent(in) :: max_dT_chem, max_dY_chem, max_rel_drho_chem real(rk), intent(in) :: ke_before, ke_after, cfl_after type(flow_fields_t), intent(in) :: fields type(energy_fields_t), intent(in) :: energy type(transport_properties_t), intent(in) :: transport logical, intent(out) :: reject character(len=*), intent(out) :: reason real(rk) :: ke_growth, cfl_limit reject = .false. reason = 'accepted' if (.not. params%enable_step_rejection) return ! Option B.1: NaN and Infinity Sanitation if (params%reject_on_nan) then if (any(ieee_is_nan(fields%u)) .or. any(abs(fields%u) > 1.0e8_rk)) then reject = .true. reason = 'NaN or Inf detected in velocity field u' return end if if (any(ieee_is_nan(fields%p)) .or. any(abs(fields%p) > 1.0e12_rk)) then reject = .true. reason = 'NaN or Inf detected in pressure field p' return end if if (allocated(energy%T)) then if (any(ieee_is_nan(energy%T)) .or. any(energy%T > 1.0e5_rk)) then reject = .true. reason = 'NaN or Inf detected in temperature field T' return end if end if if (allocated(transport%rho)) then if (any(ieee_is_nan(transport%rho))) then reject = .true. reason = 'NaN detected in density field rho' return end if end if end if ! Option B.2: Chemistry Limits if (params%chemistry_max_dT_per_step > zero .and. max_dT_chem > params%chemistry_max_dT_per_step) then reject = .true. reason = 'chemistry dT exceeded committed limit' return end if if (params%chemistry_max_dY_per_step > zero .and. max_dY_chem > params%chemistry_max_dY_per_step) then reject = .true. reason = 'chemistry dY exceeded committed limit' return end if if (params%chemistry_max_rel_rho_change_per_step > zero .and. & max_rel_drho_chem > params%chemistry_max_rel_rho_change_per_step) then reject = .true. reason = 'chemistry density change exceeded committed limit' return end if ! Option B.3: CFL Overshoot if (params%max_cfl_overshoot_factor > zero) then cfl_limit = params%max_cfl * params%max_cfl_overshoot_factor if (cfl_after > cfl_limit) then reject = .true. reason = 'post-step CFL overshoot' return end if end if ! Option B.4: Kinetic Energy Growth if (params%max_KE_growth_per_step > zero .and. ke_before > params%ke_reject_floor) then ke_growth = ke_after / ke_before if (ke_growth > params%max_KE_growth_per_step) then reject = .true. reason = 'kinetic energy growth exceeded limit' return end if end if ! Option B.5: PCG Pressure Iterations Ceiling if (params%reject_on_pressure_max_iter .and. stats%pressure_iterations >= params%pressure_max_iter) then reject = .true. reason = 'pressure solve reached max iterations' return end if ! Option B.6: PCG Convergence Health Audit if (stats%pressure_residual > params%max_pressure_residual_limit .and. & stats%pressure_residual_abs > params%pressure_abs_tol) then reject = .true. reason = 'pressure solver residual exceeds stability tolerance' return end if ! Option B.7: Density Floor/Vacuum Guard if (allocated(transport%rho)) then if (minval(transport%rho) < params%rho * params%min_density_ratio_limit) then reject = .true. reason = 'density vacuum limit violated' return end if end if ! Option B.8: Temperature Floor and Ceiling Excursions if (allocated(energy%T)) then if (minval(energy%T) < params%min_temperature_limit) then reject = .true. reason = 'temperature fell below safety floor' return end if if (maxval(energy%T) > params%max_temperature_limit) then reject = .true. reason = 'temperature exceeded safety ceiling' return end if end if end subroutine evaluate_step_rejection function int_to_string(value) result(text) integer, intent(in) :: value character(len=32) :: text write(text,'(i0)') value end function int_to_string !> Cut the next dynamic timestep when chemistry produced a large state jump. !! This does not reject the current step; it prevents the next operator-split !! chemistry/projection update from repeating an overly large heat-release jump. subroutine apply_chemistry_timestep_control(flow, params, max_dT, max_dY, max_rel_drho, min_alpha) type(flow_mpi_t), intent(in) :: flow type(case_params_t), intent(inout) :: params real(rk), intent(in) :: max_dT, max_dY, max_rel_drho, min_alpha real(rk) :: factor, old_dt, threshold if (.not. params%use_dynamic_dt) return factor = one threshold = params%chemistry_max_dT_per_step if (threshold > zero .and. max_dT > threshold) then factor = min(factor, params%chemistry_dt_safety * threshold / max(max_dT, tiny(1.0_rk))) end if threshold = params%chemistry_max_dY_per_step if (threshold > zero .and. max_dY > threshold) then factor = min(factor, params%chemistry_dt_safety * threshold / max(max_dY, tiny(1.0_rk))) end if threshold = params%chemistry_max_rel_rho_change_per_step if (threshold > zero .and. max_rel_drho > threshold) then factor = min(factor, params%chemistry_dt_safety * threshold / max(max_rel_drho, tiny(1.0_rk))) end if factor = max(params%chemistry_dt_min_factor, min(one, factor)) if (factor < one) then old_dt = params%dt params%dt = max(params%min_dt, min(params%max_dt, params%dt * factor)) if (flow%rank == 0) then write(output_unit,'(a,es11.4,a,es11.4,a,es11.4,a,es11.4,a,es11.4,a,es11.4,a,es11.4)') & 'chemistry dt control: dt ', old_dt, ' -> ', params%dt, & ' max_dT=', max_dT, ' max_dY=', max_dY, & ' max_rel_drho=', max_rel_drho, ' min_alpha=', min_alpha end if end if end subroutine apply_chemistry_timestep_control !> Apply a one-time spherical ignition kernel during initialization. !! !! The current first version supports one spherical overwrite kernel. It is !! skipped automatically for restart runs so restart state is not modified by !! fresh-run initialization controls. subroutine apply_ignition_kernel(mesh, flow, params, species, energy) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(case_params_t), intent(in) :: params type(species_fields_t), intent(inout) :: species type(energy_fields_t), intent(inout) :: energy integer :: c, affected_cells real(rk) :: delta(3), radius2 real(rk) :: kernel_Y(max_species) logical :: set_composition if (.not. params%enable_ignition_kernel) return if (.not. params%enable_energy) return if (.not. energy%initialized) return affected_cells = 0 radius2 = params%ignition_kernel_radius * params%ignition_kernel_radius set_composition = .false. kernel_Y = zero if (params%enable_species .and. species%nspecies > 0 .and. & composition_string_is_set(params%ignition_kernel_composition)) then call parse_named_composition(params%ignition_kernel_composition, species%names, species%nspecies, & kernel_Y, 'ignition_kernel_composition', & normalize=.true., require_sum_positive=.true.) set_composition = .true. end if do c = 1, mesh%ncells delta = mesh%cells(c)%center - params%ignition_kernel_center if (dot_product(delta, delta) <= radius2) then energy%T(c) = params%ignition_kernel_T if (set_composition) species%Y(1:species%nspecies, c) = kernel_Y(1:species%nspecies) affected_cells = affected_cells + 1 end if end do if (set_composition) then species%Y_old = species%Y species%rhs_old = zero species%rhs_old_valid = .false. end if if (flow%rank == 0) then write(output_unit,'(a,i0,a,es12.5,a,3(1x,es12.5))') & 'Ignition kernel initialized: cells=', affected_cells, & ' T=', params%ignition_kernel_T, ' center=', params%ignition_kernel_center end if end subroutine apply_ignition_kernel character(len=1) function yn(flag) result(mark) logical, intent(in) :: flag if (flag) then mark = 'Y' else mark = 'N' end if end function yn logical function is_restart_write_step(params, step) result(write_now) type(case_params_t), intent(in) :: params integer, intent(in) :: step write_now = .false. if (.not. params%write_restart) return if (params%restart_interval <= 0) return write_now = (mod(step, params%restart_interval) == 0) .or. (step == params%nsteps) end function is_restart_write_step logical function is_radiation_update_step(params, step) result(update_now) type(case_params_t), intent(in) :: params integer, intent(in) :: step update_now = .false. if (.not. params%enable_energy) return if (.not. params%enable_radiation) return if (params%radiation_update_interval <= 0) return update_now = (mod(step - 1, params%radiation_update_interval) == 0) .or. (step == 1) end function is_radiation_update_step subroutine compute_terminal_health(mesh, flow, params, fields, transport, energy, & projection_metric, kinetic_energy, tmin, tmax) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(case_params_t), intent(in) :: params type(flow_fields_t), intent(in) :: fields type(transport_properties_t), intent(in) :: transport type(energy_fields_t), intent(in) :: energy real(rk), intent(out) :: projection_metric real(rk), intent(out) :: kinetic_energy real(rk), intent(out) :: tmin real(rk), intent(out) :: tmax integer :: c, lf, fid, ierr real(rk) :: vol, divu, flux_out, target, rho_cell real(rk) :: local_metric, local_ke, global_ke real(rk) :: local_temp_vec(2), global_temp_vec(2) logical :: have_temperature local_metric = zero local_ke = zero local_temp_vec = [huge(one), huge(one)] have_temperature = params%enable_energy .and. allocated(energy%T) do c = 1, mesh%ncells if (.not. flow%owned(c)) cycle vol = mesh%cells(c)%volume if (vol <= zero) cycle divu = zero if (allocated(fields%face_flux)) then do lf = 1, mesh%ncell_faces(c) fid = mesh%cell_faces(lf, c) if (mesh%faces(fid)%owner == c) then flux_out = fields%face_flux(fid) else flux_out = -fields%face_flux(fid) end if divu = divu + flux_out / vol end do end if target = zero if (params%enable_variable_density .and. allocated(fields%projection_divergence_source)) then if (c <= size(fields%projection_divergence_source)) target = fields%projection_divergence_source(c) end if local_metric = max(local_metric, abs(divu - target)) rho_cell = params%rho if (params%enable_variable_density .and. allocated(transport%rho)) then if (c <= size(transport%rho)) rho_cell = max(transport%rho(c), tiny(one)) end if local_ke = local_ke + 0.5_rk * rho_cell * vol * dot_product(fields%u(:, c), fields%u(:, c)) if (have_temperature) then if (c <= size(energy%T)) then local_temp_vec(1) = min(local_temp_vec(1), energy%T(c)) local_temp_vec(2) = min(local_temp_vec(2), -energy%T(c)) end if end if end do call MPI_Allreduce(local_metric, projection_metric, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing terminal projection health') call MPI_Allreduce(local_ke, global_ke, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing terminal kinetic energy') kinetic_energy = global_ke if (have_temperature) then call MPI_Allreduce(local_temp_vec, global_temp_vec, 2, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing terminal temperature range') tmin = global_temp_vec(1) tmax = -global_temp_vec(2) else tmin = zero tmax = zero end if end subroutine compute_terminal_health !> Write a concise, reproducible summary of physics and numerical settings. subroutine write_run_config_summary(mesh, flow, rad, params, case_file) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(radiation_mpi_t), intent(in) :: rad type(case_params_t), intent(in) :: params character(len=*), intent(in) :: case_file integer :: unit_id, ios character(len=1024) :: filename if (flow%rank /= 0) return call write_run_config_summary_to_unit(output_unit, mesh, flow, rad, params, case_file) filename = trim(params%output_dir) // '/run_config_summary.txt' open(newunit=unit_id, file=trim(filename), status='replace', action='write', iostat=ios) if (ios /= 0) then write(output_unit,'(a,a)') 'warning: could not write run configuration summary: ', trim(filename) return end if call write_run_config_summary_to_unit(unit_id, mesh, flow, rad, params, case_file) close(unit_id) end subroutine write_run_config_summary subroutine write_run_config_summary_to_unit(unit_id, mesh, flow, rad, params, case_file) integer, intent(in) :: unit_id type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(radiation_mpi_t), intent(in) :: rad type(case_params_t), intent(in) :: params character(len=*), intent(in) :: case_file integer :: i write(unit_id,'(a)') '======================================================================' write(unit_id,'(a)') 'RUN CONFIGURATION SUMMARY' write(unit_id,'(a)') '======================================================================' write(unit_id,'(a,a)') 'case: ', trim(case_file) write(unit_id,'(a,a)') 'mesh_dir: ', trim(params%mesh_dir) write(unit_id,'(a,i0)') 'cells: ', mesh%ncells write(unit_id,'(a,i0)') 'faces: ', mesh%nfaces write(unit_id,'(a,i0)') 'flow/chemistry MPI ranks: ', flow%nprocs write(unit_id,'(a,i0)') 'radiation MPI ranks: ', rad%nprocs write(unit_id,'(a,i0,a,es12.5,a,i0,a,i0,a,l1)') 'time: nsteps=', params%nsteps, ' dt=', params%dt, & ' output_interval=', params%output_interval, ' terminal_interval=', params%terminal_interval, & ' terminal_detail='//trim(params%terminal_detail)//' dynamic_dt=', params%use_dynamic_dt write(unit_id,'(a,es12.5,a,es12.5,a,es12.5)') 'time limits: max_dt=', params%max_dt, & ' min_dt=', params%min_dt, ' dt_growth_limit=', params%dt_growth_limit write(unit_id,'(a,l1)') ' step rejection enabled : ', params%enable_step_rejection write(unit_id,'(a,i0)') ' max step retries : ', params%max_step_retries write(unit_id,'(a,es12.5)') ' step reject cut factor : ', params%step_reject_cut_factor write(unit_id,'(a,es12.5)') ' step reject minimum dt : ', params%step_reject_min_dt write(unit_id,'(a,es12.5)') ' max KE growth per step : ', params%max_KE_growth_per_step write(unit_id,'(a,es12.5)') ' KE reject floor : ', params%ke_reject_floor write(unit_id,'(a,es12.5)') ' max CFL overshoot factor : ', params%max_cfl_overshoot_factor write(unit_id,'(a,l1)') ' reject on pressure max iter : ', params%reject_on_pressure_max_iter write(unit_id,'(a,l1)') ' stop on min-dt failure : ', params%stop_on_min_dt_failure write(unit_id,'(a)') '' write(unit_id,'(a)') 'Numerical schemes:' write(unit_id,'(a,a)') ' legacy convection_scheme : ', trim(params%convection_scheme) write(unit_id,'(a,a)') ' momentum convection : ', trim(params%momentum_convection_scheme) write(unit_id,'(a,a)') ' species convection : ', trim(params%species_convection_scheme) write(unit_id,'(a,a)') ' energy convection : ', trim(params%energy_convection_scheme) write(unit_id,'(a,a)') ' scalar limiter : ', trim(params%scalar_limiter) write(unit_id,'(a,a)') ' species time scheme : ', trim(params%species_time_scheme) write(unit_id,'(a,a)') ' energy time scheme : ', trim(params%energy_time_scheme) write(unit_id,'(a,i0)') ' pressure max iterations : ', params%pressure_max_iter write(unit_id,'(a,es12.5)') ' pressure tolerance : ', params%pressure_tol write(unit_id,'(a)') '' write(unit_id,'(a)') 'Flow and thermodynamics:' write(unit_id,'(a,l1)') ' variable density : ', params%enable_variable_density write(unit_id,'(a,a)') ' density EOS : ', trim(params%density_eos) write(unit_id,'(a,l1)') ' variable viscosity : ', params%enable_variable_nu write(unit_id,'(a,i0)') ' transport update interval : ', params%transport_update_interval write(unit_id,'(a,l1)') ' energy equation : ', params%enable_energy write(unit_id,'(a,l1)') ' Cantera thermo : ', params%enable_cantera_thermo write(unit_id,'(a,l1)') ' species enthalpy diffusion : ', params%enable_species_enthalpy_diffusion write(unit_id,'(a,l1)') ' upwind dilatation density : ', params%upwind_dilatation_density write(unit_id,'(a,es12.5)') ' background pressure [Pa] : ', params%background_press write(unit_id,'(a,es12.5)') ' background temperature [K] : ', params%background_temp write(unit_id,'(a)') '' write(unit_id,'(a)') 'Chemistry timestep/source controls:' write(unit_id,'(a,es12.5)') ' chemistry max dT per step [K] : ', params%chemistry_max_dT_per_step write(unit_id,'(a,es12.5)') ' chemistry max rel rho change : ', params%chemistry_max_rel_rho_change_per_step write(unit_id,'(a,es12.5)') ' chemistry max dY per step : ', params%chemistry_max_dY_per_step write(unit_id,'(a,es12.5)') ' chemistry dt safety : ', params%chemistry_dt_safety write(unit_id,'(a,l1)') ' chemistry source limiter : ', params%chemistry_limit_source_update write(unit_id,'(a,es12.5)') ' chemistry source relaxation : ', params%chemistry_source_relaxation write(unit_id,'(a)') '' write(unit_id,'(a)') 'Species and chemistry:' write(unit_id,'(a,l1)') ' species transport : ', params%enable_species write(unit_id,'(a,i0)') ' transported species count : ', params%nspecies write(unit_id,'(a,l1)') ' reactions enabled : ', params%enable_reactions write(unit_id,'(a,a)') ' Cantera mechanism : ', trim(params%cantera_mech_file) if (len_trim(params%cantera_phase_name) > 0) then write(unit_id,'(a,a)') ' Cantera phase : ', trim(params%cantera_phase_name) else write(unit_id,'(a)') ' Cantera phase : ' end if write(unit_id,'(a,i0)') ' Cantera reaction count : ', params%cantera_reaction_count write(unit_id,'(a,i0)') ' chemistry update interval : ', params%chemistry_update_interval write(unit_id,'(a,es12.5)') ' chemistry rtol : ', params%chemistry_rtol write(unit_id,'(a,es12.5)') ' chemistry atol : ', params%chemistry_atol write(unit_id,'(a,i0)') ' chemistry max steps : ', params%chemistry_max_steps write(unit_id,'(a,l1)') ' chemistry energy enabled : ', params%chemistry_energy_enabled write(unit_id,'(a,es12.5)') ' chemistry T cutoff [K] : ', params%chemistry_temperature_cutoff write(unit_id,'(a,es12.5)') ' chemistry min reactive Y : ', params%chemistry_min_reactive_mass_fraction write(unit_id,'(a,l1)') ' chemistry_active_screening : ', & params%chemistry_active_species_threshold > zero .and. params%chemistry_n_active_species > 0 write(unit_id,'(a,i0)') ' chemistry_n_active_species : ', params%chemistry_n_active_species write(unit_id,'(a,es12.5)') ' chemistry_active_species_threshold : ', params%chemistry_active_species_threshold if (params%chemistry_n_active_species > 0) then do i = 1, params%chemistry_n_active_species if (len_trim(params%chemistry_active_species_name(i)) > 0) then write(unit_id,'(a,i0,a,a)') ' chemistry_active_species_name(', i, ') : ', & trim(params%chemistry_active_species_name(i)) end if end do end if if (len_trim(params%initial_composition) > 0) then write(unit_id,'(a,a)') ' initial composition : ', trim(params%initial_composition) end if call write_ignition_config_summary_to_unit(unit_id, mesh, params) call write_boundary_config_summary_to_unit(unit_id, params) write(unit_id,'(a)') '' write(unit_id,'(a)') 'Radiation coupling:' write(unit_id,'(a,l1)') ' radiation enabled : ', params%enable_radiation write(unit_id,'(a,a)') ' radiation pressure source : ', trim(params%radiation_pressure_source) write(unit_id,'(a,a)') ' radiation model : ', trim(params%radiation_source_model) write(unit_id,'(a,i0)') ' radiation update interval : ', params%radiation_update_interval write(unit_id,'(a,i0)') ' n wavenumbers : ', params%radiation_n_wavenumbers write(unit_id,'(a,i0,a,i0)') ' local wavenumber range : ', rad%first_wavenumber, '-', rad%last_wavenumber write(unit_id,'(a,i0)') ' selected radiation species : ', params%radiation_n_species write(unit_id,'(a,l1)') ' radiation debug checks : ', params%radiation_debug write(unit_id,'(a,l1)') ' radiation diagnostics : ', params%write_radiation_diagnostics write(unit_id,'(a)') '' write(unit_id,'(a)') 'Restart:' write(unit_id,'(a,l1)') ' restart from file : ', params%restart_from_file write(unit_id,'(a,a)') ' restart file : ', trim(params%restart_file) write(unit_id,'(a,l1)') ' write restart : ', params%write_restart write(unit_id,'(a,i0)') ' restart interval : ', params%restart_interval write(unit_id,'(a,a)') ' restart output dir : ', trim(params%restart_output_dir) write(unit_id,'(a,i0)') ' restart keep : ', params%restart_keep write(unit_id,'(a,l1)') ' restart use file dt : ', params%restart_use_file_dt write(unit_id,'(a,l1)') ' restart nsteps additional : ', params%restart_nsteps_are_additional write(unit_id,'(a)') '======================================================================' end subroutine write_run_config_summary_to_unit !> Write ignition-kernel settings to the run configuration summary. subroutine write_ignition_config_summary_to_unit(unit_id, mesh, params) integer, intent(in) :: unit_id type(mesh_t), intent(in) :: mesh type(case_params_t), intent(in) :: params integer :: c, affected_cells real(rk) :: delta(3), radius2 write(unit_id,'(a)') '' write(unit_id,'(a)') 'Ignition initialization:' write(unit_id,'(a,l1)') ' enable_ignition_kernel : ', params%enable_ignition_kernel if (.not. params%enable_ignition_kernel) return affected_cells = 0 radius2 = params%ignition_kernel_radius * params%ignition_kernel_radius do c = 1, mesh%ncells delta = mesh%cells(c)%center - params%ignition_kernel_center if (dot_product(delta, delta) <= radius2) affected_cells = affected_cells + 1 end do write(unit_id,'(a,a)') ' ignition_kernel_shape : ', trim(params%ignition_kernel_shape) write(unit_id,'(a,3(1x,es12.5))') ' ignition_kernel_center : ', params%ignition_kernel_center write(unit_id,'(a,es12.5)') ' ignition_kernel_radius : ', params%ignition_kernel_radius write(unit_id,'(a,es12.5)') ' ignition_kernel_T [K] : ', params%ignition_kernel_T write(unit_id,'(a,a)') ' ignition_kernel_blend : ', trim(params%ignition_kernel_blend) write(unit_id,'(a,i0)') ' ignition_kernel_cells : ', affected_cells if (len_trim(params%ignition_kernel_composition) > 0) then write(unit_id,'(a,a)') ' ignition_kernel_composition : ', trim(params%ignition_kernel_composition) else write(unit_id,'(a)') ' ignition_kernel_composition : ' end if if (params%restart_from_file) then write(unit_id,'(a)') ' note : restart_from_file=.true.; fresh-run kernel is skipped' end if end subroutine write_ignition_config_summary_to_unit !> Write boundary-condition input settings to the run configuration summary. subroutine write_boundary_config_summary_to_unit(unit_id, params) integer, intent(in) :: unit_id type(case_params_t), intent(in) :: params integer :: p write(unit_id,'(a)') '' write(unit_id,'(a)') 'Boundary conditions:' write(unit_id,'(a,i0)') ' number of patches : ', params%n_patches write(unit_id,'(a)') ' mass-flux velocity convention:' write(unit_id,'(a)') ' mass_flux/signed_mass_flux : signed rho*dot(U,n_out) [kg/m^2/s]' write(unit_id,'(a)') ' inlet_mass_flux : positive magnitude directed into the domain' write(unit_id,'(a)') ' outlet_mass_flux : positive magnitude directed out of the domain' if (params%n_patches <= 0) return do p = 1, params%n_patches write(unit_id,'(a,i0,a,a)') ' patch[', p, '] ', trim(params%patch_name(p)) write(unit_id,'(a,a)') ' patch type : ', trim(params%patch_type(p)) write(unit_id,'(a,a)') ' velocity type : ', trim(params%patch_velocity_type(p)) write(unit_id,'(a,3(es12.5,1x))') ' velocity [m/s] : ', & params%patch_u(p), params%patch_v(p), params%patch_w(p) write(unit_id,'(a,es12.5)') ' mass flux [kg/m^2/s] : ', params%patch_mass_flux(p) write(unit_id,'(a,a)') ' pressure type : ', trim(params%patch_pressure_type(p)) write(unit_id,'(a,es12.5)') ' pressure value : ', params%patch_p(p) write(unit_id,'(a,es12.5)') ' pressure gradient : ', params%patch_dpdn(p) write(unit_id,'(a,a)') ' temperature type : ', trim(params%patch_temperature_type(p)) write(unit_id,'(a,es12.5)') ' temperature [K] : ', params%patch_T(p) write(unit_id,'(a,a)') ' species type : ', trim(params%patch_species_type(p)) if (len_trim(params%patch_composition(p)) > 0) then write(unit_id,'(a,a)') ' composition : ', trim(params%patch_composition(p)) end if end do end subroutine write_boundary_config_summary_to_unit !> Computes global integral quantities for diagnostic reporting. !! !! Performs Allreduce operations to find the maximum velocity magnitude, !! total system mass, and minimum species mass fraction. subroutine compute_global_observables(mesh, flow, fields, species, transport, stats) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(flow_fields_t), intent(in) :: fields type(species_fields_t), intent(in) :: species type(transport_properties_t), intent(in) :: transport type(solver_stats_t), intent(inout) :: stats integer :: c real(rk) :: local_max_vel, local_min_y, local_total_mass real(rk) :: global_max_vel, global_min_y, global_total_mass real(rk) :: u, v, w, vel_mag local_max_vel = zero local_min_y = huge(one) local_total_mass = zero do c = 1, mesh%ncells if (.not. flow%owned(c)) cycle u = fields%u(1, c) v = fields%u(2, c) w = fields%u(3, c) vel_mag = sqrt(u**2 + v**2 + w**2) if (vel_mag > local_max_vel) local_max_vel = vel_mag local_total_mass = local_total_mass + (transport%rho(c) * mesh%cells(c)%volume) if (species%nspecies > 0) then local_min_y = min(local_min_y, minval(species%Y(:, c))) end if end do call MPI_Allreduce(local_max_vel, global_max_vel, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm) call MPI_Allreduce(local_total_mass, global_total_mass, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm) if (species%nspecies > 0) then call MPI_Allreduce(local_min_y, global_min_y, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm) else global_min_y = zero end if stats%max_velocity = global_max_vel stats%total_mass = global_total_mass stats%min_species_y = global_min_y end subroutine compute_global_observables !> Compute concise variable-density low-Mach diagnostics for the terminal. !! !! In variable-density mode raw div(u) is not a projection-error metric !! because the low-Mach target is div(u)=S. This helper reports both raw !! div(u) and residuals against the source actually used by projection. subroutine compute_variable_density_console_diagnostics(mesh, flow, params, fields, transport, time, & max_divu, max_s_projection, & max_lm_res_projection, max_lm_res_current, & mass_balance_defect) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(case_params_t), intent(in) :: params type(flow_fields_t), intent(in) :: fields type(transport_properties_t), intent(in) :: transport real(rk), intent(in) :: time real(rk), intent(out) :: max_divu real(rk), intent(out) :: max_s_projection real(rk), intent(out) :: max_lm_res_projection real(rk), intent(out) :: max_lm_res_current real(rk), intent(out) :: mass_balance_defect integer :: c, lf, fid, f, owner, nb, ierr real(rk) :: vol, divu, flux_out, s_current, s_projection real(rk) :: local_max_divu, local_max_s_projection real(rk) :: local_max_lm_res_projection, local_max_lm_res_current real(rk) :: local_mass, global_mass, local_net_mass_flux, net_mass_flux real(rk) :: dt_since_previous, mass_rate real(rk) :: local_max_vec(4), global_max_vec(4) real(rk) :: local_sum_vec(2), global_sum_vec(2) logical, save :: initialized = .false. real(rk), save :: previous_mass = 0.0_rk real(rk), save :: previous_time = 0.0_rk max_divu = 0.0_rk max_s_projection = 0.0_rk max_lm_res_projection = 0.0_rk max_lm_res_current = 0.0_rk mass_balance_defect = 0.0_rk if (.not. params%enable_variable_density) return if (.not. allocated(fields%face_flux)) return local_max_divu = 0.0_rk local_max_s_projection = 0.0_rk local_max_lm_res_projection = 0.0_rk local_max_lm_res_current = 0.0_rk local_mass = 0.0_rk local_net_mass_flux = 0.0_rk do c = 1, mesh%ncells if (.not. flow%owned(c)) cycle vol = mesh%cells(c)%volume if (vol <= 0.0_rk) cycle divu = 0.0_rk do lf = 1, mesh%ncell_faces(c) fid = mesh%cell_faces(lf, c) if (mesh%faces(fid)%owner == c) then flux_out = fields%face_flux(fid) else flux_out = -fields%face_flux(fid) end if divu = divu + flux_out / vol end do s_current = 0.0_rk if (allocated(fields%divergence_source)) then if (size(fields%divergence_source) >= c) s_current = fields%divergence_source(c) end if s_projection = s_current if (allocated(fields%projection_divergence_source)) then if (size(fields%projection_divergence_source) >= c) s_projection = fields%projection_divergence_source(c) end if local_max_divu = max(local_max_divu, abs(divu)) local_max_s_projection = max(local_max_s_projection, abs(s_projection)) local_max_lm_res_projection = max(local_max_lm_res_projection, abs(divu - s_projection)) local_max_lm_res_current = max(local_max_lm_res_current, abs(divu - s_current)) if (allocated(transport%rho)) then if (size(transport%rho) >= c) local_mass = local_mass + transport%rho(c) * vol end if end do if (allocated(fields%mass_flux)) then do f = 1, mesh%nfaces owner = mesh%faces(f)%owner nb = mesh%faces(f)%neighbor if (nb <= 0) nb = mesh%faces(f)%periodic_neighbor if (nb <= 0) then if (owner >= 1 .and. owner <= mesh%ncells) then if (flow%owned(owner)) local_net_mass_flux = local_net_mass_flux + fields%mass_flux(f) end if end if end do else if (allocated(fields%face_flux) .and. allocated(transport%rho)) then do f = 1, mesh%nfaces owner = mesh%faces(f)%owner nb = mesh%faces(f)%neighbor if (nb <= 0) nb = mesh%faces(f)%periodic_neighbor if (nb <= 0) then if (owner >= 1 .and. owner <= mesh%ncells) then if (flow%owned(owner)) local_net_mass_flux = local_net_mass_flux + transport%rho(owner) * fields%face_flux(f) end if end if end do end if local_max_vec = [local_max_divu, local_max_s_projection, & local_max_lm_res_projection, local_max_lm_res_current] call MPI_Allreduce(local_max_vec, global_max_vec, 4, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing variable-density max diagnostics') max_divu = global_max_vec(1) max_s_projection = global_max_vec(2) max_lm_res_projection = global_max_vec(3) max_lm_res_current = global_max_vec(4) local_sum_vec = [local_mass, local_net_mass_flux] call MPI_Allreduce(local_sum_vec, global_sum_vec, 2, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing variable-density sum diagnostics') global_mass = global_sum_vec(1) net_mass_flux = global_sum_vec(2) if (initialized) then dt_since_previous = time - previous_time if (dt_since_previous > tiny(1.0_rk)) then mass_rate = (global_mass - previous_mass) / dt_since_previous mass_balance_defect = mass_rate + net_mass_flux else mass_balance_defect = 0.0_rk end if else mass_balance_defect = 0.0_rk initialized = .true. end if previous_mass = global_mass previous_time = time end subroutine compute_variable_density_console_diagnostics !> Parses command line arguments to find the case configuration file. !! !! If no argument is provided, defaults to `case.nml`. subroutine get_case_filename(filename) character(len=*), intent(out) :: filename integer :: argc argc = command_argument_count() if (argc >= 1) then call get_command_argument(1, filename) else filename = 'case.nml' end if if (len_trim(filename) == 0) call fatal_error('main', 'empty case filename') end subroutine get_case_filename end program lowmach_react_hex