!> Energy field storage, diagnostics, and sensible-enthalpy transport. !! !! This module owns the enthalpy/temperature/radiation-source fields used by !! the current constant-density energy path. The transported thermodynamic !! state is `h`; temperature is recovered from `h`, composition `Y`, and the !! background thermodynamic pressure `params%background_press`. !! !! In Cantera mode the stored sensible enthalpy is !! `h_abs(T,Y,p0) - h_abs(T_ref,Y,p0)`, where `p0` is the background pressure !! and `T_ref` is `params%energy_reference_T`. After species transport changes !! `Y`, Option A is enforced: preserve transported `h` and recover !! `T(h,Y,p0)`. Do not rebuild `h` from the old temperature and new !! composition. !! !! Density ownership depends on mode. In constant-density mode !! `transport%rho = params%rho` and `energy%rho_thermo` is diagnostic. In !! guarded variable-density mode with `density_eos="cantera"`, !! `transport%rho <- energy%rho_thermo` after Cantera thermo sync, so density !! comes from the selected Cantera phase and its YAML EOS. Heat conduction uses !! `grad(T)`, not `grad(h)`. `qrad` is a volumetric source hook with `qrad > 0` !! adding energy to the gas; external radiation coupling is future work. !! !! Reactions and reaction heat release are not implemented in the supported !! path. The optional species-enthalpy diffusion correction !! `-div(sum_k h_k J_k)` is controlled by !! `enable_species_enthalpy_diffusion`. module mod_energy use mpi_f08 use iso_c_binding, only : c_char, c_double, c_int use mod_precision, only : rk, zero, one, tiny_safe, fatal_error, output_unit, lowercase, name_len use mod_mesh, only : mesh_t use mod_input, only : case_params_t, max_species use mod_mpi_flow, only : flow_mpi_t, flow_exchange_cell_scalar, flow_exchange_cell_scalars, flow_exchange_cell_matrix use mod_cantera, only : transport_properties_t use mod_reconstruction, only : compute_scalar_gradients, scalar_face_value use mod_bc, only : bc_set_t, bc_periodic, face_effective_neighbor, & patch_type_for_face, boundary_temperature, boundary_species use mod_flow_fields, only : flow_fields_t use mod_profiling, only : profiler_start, profiler_stop implicit none private !> Cell-centered energy variables. type, public :: energy_fields_t real(rk), allocatable :: T(:) !< Temperature [K]. real(rk), allocatable :: T_old(:) !< Previous temperature state [K]. real(rk), allocatable :: h(:) !< Transported mixture sensible enthalpy [J/kg]. real(rk), allocatable :: h_old(:) !< Previous enthalpy state [J/kg]. real(rk), allocatable :: rhs_old(:) !< Previous enthalpy transport RHS for optional AB2 time integration. logical :: rhs_old_valid = .false. !< True once rhs_old contains a previous transport RHS. real(rk), allocatable :: qrad(:) !< Volumetric source term [W/m^3]; positive values heat the gas. real(rk), allocatable :: cp(:) !< Mixture heat capacity at constant pressure [J/kg/K]. real(rk), allocatable :: lambda(:) !< Thermal conductivity [W/m/K]. real(rk), allocatable :: rho_thermo(:) !< Cantera thermodynamic density [kg/m^3]; diagnostic unless variable-density mode syncs it into `transport%rho`. real(rk), allocatable :: operator_consistent_rho_h(:) !< Cellwise update-consistent rho*h [J/m^3]. integer :: operator_consistent_rho_h_available = 0 !< 1 when cellwise update-consistent rho*h is valid. real(rk), allocatable :: species_enthalpy_diffusion(:) !< Specific enthalpy source from species diffusion [J/kg/s]. real(rk), allocatable :: chem_heat_release_rate(:) !< Approx. committed chemistry heat release rho*cp*dT/dt [W/m^3]. real(rk), allocatable :: chem_dT_dt(:) !< Committed chemistry temperature source [K/s]. real(rk), allocatable :: chem_dY_max_dt(:) !< Max committed chemistry species source magnitude [1/s]. real(rk), allocatable :: chem_rel_rho_change_dt(:) !< Relative chemistry density-change rate [1/s]. real(rk), allocatable :: chem_source_alpha(:) !< Chemistry limiter alpha committed in the last update [-]. real(rk), allocatable :: chem_active(:) !< 1 if chemistry was called in the last update, else 0. real(rk) :: last_conductive_boundary_flux_out = zero !< Last local outward conductive boundary flux [W]. integer :: last_conductive_boundary_flux_available = 0 !< 1 after conductive-boundary diagnostic is available. real(rk) :: cumulative_boundary_rho_h_advective_flux_out = zero !< Time-integrated outward advective rho*h boundary flux [J]. real(rk) :: cumulative_boundary_rho_h_conductive_flux_out = zero !< Time-integrated outward conductive boundary flux [J]. real(rk) :: cumulative_qrad_integral = zero !< Time-integrated qrad volume source [J]. real(rk) :: cumulative_rho_species_hdiff_integral = zero !< Time-integrated rho*species enthalpy diffusion source [J]. integer :: cumulative_energy_budget_available = 0 !< 1 after cumulative energy budget terms are available. real(rk) :: last_energy_update_delta_rate_integral = zero !< Last global d/dt integral(rho*h) from update [W]. real(rk) :: last_energy_update_rhs_integral = zero !< Last global energy RHS integral from update [W]. real(rk) :: last_energy_update_balance_defect = zero !< Last direct energy update closure defect [W]. real(rk) :: relative_last_energy_update_balance_defect = zero !< Relative direct update closure defect. real(rk) :: last_operator_consistent_rho_h_integral = zero !< Last global rho*h integral using update density/time level. real(rk) :: cumulative_energy_update_delta_integral = zero !< Time-integrated conservative energy-update delta [J]. real(rk) :: cumulative_energy_update_rhs_integral = zero !< Time-integrated conservative energy-update RHS [J]. logical :: initialized = .false. !< True after allocation/initialization. end type energy_fields_t public :: allocate_energy, initialize_energy, finalize_energy public :: refresh_energy_from_temperature public :: update_enthalpy_from_temperature_constant_cp public :: recover_temperature_constant_cp public :: recover_temperature_cantera public :: zero_radiation_source public :: advance_cantera_chemistry public :: advance_energy_transport public :: write_energy_diagnostics_header, write_energy_diagnostics_row public :: write_chemistry_diagnostics_header public :: write_cantera_cache_stats, write_chemistry_screening_stats public :: recover_temperature_and_update_thermo_cantera public :: compute_species_sensible_enthalpies_cantera !> Local accumulated chemistry-screening statistics. These are reduced !! across flow ranks at shutdown by write_chemistry_screening_stats. real(rk), save :: chemistry_stats_updates_local = zero real(rk), save :: chemistry_stats_owned_cells_local = zero real(rk), save :: chemistry_stats_active_cells_local = zero real(rk), save :: chemistry_stats_skipped_cells_local = zero real(rk), save :: chemistry_stats_reactor_calls_local = zero real(rk), save :: chemistry_stats_skipped_temperature_local = zero real(rk), save :: chemistry_stats_skipped_named_species_local = zero real(rk), save :: chemistry_stats_skipped_legacy_mass_fraction_local = zero real(rk), save :: chemistry_stats_all_skipped_updates_local = zero real(rk), save :: chemistry_stats_total_time_local = zero real(rk), save :: chemistry_stats_reactor_time_local = zero !> C-Binding interface for Cantera thermodynamics. !! !! The C++ bridge expects temperatures \(T\), thermodynamic pressures \(p_0\), !! and normalized solver species mass fractions \(Y_k\). It returns sensible !! enthalpy using !! $$ !! h_{sens}(T,Y,p_0) = !! h_{abs}(T,Y,p_0) - h_{abs}(T_{ref},Y,p_0) !! $$ !! and recovers temperature by adding the same-composition reference !! enthalpy before Cantera HP inversion. Projection pressure is never passed !! to these calls. interface subroutine cantera_get_cache_stats_c(stats_out, nstats) bind(c, name="cantera_get_cache_stats_c") import :: c_double, c_int real(c_double), intent(out) :: stats_out(*) integer(c_int), value :: nstats end subroutine cantera_get_cache_stats_c subroutine cantera_update_thermo_c(ncells, T, P, nspecies, Y_in, h_out, cp_out, lambda_out, & rho_thermo_out, T_ref, species_names_flat, name_len) bind(c, name="cantera_update_thermo_c") import :: c_char, c_double, c_int integer(c_int), value :: ncells real(c_double), intent(in) :: T(ncells), P(ncells) integer(c_int), value :: nspecies real(c_double), intent(in) :: Y_in(*) real(c_double), intent(out) :: h_out(ncells), cp_out(ncells), lambda_out(ncells), rho_thermo_out(ncells) real(c_double), value :: T_ref character(kind=c_char), intent(in) :: species_names_flat(*) integer(c_int), value :: name_len end subroutine cantera_update_thermo_c subroutine cantera_recover_temperature_from_h_c(ncells, h_in, P, nspecies, Y_in, T_out, & T_ref, species_names_flat, name_len) bind(c, name="cantera_recover_temperature_from_h_c") import :: c_char, c_double, c_int integer(c_int), value :: ncells real(c_double), intent(in) :: h_in(ncells), P(ncells) integer(c_int), value :: nspecies real(c_double), intent(in) :: Y_in(*) real(c_double), intent(out) :: T_out(ncells) real(c_double), value :: T_ref character(kind=c_char), intent(in) :: species_names_flat(*) integer(c_int), value :: name_len end subroutine cantera_recover_temperature_from_h_c subroutine cantera_species_sensible_enthalpies_c(npoints, T, P, nspecies, Y_in, hk_out, & T_ref, species_names_flat, name_len) & bind(c, name="cantera_species_sensible_enthalpies_c") import :: c_char, c_double, c_int integer(c_int), value :: npoints real(c_double), intent(in) :: T(npoints), P(npoints) integer(c_int), value :: nspecies real(c_double), intent(in) :: Y_in(*) real(c_double), intent(out) :: hk_out(*) real(c_double), value :: T_ref character(kind=c_char), intent(in) :: species_names_flat(*) integer(c_int), value :: name_len end subroutine cantera_species_sensible_enthalpies_c subroutine cantera_recover_temperature_and_update_thermo_c(ncells, h_in, P, nspecies, Y_in, T_out, & cp_out, lambda_out, rho_thermo_out, & T_ref, species_names_flat, name_len) & bind(c, name="cantera_recover_temperature_and_update_thermo_c") import :: c_char, c_double, c_int integer(c_int), value :: ncells real(c_double), intent(in) :: h_in(ncells), P(ncells) integer(c_int), value :: nspecies real(c_double), intent(in) :: Y_in(*) real(c_double), intent(out) :: T_out(ncells), cp_out(ncells), lambda_out(ncells), rho_thermo_out(ncells) real(c_double), value :: T_ref character(kind=c_char), intent(in) :: species_names_flat(*) integer(c_int), value :: name_len end subroutine cantera_recover_temperature_and_update_thermo_c subroutine cantera_recover_temperature_update_thermo_and_species_h_c(npoints, h_in, P, nspecies, Y_in, T_out, & cp_out, lambda_out, rho_thermo_out, hk_out, & T_ref, species_names_flat, name_len) & bind(c, name="cantera_recover_temperature_update_thermo_and_species_h_c") import :: c_char, c_double, c_int integer(c_int), value :: npoints real(c_double), intent(in) :: h_in(npoints), P(npoints) integer(c_int), value :: nspecies real(c_double), intent(in) :: Y_in(*) real(c_double), intent(out) :: T_out(npoints), cp_out(npoints), lambda_out(npoints), rho_thermo_out(npoints) real(c_double), intent(out) :: hk_out(*) real(c_double), value :: T_ref character(kind=c_char), intent(in) :: species_names_flat(*) integer(c_int), value :: name_len end subroutine cantera_recover_temperature_update_thermo_and_species_h_c subroutine cantera_heat_release_rate_c(npoints, T, P, nspecies, Y_in, qdot_out, & species_names_flat, name_len) & bind(c, name="cantera_heat_release_rate_c") import :: c_char, c_double, c_int integer(c_int), value :: npoints real(c_double), intent(in) :: T(npoints), P(npoints) integer(c_int), value :: nspecies real(c_double), intent(in) :: Y_in(*) real(c_double), intent(out) :: qdot_out(npoints) character(kind=c_char), intent(in) :: species_names_flat(*) integer(c_int), value :: name_len end subroutine cantera_heat_release_rate_c subroutine cantera_integrate_const_p_chemistry_c(npoints, dt, T_inout, P, nspecies, Y_inout, & h_out, cp_out, lambda_out, rho_out, T_ref, & species_names_flat, name_len, rtol, atol, & max_steps, energy_enabled) & bind(c, name="cantera_integrate_const_p_chemistry_c") import :: c_char, c_double, c_int integer(c_int), value :: npoints real(c_double), value :: dt real(c_double), intent(inout) :: T_inout(npoints) real(c_double), intent(in) :: P(npoints) integer(c_int), value :: nspecies real(c_double), intent(inout) :: Y_inout(*) real(c_double), intent(out) :: h_out(npoints), cp_out(npoints), lambda_out(npoints), rho_out(npoints) real(c_double), value :: T_ref character(kind=c_char), intent(in) :: species_names_flat(*) integer(c_int), value :: name_len real(c_double), value :: rtol, atol integer(c_int), value :: max_steps, energy_enabled end subroutine cantera_integrate_const_p_chemistry_c end interface contains !> Allocate all energy arrays for the mesh. subroutine allocate_energy(mesh, energy) type(mesh_t), intent(in) :: mesh type(energy_fields_t), intent(inout) :: energy call finalize_energy(energy) allocate(energy%T(mesh%ncells)) allocate(energy%T_old(mesh%ncells)) allocate(energy%h(mesh%ncells)) allocate(energy%h_old(mesh%ncells)) allocate(energy%rhs_old(mesh%ncells)) allocate(energy%qrad(mesh%ncells)) allocate(energy%cp(mesh%ncells)) allocate(energy%lambda(mesh%ncells)) allocate(energy%rho_thermo(mesh%ncells)) allocate(energy%operator_consistent_rho_h(mesh%ncells)) allocate(energy%species_enthalpy_diffusion(mesh%ncells)) allocate(energy%chem_heat_release_rate(mesh%ncells)) allocate(energy%chem_dT_dt(mesh%ncells)) allocate(energy%chem_dY_max_dt(mesh%ncells)) allocate(energy%chem_rel_rho_change_dt(mesh%ncells)) allocate(energy%chem_source_alpha(mesh%ncells)) allocate(energy%chem_active(mesh%ncells)) energy%T = zero energy%T_old = zero energy%h = zero energy%h_old = zero energy%rhs_old = zero energy%rhs_old_valid = .false. energy%qrad = zero energy%cp = zero energy%lambda = zero energy%rho_thermo = zero energy%operator_consistent_rho_h = zero energy%operator_consistent_rho_h_available = 0 energy%species_enthalpy_diffusion = zero energy%chem_heat_release_rate = zero energy%chem_dT_dt = zero energy%chem_dY_max_dt = zero energy%chem_rel_rho_change_dt = zero energy%chem_source_alpha = one energy%chem_active = zero energy%last_conductive_boundary_flux_out = zero energy%last_conductive_boundary_flux_available = 0 energy%cumulative_boundary_rho_h_advective_flux_out = zero energy%cumulative_boundary_rho_h_conductive_flux_out = zero energy%cumulative_qrad_integral = zero energy%cumulative_rho_species_hdiff_integral = zero energy%cumulative_energy_budget_available = 0 energy%last_energy_update_delta_rate_integral = zero energy%last_energy_update_rhs_integral = zero energy%last_energy_update_balance_defect = zero energy%relative_last_energy_update_balance_defect = zero energy%last_operator_consistent_rho_h_integral = zero energy%cumulative_energy_update_delta_integral = zero energy%cumulative_energy_update_rhs_integral = zero energy%initialized = .true. end subroutine allocate_energy !> Initialize energy fields from case parameters. subroutine initialize_energy(mesh, params, energy, species_Y) type(mesh_t), intent(in) :: mesh type(case_params_t), intent(in) :: params type(energy_fields_t), intent(inout) :: energy real(rk), intent(in), optional :: species_Y(:,:) call allocate_energy(mesh, energy) energy%T = params%initial_T energy%T_old = energy%T if (params%enable_cantera_thermo) then call update_thermo_from_temperature_cantera(mesh, params, energy, species_Y) else call update_enthalpy_from_temperature_constant_cp(params, energy) energy%cp = params%energy_cp energy%lambda = params%energy_lambda energy%rho_thermo = params%rho end if energy%h_old = energy%h call zero_radiation_source(energy) end subroutine initialize_energy !> Refresh h/cp/lambda/rho_thermo from the current temperature field. !! !! This is used by initialization utilities such as the ignition kernel that !! overwrite energy%T before the first timestep. It keeps the transported !! sensible enthalpy and derived thermodynamic fields consistent with the !! modified temperature and composition. subroutine refresh_energy_from_temperature(mesh, params, energy, species_Y) type(mesh_t), intent(in) :: mesh type(case_params_t), intent(in) :: params type(energy_fields_t), intent(inout) :: energy real(rk), intent(in), optional :: species_Y(:,:) if (.not. energy%initialized) then call fatal_error('energy', 'refresh requested before energy initialization') end if if (params%enable_cantera_thermo) then call update_thermo_from_temperature_cantera(mesh, params, energy, species_Y) else call update_enthalpy_from_temperature_constant_cp(params, energy) end if if (allocated(energy%T_old)) energy%T_old = energy%T if (allocated(energy%h_old)) energy%h_old = energy%h if (allocated(energy%rhs_old)) energy%rhs_old = zero energy%rhs_old_valid = .false. end subroutine refresh_energy_from_temperature !> Deallocate all energy arrays. subroutine finalize_energy(energy) type(energy_fields_t), intent(inout) :: energy if (allocated(energy%T)) deallocate(energy%T) if (allocated(energy%T_old)) deallocate(energy%T_old) if (allocated(energy%h)) deallocate(energy%h) if (allocated(energy%h_old)) deallocate(energy%h_old) if (allocated(energy%rhs_old)) deallocate(energy%rhs_old) if (allocated(energy%qrad)) deallocate(energy%qrad) if (allocated(energy%cp)) deallocate(energy%cp) if (allocated(energy%lambda)) deallocate(energy%lambda) if (allocated(energy%rho_thermo)) deallocate(energy%rho_thermo) if (allocated(energy%operator_consistent_rho_h)) deallocate(energy%operator_consistent_rho_h) if (allocated(energy%species_enthalpy_diffusion)) deallocate(energy%species_enthalpy_diffusion) if (allocated(energy%chem_heat_release_rate)) deallocate(energy%chem_heat_release_rate) if (allocated(energy%chem_dT_dt)) deallocate(energy%chem_dT_dt) if (allocated(energy%chem_dY_max_dt)) deallocate(energy%chem_dY_max_dt) if (allocated(energy%chem_rel_rho_change_dt)) deallocate(energy%chem_rel_rho_change_dt) if (allocated(energy%chem_source_alpha)) deallocate(energy%chem_source_alpha) if (allocated(energy%chem_active)) deallocate(energy%chem_active) energy%operator_consistent_rho_h_available = 0 energy%rhs_old_valid = .false. energy%initialized = .false. end subroutine finalize_energy !> Update h from T using the constant-cp thermodynamic model. subroutine update_enthalpy_from_temperature_constant_cp(params, energy) type(case_params_t), intent(in) :: params type(energy_fields_t), intent(inout) :: energy if (.not. allocated(energy%T) .or. .not. allocated(energy%h)) then call fatal_error('energy', 'energy arrays are not allocated') end if if (params%energy_cp <= tiny_safe) then call fatal_error('energy', 'energy_cp must be positive') end if energy%h = params%energy_reference_h + & params%energy_cp * (energy%T - params%energy_reference_T) if (allocated(energy%cp)) energy%cp = params%energy_cp if (allocated(energy%lambda)) energy%lambda = params%energy_lambda if (allocated(energy%rho_thermo)) energy%rho_thermo = params%rho end subroutine update_enthalpy_from_temperature_constant_cp !> Recover T from h using the constant-cp thermodynamic model. subroutine recover_temperature_constant_cp(params, energy) type(case_params_t), intent(in) :: params type(energy_fields_t), intent(inout) :: energy if (.not. allocated(energy%T) .or. .not. allocated(energy%h)) then call fatal_error('energy', 'energy arrays are not allocated') end if if (params%energy_cp <= tiny_safe) then call fatal_error('energy', 'energy_cp must be positive') end if energy%T = params%energy_reference_T + & (energy%h - params%energy_reference_h) / params%energy_cp end subroutine recover_temperature_constant_cp !> Constant-cp helper for a single boundary/face temperature value. pure function enthalpy_from_temperature_value(params, temperature) result(h_value) type(case_params_t), intent(in) :: params real(rk), intent(in) :: temperature real(rk) :: h_value h_value = params%energy_reference_h + & params%energy_cp * (temperature - params%energy_reference_T) end function enthalpy_from_temperature_value !> Build a flattened C-compatible species-name buffer. subroutine build_c_species_names(params, c_names_flat, n_len) type(case_params_t), intent(in) :: params character(kind=c_char), allocatable, intent(out) :: c_names_flat(:) integer, intent(out) :: n_len integer :: k, j, nsp nsp = max(1, params%nspecies) n_len = len(params%species_name(1)) allocate(c_names_flat(n_len * nsp)) c_names_flat = ' ' do k = 1, nsp do j = 1, n_len c_names_flat((k-1)*n_len + j) = char(iachar(params%species_name(k)(j:j)), kind=c_char) end do end do end subroutine build_c_species_names !> Resolve optional named chemistry active-species screening names to current species indices. !! !! This is intentionally done at chemistry-update time rather than during !! namelist validation because Cantera mechanism discovery may replace the !! user-specified species list after the namelist has been read. subroutine resolve_chemistry_active_species(params, nactive_species, active_species_indices) type(case_params_t), intent(in) :: params integer, intent(out) :: nactive_species integer, intent(out) :: active_species_indices(:) integer :: requested, i, k, found, nslots character(len=name_len) :: requested_name, mechanism_name logical :: duplicate nactive_species = 0 active_species_indices = 0 if (params%chemistry_active_species_threshold <= zero) return requested = params%chemistry_n_active_species if (requested <= 0) then do i = 1, max_species if (len_trim(params%chemistry_active_species_name(i)) > 0) requested = i end do end if if (requested <= 0) return nslots = min(size(active_species_indices), requested) do i = 1, nslots if (len_trim(params%chemistry_active_species_name(i)) <= 0) cycle requested_name = lowercase(trim(params%chemistry_active_species_name(i))) found = 0 do k = 1, params%nspecies mechanism_name = lowercase(trim(params%species_name(k))) if (trim(mechanism_name) == trim(requested_name)) then found = k exit end if end do if (found <= 0) then call fatal_error('chemistry', 'chemistry_active_species_name not found in active mechanism: '// & trim(params%chemistry_active_species_name(i))) end if duplicate = .false. do k = 1, nactive_species if (active_species_indices(k) == found) then duplicate = .true. exit end if end do if (.not. duplicate) then nactive_species = nactive_species + 1 active_species_indices(nactive_species) = found end if end do end subroutine resolve_chemistry_active_species !> Build a cellwise thermodynamic composition array for Cantera calls. !! !! If transported species are available, they are used directly and clipped !! to non-negative values. If their sum is positive, the local composition is !! normalized before passing to Cantera. If no transported species are !! available, use the single-species/default mixture prepared during Cantera !! initialization, usually thermo_default_species. subroutine build_thermo_Y(params, ncells, Y_local, species_Y) type(case_params_t), intent(in) :: params integer, intent(in) :: ncells real(rk), allocatable, intent(out) :: Y_local(:,:) real(rk), intent(in), optional :: species_Y(:,:) integer :: c, k, nsp real(rk) :: sum_Y nsp = max(1, params%nspecies) allocate(Y_local(nsp, ncells)) Y_local = zero if (present(species_Y) .and. params%enable_species .and. params%nspecies > 0) then if (size(species_Y, 1) < params%nspecies .or. size(species_Y, 2) < ncells) then call fatal_error('energy', 'species_Y has incompatible shape for Cantera thermo update') end if do c = 1, ncells sum_Y = zero do k = 1, params%nspecies Y_local(k, c) = max(zero, species_Y(k, c)) sum_Y = sum_Y + Y_local(k, c) end do if (sum_Y > tiny_safe) then Y_local(1:params%nspecies, c) = Y_local(1:params%nspecies, c) / sum_Y else Y_local(1, c) = 1.0_rk end if end do else Y_local(1, :) = 1.0_rk end if end subroutine build_thermo_Y !> Build a boundary thermodynamic composition vector for fixed-T boundaries. !! !! Fixed-composition species boundaries use `patch_Y` through !! `boundary_species`; non-Dirichlet species boundaries use the adjacent !! interior composition. Fixed-temperature enthalpy is then evaluated as !! \(h_b = h(T_b,Y_b,p_0)\). subroutine build_boundary_thermo_Y(mesh, bc, params, face_id, cell_id, species_Y, Y_point) type(mesh_t), intent(in) :: mesh type(bc_set_t), intent(in) :: bc type(case_params_t), intent(in) :: params integer, intent(in) :: face_id integer, intent(in) :: cell_id real(rk), intent(in) :: species_Y(:,:) real(rk), intent(out) :: Y_point(:) integer :: k real(rk) :: ext_Y, sum_Y logical :: is_dirichlet if (size(Y_point) < params%nspecies) then call fatal_error('energy', 'Y_point too small in build_boundary_thermo_Y') end if if (size(species_Y, 1) < params%nspecies .or. size(species_Y, 2) < cell_id) then call fatal_error('energy', 'species_Y has incompatible shape in boundary thermo composition') end if sum_Y = zero do k = 1, params%nspecies call boundary_species(mesh, bc, face_id, k, species_Y(k, cell_id), ext_Y, is_dirichlet) Y_point(k) = max(zero, ext_Y) sum_Y = sum_Y + Y_point(k) end do if (sum_Y > tiny_safe) then Y_point(1:params%nspecies) = Y_point(1:params%nspecies) / sum_Y else Y_point = zero Y_point(1) = 1.0_rk end if end subroutine build_boundary_thermo_Y !> Update h, cp, lambda, and diagnostic thermo density from current T. subroutine update_thermo_from_temperature_cantera(mesh, params, energy, species_Y) type(mesh_t), intent(in) :: mesh type(case_params_t), intent(in) :: params type(energy_fields_t), intent(inout) :: energy real(rk), intent(in), optional :: species_Y(:,:) integer :: n_len real(rk), allocatable :: P_arr(:), Y_local(:,:) character(kind=c_char), allocatable :: c_names_flat(:) if (.not. params%enable_cantera_thermo) return if (params%nspecies <= 0) then call fatal_error('energy', 'enable_cantera_thermo requires at least one inert thermo species') end if allocate(P_arr(mesh%ncells)) P_arr = params%background_press call build_thermo_Y(params, mesh%ncells, Y_local, species_Y) call build_c_species_names(params, c_names_flat, n_len) call cantera_update_thermo_c(mesh%ncells, energy%T, P_arr, params%nspecies, Y_local, & energy%h, energy%cp, energy%lambda, energy%rho_thermo, & params%energy_reference_T, c_names_flat, n_len) deallocate(P_arr) deallocate(Y_local) deallocate(c_names_flat) end subroutine update_thermo_from_temperature_cantera !> Recover T from h using Cantera HPY inversion. subroutine recover_temperature_cantera(mesh, params, energy, species_Y) type(mesh_t), intent(in) :: mesh type(case_params_t), intent(in) :: params type(energy_fields_t), intent(inout) :: energy real(rk), intent(in), optional :: species_Y(:,:) integer :: n_len real(rk), allocatable :: P_arr(:), Y_local(:,:) character(kind=c_char), allocatable :: c_names_flat(:) if (.not. params%enable_cantera_thermo) return if (params%nspecies <= 0) then call fatal_error('energy', 'enable_cantera_thermo requires at least one inert thermo species') end if allocate(P_arr(mesh%ncells)) P_arr = params%background_press call build_thermo_Y(params, mesh%ncells, Y_local, species_Y) call build_c_species_names(params, c_names_flat, n_len) call cantera_recover_temperature_from_h_c(mesh%ncells, & energy%h, & P_arr, & params%nspecies, & Y_local, & energy%T, & params%energy_reference_T, & c_names_flat, & n_len) deallocate(P_arr) deallocate(Y_local) deallocate(c_names_flat) end subroutine recover_temperature_cantera !> Recover T from h and refresh cp/lambda/rho_thermo in one Cantera sync. !! !! This is the preferred energy-step thermo sync routine: !! `(T, cp, lambda, rho_thermo) = sync(h, Y, p0)`. It preserves the !! transported sensible enthalpy field and updates only derived thermo !! state: T, cp, lambda, and rho_thermo. The thermo-sync cache key in the !! C++ bridge is `h`, `p0`, and `Y`; `rho_thermo` becomes active density only !! after `sync_active_density_from_thermo` in variable-density mode. subroutine recover_temperature_and_update_thermo_cantera(mesh, params, energy, species_Y) type(mesh_t), intent(in) :: mesh type(case_params_t), intent(in) :: params type(energy_fields_t), intent(inout) :: energy real(rk), intent(in), optional :: species_Y(:,:) integer :: n_len real(rk), allocatable :: P_arr(:), Y_local(:,:) character(kind=c_char), allocatable :: c_names_flat(:) if (.not. params%enable_cantera_thermo) return if (params%nspecies <= 0) then call fatal_error('energy', 'enable_cantera_thermo requires at least one inert thermo species') end if allocate(P_arr(mesh%ncells)) P_arr = params%background_press call build_thermo_Y(params, mesh%ncells, Y_local, species_Y) call build_c_species_names(params, c_names_flat, n_len) call cantera_recover_temperature_and_update_thermo_c(mesh%ncells, & energy%h, & P_arr, & params%nspecies, & Y_local, & energy%T, & energy%cp, & energy%lambda, & energy%rho_thermo, & params%energy_reference_T, & c_names_flat, & n_len) deallocate(P_arr) deallocate(Y_local) deallocate(c_names_flat) end subroutine recover_temperature_and_update_thermo_cantera !> Owned-cell thermo sync for the timestep hot path. !! !! The legacy helper above updates the replicated full mesh. During a !! timestep each MPI rank only needs to evaluate Cantera on its owned cells; !! the caller performs the normal halo exchanges for stencil fields. This !! preserves the separate radiation MPI communicator because all exchanges !! still go through flow_mpi_t. subroutine recover_temperature_and_update_thermo_cantera_owned(mesh, flow, params, energy, species_Y) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(case_params_t), intent(in) :: params type(energy_fields_t), intent(inout) :: energy real(rk), intent(in), optional :: species_Y(:,:) integer :: n_len, nloc, i, c, k real(rk) :: sum_Y real(rk), allocatable :: h_local(:), P_arr(:), Y_local(:,:) real(rk), allocatable :: T_local(:), cp_local(:), lambda_local(:), rho_local(:) character(kind=c_char), allocatable :: c_names_flat(:) if (.not. params%enable_cantera_thermo) return if (params%nspecies <= 0) then call fatal_error('energy', 'enable_cantera_thermo requires at least one inert thermo species') end if if (mesh%ncells < flow%last_cell) then call fatal_error('energy', 'flow ownership exceeds mesh size in owned Cantera thermo sync') end if nloc = max(0, flow%last_cell - flow%first_cell + 1) if (nloc <= 0) return allocate(h_local(nloc), P_arr(nloc), Y_local(max(1, params%nspecies), nloc)) allocate(T_local(nloc), cp_local(nloc), lambda_local(nloc), rho_local(nloc)) P_arr = params%background_press do i = 1, nloc c = flow%first_cell + i - 1 h_local(i) = energy%h(c) if (present(species_Y) .and. params%nspecies > 0) then sum_Y = zero do k = 1, params%nspecies Y_local(k, i) = max(zero, species_Y(k, c)) sum_Y = sum_Y + Y_local(k, i) end do if (sum_Y > tiny_safe) then Y_local(1:params%nspecies, i) = Y_local(1:params%nspecies, i) / sum_Y else Y_local(:, i) = zero Y_local(1, i) = 1.0_rk end if else Y_local(:, i) = zero Y_local(1, i) = 1.0_rk end if end do call build_c_species_names(params, c_names_flat, n_len) call cantera_recover_temperature_and_update_thermo_c(nloc, h_local, P_arr, params%nspecies, & Y_local, T_local, cp_local, lambda_local, rho_local, & params%energy_reference_T, c_names_flat, n_len) do i = 1, nloc c = flow%first_cell + i - 1 energy%T(c) = T_local(i) energy%cp(c) = cp_local(i) energy%lambda(c) = lambda_local(i) energy%rho_thermo(c) = rho_local(i) end do deallocate(h_local, P_arr, Y_local, T_local, cp_local, lambda_local, rho_local, c_names_flat) end subroutine recover_temperature_and_update_thermo_cantera_owned !> Owned-cell fused thermo sync plus species sensible enthalpies. subroutine recover_temperature_update_thermo_and_species_h_cantera_owned(mesh, flow, params, energy, species_Y, hk_species) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(inout) :: flow type(case_params_t), intent(in) :: params type(energy_fields_t), intent(inout) :: energy real(rk), intent(in) :: species_Y(:,:) real(rk), intent(out) :: hk_species(:,:) integer :: n_len, nloc, i, c, k real(rk) :: sum_Y real(rk), allocatable :: h_local(:), P_arr(:), Y_local(:,:) real(rk), allocatable :: T_local(:), cp_local(:), lambda_local(:), rho_local(:), hk_local(:,:) character(kind=c_char), allocatable :: c_names_flat(:) if (.not. params%enable_cantera_thermo) then call fatal_error('energy', 'fused thermo/species enthalpy sync requires Cantera thermo') end if if (params%nspecies <= 0) then call fatal_error('energy', 'fused thermo/species enthalpy sync requires at least one species') end if if (size(hk_species, 1) < params%nspecies .or. size(hk_species, 2) < mesh%ncells) then call fatal_error('energy', 'hk_species has incompatible shape in fused thermo/species enthalpy sync') end if nloc = max(0, flow%last_cell - flow%first_cell + 1) if (nloc <= 0) return allocate(h_local(nloc), P_arr(nloc), Y_local(max(1, params%nspecies), nloc)) allocate(T_local(nloc), cp_local(nloc), lambda_local(nloc), rho_local(nloc)) allocate(hk_local(params%nspecies, nloc)) P_arr = params%background_press do i = 1, nloc c = flow%first_cell + i - 1 h_local(i) = energy%h(c) sum_Y = zero do k = 1, params%nspecies Y_local(k, i) = max(zero, species_Y(k, c)) sum_Y = sum_Y + Y_local(k, i) end do if (sum_Y > tiny_safe) then Y_local(1:params%nspecies, i) = Y_local(1:params%nspecies, i) / sum_Y else Y_local(:, i) = zero Y_local(1, i) = 1.0_rk end if end do call build_c_species_names(params, c_names_flat, n_len) call cantera_recover_temperature_update_thermo_and_species_h_c(nloc, h_local, P_arr, params%nspecies, & Y_local, T_local, cp_local, lambda_local, rho_local, hk_local, & params%energy_reference_T, c_names_flat, n_len) do i = 1, nloc c = flow%first_cell + i - 1 energy%T(c) = T_local(i) energy%cp(c) = cp_local(i) energy%lambda(c) = lambda_local(i) energy%rho_thermo(c) = rho_local(i) do k = 1, params%nspecies hk_species(k, c) = hk_local(k, i) end do end do call flow_exchange_cell_matrix(flow, hk_species) deallocate(h_local, P_arr, Y_local, T_local, cp_local, lambda_local, rho_local, hk_local, c_names_flat) end subroutine recover_temperature_update_thermo_and_species_h_cantera_owned !> Compute species sensible enthalpies h_k(T)-h_k(T_ref) [J/kg_k] for all cells. subroutine compute_species_sensible_enthalpies_cantera(mesh, params, T_state, species_Y, hk_species) type(mesh_t), intent(in) :: mesh type(case_params_t), intent(in) :: params real(rk), intent(in) :: T_state(:) real(rk), intent(in) :: species_Y(:,:) real(rk), intent(out) :: hk_species(:,:) integer :: n_len real(rk), allocatable :: P_arr(:), Y_local(:,:) character(kind=c_char), allocatable :: c_names_flat(:) if (.not. params%enable_cantera_thermo) then call fatal_error('energy', 'species sensible enthalpies require Cantera thermo') end if if (params%nspecies <= 0) then call fatal_error('energy', 'species sensible enthalpies require at least one species') end if if (size(T_state) < mesh%ncells) then call fatal_error('energy', 'T_state has incompatible size for species sensible enthalpies') end if if (size(hk_species, 1) < params%nspecies .or. size(hk_species, 2) < mesh%ncells) then call fatal_error('energy', 'hk_species has incompatible shape') end if allocate(P_arr(mesh%ncells)) P_arr = params%background_press call build_thermo_Y(params, mesh%ncells, Y_local, species_Y) call build_c_species_names(params, c_names_flat, n_len) call cantera_species_sensible_enthalpies_c(mesh%ncells, T_state, P_arr, params%nspecies, & Y_local, hk_species, params%energy_reference_T, & c_names_flat, n_len) deallocate(P_arr) deallocate(Y_local) deallocate(c_names_flat) end subroutine compute_species_sensible_enthalpies_cantera !> Compute species sensible enthalpies for one boundary/face state. subroutine species_sensible_enthalpies_from_temperature_cantera_point(params, temperature, hk_value, Y_point) type(case_params_t), intent(in) :: params real(rk), intent(in) :: temperature real(rk), intent(out) :: hk_value(:) real(rk), intent(in) :: Y_point(:) integer :: n_len, nsp real(rk) :: T_arr(1), P_arr(1) real(rk), allocatable :: Y_local(:,:) character(kind=c_char), allocatable :: c_names_flat(:) if (params%nspecies <= 0) then call fatal_error('energy', 'Cantera boundary species enthalpies require at least one species') end if if (size(hk_value) < params%nspecies) then call fatal_error('energy', 'species enthalpy point hk array has incompatible size') end if if (size(Y_point) < params%nspecies) then call fatal_error('energy', 'species enthalpy point Y array has incompatible size') end if nsp = max(1, params%nspecies) allocate(Y_local(nsp, 1)) Y_local = zero Y_local(1:params%nspecies, 1) = max(zero, Y_point(1:params%nspecies)) if (sum(Y_local(1:params%nspecies, 1)) > tiny_safe) then Y_local(1:params%nspecies, 1) = Y_local(1:params%nspecies, 1) / sum(Y_local(1:params%nspecies, 1)) else Y_local(1, 1) = 1.0_rk end if T_arr(1) = temperature P_arr(1) = params%background_press call build_c_species_names(params, c_names_flat, n_len) call cantera_species_sensible_enthalpies_c(1, T_arr, P_arr, params%nspecies, Y_local, & hk_value, params%energy_reference_T, & c_names_flat, n_len) deallocate(Y_local) deallocate(c_names_flat) end subroutine species_sensible_enthalpies_from_temperature_cantera_point !> Compute h(T,Y,p0) for one boundary state using Cantera. subroutine enthalpy_from_temperature_cantera_point(params, temperature, h_value, Y_point) type(case_params_t), intent(in) :: params real(rk), intent(in) :: temperature real(rk), intent(out) :: h_value real(rk), intent(in), optional :: Y_point(:) integer :: n_len, nsp real(rk) :: T_arr(1), P_arr(1) real(rk) :: h_arr(1), cp_arr(1), lambda_arr(1), rho_arr(1) real(rk), allocatable :: Y_local(:,:) character(kind=c_char), allocatable :: c_names_flat(:) if (params%nspecies <= 0) then call fatal_error('energy', 'Cantera boundary enthalpy requires at least one species') end if nsp = max(1, params%nspecies) allocate(Y_local(nsp, 1)) Y_local = zero if (present(Y_point)) then if (size(Y_point) < params%nspecies) then call fatal_error('energy', 'Y_point has incompatible size for Cantera boundary enthalpy') end if Y_local(1:params%nspecies, 1) = max(zero, Y_point(1:params%nspecies)) if (sum(Y_local(1:params%nspecies, 1)) > tiny_safe) then Y_local(1:params%nspecies, 1) = Y_local(1:params%nspecies, 1) / & sum(Y_local(1:params%nspecies, 1)) else Y_local(1, 1) = 1.0_rk end if else Y_local(1, 1) = 1.0_rk end if T_arr(1) = temperature P_arr(1) = params%background_press call build_c_species_names(params, c_names_flat, n_len) call cantera_update_thermo_c(1, T_arr, P_arr, params%nspecies, Y_local, & h_arr, cp_arr, lambda_arr, rho_arr, & params%energy_reference_T, c_names_flat, n_len) h_value = h_arr(1) deallocate(Y_local) deallocate(c_names_flat) end subroutine enthalpy_from_temperature_cantera_point !> Convert a boundary temperature to enthalpy using the active thermo model. !> Boundary enthalpy from a boundary temperature. !! !! When a transported boundary composition is available, the fixed-temperature !! boundary enthalpy is evaluated as h(T_b,Y_b,p0). This keeps hot/cold !! species inlets thermodynamically consistent. function boundary_enthalpy_from_temperature(mesh, params, temperature, Y_point) result(h_value) type(mesh_t), intent(in) :: mesh type(case_params_t), intent(in) :: params real(rk), intent(in) :: temperature real(rk), intent(in), optional :: Y_point(:) real(rk) :: h_value if (.not. params%enable_cantera_thermo) then h_value = enthalpy_from_temperature_value(params, temperature) return end if call enthalpy_from_temperature_cantera_point(params, temperature, h_value, Y_point) associate(dummy => mesh%ncells) end associate end function boundary_enthalpy_from_temperature !> Reset the volumetric energy source to zero. subroutine zero_radiation_source(energy) type(energy_fields_t), intent(inout) :: energy if (allocated(energy%qrad)) energy%qrad = zero end subroutine zero_radiation_source !> Initialize chemistry workload diagnostics CSV. subroutine write_chemistry_diagnostics_header(params, flow) type(case_params_t), intent(in) :: params type(flow_mpi_t), intent(in) :: flow integer :: unit_id, ios character(len=1024) :: filename if (.not. params%enable_reactions) return if (.not. params%write_diagnostics) return if (flow%rank /= 0) return call execute_command_line('mkdir -p ' // trim(params%output_dir) // '/diagnostics') filename = trim(params%output_dir) // '/diagnostics/chemistry_diagnostics.csv' 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 create chemistry diagnostics file: ', trim(filename) return end if write(unit_id,'(a)') 'step,chemistry_dt,nprocs,owned_cells,active_cells,skipped_cells,reactor_calls,'// & 'active_cells_min_rank,active_cells_max_rank,skipped_cells_min_rank,skipped_cells_max_rank,'// & 'total_time_min,total_time_max,total_time_avg,reactor_time_min,reactor_time_max,'// & 'reactor_time_avg,reactor_time_per_call_avg,max_dT_chem,max_dY_chem,'// & 'max_rel_drho_chem,min_source_alpha,raw_max_dT_chem,raw_max_dY_chem,raw_max_rel_drho_chem,'// & 'max_qdot_exact,integrated_qdot_exact,energy_enabled' close(unit_id) end subroutine write_chemistry_diagnostics_header subroutine write_chemistry_diagnostics_row(flow, params, step, chemistry_dt, owned_local, active_local, & skipped_local, total_time_local, reactor_time_local, & max_dT_local, max_dY_local, max_rel_drho_local, min_source_alpha_local, & raw_max_dT_local, raw_max_dY_local, raw_max_rel_drho_local, & max_qdot_local, integrated_qdot_local) type(flow_mpi_t), intent(in) :: flow type(case_params_t), intent(in) :: params integer, intent(in) :: step real(rk), intent(in) :: chemistry_dt integer, intent(in) :: owned_local, active_local, skipped_local real(rk), intent(in) :: total_time_local, reactor_time_local real(rk), intent(in), optional :: max_dT_local, max_dY_local, max_rel_drho_local, min_source_alpha_local real(rk), intent(in), optional :: raw_max_dT_local, raw_max_dY_local, raw_max_rel_drho_local real(rk), intent(in), optional :: max_qdot_local, integrated_qdot_local integer :: ierr, unit_id, ios integer :: local_counts(4), sum_counts(4) integer :: local_min_counts(2), min_counts(2), max_counts(2) real(rk) :: local_times(2), sum_times(2), min_times(2), max_times(2) real(rk) :: local_chem_metrics(3), max_chem_metrics(3) real(rk) :: local_raw_chem_metrics(3), max_raw_chem_metrics(3) real(rk) :: local_qdot_metrics(2), sum_qdot_metrics(2), max_qdot_metric real(rk) :: local_alpha, min_alpha real(rk) :: avg_total_time, avg_reactor_time, avg_reactor_time_per_call character(len=1024) :: filename if (.not. params%write_diagnostics) return if (mod(step, params%output_interval) /= 0 .and. step /= params%nsteps .and. step /= 0) return local_counts = [owned_local, active_local, skipped_local, active_local] local_min_counts = [active_local, skipped_local] local_times = [total_time_local, reactor_time_local] local_chem_metrics = zero local_raw_chem_metrics = zero local_qdot_metrics = zero if (present(max_dT_local)) local_chem_metrics(1) = max(zero, max_dT_local) if (present(max_dY_local)) local_chem_metrics(2) = max(zero, max_dY_local) if (present(max_rel_drho_local)) local_chem_metrics(3) = max(zero, max_rel_drho_local) if (present(raw_max_dT_local)) local_raw_chem_metrics(1) = max(zero, raw_max_dT_local) if (present(raw_max_dY_local)) local_raw_chem_metrics(2) = max(zero, raw_max_dY_local) if (present(raw_max_rel_drho_local)) local_raw_chem_metrics(3) = max(zero, raw_max_rel_drho_local) if (present(max_qdot_local)) local_qdot_metrics(1) = max(zero, max_qdot_local) if (present(integrated_qdot_local)) local_qdot_metrics(2) = integrated_qdot_local local_alpha = 1.0_rk if (present(min_source_alpha_local)) local_alpha = min(1.0_rk, max(zero, min_source_alpha_local)) call MPI_Reduce(local_counts, sum_counts, 4, MPI_INTEGER, MPI_SUM, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Reduce failed for chemistry count sums') call MPI_Reduce(local_min_counts, min_counts, 2, MPI_INTEGER, MPI_MIN, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Reduce failed for chemistry count minima') call MPI_Reduce(local_min_counts, max_counts, 2, MPI_INTEGER, MPI_MAX, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Reduce failed for chemistry count maxima') call MPI_Reduce(local_times, sum_times, 2, MPI_DOUBLE_PRECISION, MPI_SUM, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Reduce failed for chemistry time sums') call MPI_Reduce(local_times, min_times, 2, MPI_DOUBLE_PRECISION, MPI_MIN, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Reduce failed for chemistry time minima') call MPI_Reduce(local_times, max_times, 2, MPI_DOUBLE_PRECISION, MPI_MAX, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Reduce failed for chemistry time maxima') call MPI_Reduce(local_chem_metrics, max_chem_metrics, 3, MPI_DOUBLE_PRECISION, MPI_MAX, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Reduce failed for chemistry source maxima') call MPI_Reduce(local_raw_chem_metrics, max_raw_chem_metrics, 3, MPI_DOUBLE_PRECISION, MPI_MAX, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Reduce failed for raw chemistry source maxima') call MPI_Reduce(local_alpha, min_alpha, 1, MPI_DOUBLE_PRECISION, MPI_MIN, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Reduce failed for chemistry source alpha minimum') call MPI_Reduce(local_qdot_metrics, sum_qdot_metrics, 2, MPI_DOUBLE_PRECISION, MPI_SUM, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Reduce failed for chemistry heat-release sums') call MPI_Reduce(local_qdot_metrics(1), max_qdot_metric, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Reduce failed for chemistry heat-release maximum') if (flow%rank /= 0) return avg_total_time = sum_times(1) / real(max(1, flow%nprocs), rk) avg_reactor_time = sum_times(2) / real(max(1, flow%nprocs), rk) if (sum_counts(4) > 0) then avg_reactor_time_per_call = sum_times(2) / real(sum_counts(4), rk) else avg_reactor_time_per_call = zero end if filename = trim(params%output_dir) // '/diagnostics/chemistry_diagnostics.csv' open(newunit=unit_id, file=trim(filename), status='old', position='append', action='write', iostat=ios) if (ios /= 0) return write(unit_id,'(i0,a,es16.8,a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,a,'// & 'es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,'// & 'es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,'// & 'es16.8,a,es16.8,a,l1)') & step, ',', chemistry_dt, ',', flow%nprocs, ',', sum_counts(1), ',', sum_counts(2), ',', & sum_counts(3), ',', sum_counts(4), ',', min_counts(1), ',', max_counts(1), ',', & min_counts(2), ',', max_counts(2), ',', min_times(1), ',', max_times(1), ',', & avg_total_time, ',', min_times(2), ',', max_times(2), ',', avg_reactor_time, ',', & avg_reactor_time_per_call, ',', max_chem_metrics(1), ',', max_chem_metrics(2), ',', & max_chem_metrics(3), ',', min_alpha, ',', max_raw_chem_metrics(1), ',', & max_raw_chem_metrics(2), ',', max_raw_chem_metrics(3), ',', max_qdot_metric, ',', & sum_qdot_metrics(2), ',', params%chemistry_energy_enabled close(unit_id) end subroutine write_chemistry_diagnostics_row !> Accumulate local chemistry workload/screening statistics. subroutine accumulate_chemistry_screening_stats(nowned, nactive, skipped_temperature, & skipped_named_species, skipped_legacy_mass_fraction, & total_time, reactor_time, reactor_calls) integer, intent(in) :: nowned, nactive integer, intent(in) :: skipped_temperature, skipped_named_species, skipped_legacy_mass_fraction real(rk), intent(in) :: total_time, reactor_time integer, intent(in), optional :: reactor_calls integer :: nskipped nskipped = max(0, nowned - nactive) chemistry_stats_updates_local = chemistry_stats_updates_local + 1.0_rk chemistry_stats_owned_cells_local = chemistry_stats_owned_cells_local + real(nowned, rk) chemistry_stats_active_cells_local = chemistry_stats_active_cells_local + real(nactive, rk) chemistry_stats_skipped_cells_local = chemistry_stats_skipped_cells_local + real(nskipped, rk) if (present(reactor_calls)) then chemistry_stats_reactor_calls_local = chemistry_stats_reactor_calls_local + real(reactor_calls, rk) else chemistry_stats_reactor_calls_local = chemistry_stats_reactor_calls_local + real(nactive, rk) end if chemistry_stats_skipped_temperature_local = chemistry_stats_skipped_temperature_local + real(skipped_temperature, rk) chemistry_stats_skipped_named_species_local = chemistry_stats_skipped_named_species_local + real(skipped_named_species, rk) chemistry_stats_skipped_legacy_mass_fraction_local = chemistry_stats_skipped_legacy_mass_fraction_local + & real(skipped_legacy_mass_fraction, rk) if (nowned > 0 .and. nactive <= 0) then chemistry_stats_all_skipped_updates_local = chemistry_stats_all_skipped_updates_local + 1.0_rk end if chemistry_stats_total_time_local = chemistry_stats_total_time_local + total_time chemistry_stats_reactor_time_local = chemistry_stats_reactor_time_local + reactor_time end subroutine accumulate_chemistry_screening_stats !> Print accumulated chemistry workload and screening effectiveness statistics. !! !! These diagnostics are intentionally terminal-only, like the Cantera cache !! statistics. The per-step CSV remains the detailed machine-readable record. subroutine write_chemistry_screening_stats(flow, params) type(flow_mpi_t), intent(in) :: flow type(case_params_t), intent(in) :: params integer, parameter :: nstats = 11 real(rk) :: local_stats(nstats), global_stats(nstats) real(rk) :: owned_total, active_total, skipped_total, reactor_calls_total real(rk) :: skip_fraction, active_fraction, avg_reactor_time_per_call real(rk) :: estimated_reactor_time_avoided, update_steps logical :: screening_enabled integer :: ierr if (.not. params%enable_reactions) return local_stats = [chemistry_stats_updates_local, & chemistry_stats_owned_cells_local, & chemistry_stats_active_cells_local, & chemistry_stats_skipped_cells_local, & chemistry_stats_reactor_calls_local, & chemistry_stats_skipped_temperature_local, & chemistry_stats_skipped_named_species_local, & chemistry_stats_skipped_legacy_mass_fraction_local, & chemistry_stats_all_skipped_updates_local, & chemistry_stats_total_time_local, & chemistry_stats_reactor_time_local] global_stats = zero call MPI_Reduce(local_stats, global_stats, nstats, MPI_DOUBLE_PRECISION, MPI_SUM, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Reduce failed for chemistry screening statistics') if (flow%rank /= 0) return if (global_stats(1) <= zero) return owned_total = global_stats(2) active_total = global_stats(3) skipped_total = global_stats(4) reactor_calls_total = global_stats(5) update_steps = global_stats(1) / real(max(1, flow%nprocs), rk) if (owned_total > zero) then active_fraction = 100.0_rk * active_total / owned_total skip_fraction = 100.0_rk * skipped_total / owned_total else active_fraction = zero skip_fraction = zero end if if (reactor_calls_total > zero) then avg_reactor_time_per_call = global_stats(11) / reactor_calls_total else avg_reactor_time_per_call = zero end if estimated_reactor_time_avoided = skipped_total * avg_reactor_time_per_call screening_enabled = (params%chemistry_temperature_cutoff > zero) .or. & (params%chemistry_min_reactive_mass_fraction > zero) .or. & (params%chemistry_active_species_threshold > zero .and. & params%chemistry_n_active_species > 0) write(output_unit,'(a)') ' ======================================================================' write(output_unit,'(a)') ' CHEMISTRY SCREENING STATISTICS' write(output_unit,'(a)') ' Counts are summed over flow ranks and chemistry-update calls.' write(output_unit,'(a)') ' ======================================================================' write(output_unit,'(a,l1)') ' screening controls enabled : ', screening_enabled write(output_unit,'(a,f12.0)') ' chemistry update steps : ', update_steps write(output_unit,'(a,f12.0)') ' rank-update samples : ', global_stats(1) write(output_unit,'(a,f14.0)') ' owned cell-visits : ', owned_total write(output_unit,'(a,f14.0)') ' active reactor calls : ', reactor_calls_total write(output_unit,'(a,f14.0)') ' skipped cell-visits : ', skipped_total write(output_unit,'(a,f10.2)') ' active fraction [%] : ', active_fraction write(output_unit,'(a,f10.2)') ' skipped fraction [%] : ', skip_fraction write(output_unit,'(a,f14.0)') ' skipped by temperature cutoff : ', global_stats(6) write(output_unit,'(a,f14.0)') ' skipped by active-species screen : ', global_stats(7) write(output_unit,'(a,f14.0)') ' skipped by legacy max-Y screen : ', global_stats(8) write(output_unit,'(a,f14.0)') ' updates with all owned cells skip : ', global_stats(9) write(output_unit,'(a,es12.5)') ' avg reactor time per call [s] : ', avg_reactor_time_per_call write(output_unit,'(a,es12.5)') ' estimated ReactorNet time avoided : ', estimated_reactor_time_avoided write(output_unit,'(a,es12.5)') ' accumulated chemistry wall time : ', global_stats(10) write(output_unit,'(a,es12.5)') ' accumulated ReactorNet wall time : ', global_stats(11) write(output_unit,'(a)') ' ======================================================================' end subroutine write_chemistry_screening_stats !> Advance cell-local Cantera chemistry on owned cells. !! !! This is an operator-split chemistry hook. Each owned cell is integrated !! as an adiabatic constant-pressure reactor for the accumulated chemistry !! interval using Cantera's !! ReactorNet path. The routine updates transported species mass fractions !! and the energy dependent state (`T`, sensible `h`, `cp`, `lambda`, !! `rho_thermo`) for owned cells, then exchanges the updated replicated/ghost !! values using the flow communicator. Radiation MPI remains untouched. subroutine advance_cantera_chemistry(mesh, flow, params, species_Y, energy, step, chemistry_dt, & max_dT_chem, max_dY_chem, max_rel_drho_chem, min_source_alpha, & raw_max_dT_chem, raw_max_dY_chem, raw_max_rel_drho_chem) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(inout) :: flow type(case_params_t), intent(in) :: params real(rk), intent(inout) :: species_Y(:,:) type(energy_fields_t), intent(inout) :: energy integer, intent(in) :: step real(rk), intent(in), optional :: chemistry_dt real(rk), intent(out), optional :: max_dT_chem, max_dY_chem, max_rel_drho_chem, min_source_alpha real(rk), intent(out), optional :: raw_max_dT_chem, raw_max_dY_chem, raw_max_rel_drho_chem integer :: nloc, nactive, global_active, i, c, k, n_len, ierr, j, n_sub integer :: nactive_screen_species, active_screen_species(max_species) integer :: skipped_temperature_local, skipped_named_species_local, skipped_legacy_mass_fraction_local integer, allocatable :: active_cell(:) real(rk) :: sum_Y, dt_local, reactive_indicator, dt_sub real(rk) :: dT_abs, dY_abs, rel_drho, alpha, alpha_base real(rk) :: max_dT_local, max_dY_local, max_rel_drho_local, min_alpha_local real(rk) :: raw_max_dT_local, raw_max_dY_local, raw_max_rel_drho_local real(rk) :: qdot_max_local, qdot_integral_local real(rk) :: reduce_in(3), reduce_out(3), raw_reduce_in(3), raw_reduce_out(3), alpha_reduce logical :: use_named_species_screen, use_source_limiter real(rk) :: t_total_start, t_reactor_start, total_time_local, reactor_time_local real(rk), allocatable :: T_local(:), P_local(:), Y_local(:,:) real(rk), allocatable :: T_before(:), rho_before(:), Y_before(:,:) real(rk), allocatable :: h_local(:), cp_local(:), lambda_local(:), rho_local(:), alpha_local(:), qdot_local(:) character(kind=c_char), allocatable :: c_names_flat(:) real(rk), allocatable :: prev_chem_dT_dt(:) real(rk) :: dT_max_est ! Load-balancing packed dynamic buffers integer :: nfields, r real(rk) :: qdot_val, T_prev, rho_prev, dY_max_val integer, allocatable :: pack_cell(:), all_active_cells(:) real(rk), allocatable :: pack_data(:,:), all_data_new(:,:) integer, allocatable :: active_counts(:), active_displs(:) integer, allocatable :: active_matrix_counts(:), active_matrix_displs(:) n_sub = 1 if (present(max_dT_chem)) max_dT_chem = zero if (present(max_dY_chem)) max_dY_chem = zero if (present(max_rel_drho_chem)) max_rel_drho_chem = zero if (present(min_source_alpha)) min_source_alpha = 1.0_rk if (present(raw_max_dT_chem)) raw_max_dT_chem = zero if (present(raw_max_dY_chem)) raw_max_dY_chem = zero if (present(raw_max_rel_drho_chem)) raw_max_rel_drho_chem = zero if (.not. params%enable_reactions) return if (params%chemistry_update_interval <= 0) return t_total_start = real(MPI_Wtime(), rk) reactor_time_local = zero total_time_local = zero if (.not. params%enable_cantera_thermo) then call fatal_error('chemistry', 'Cantera chemistry requires enable_cantera_thermo=.true.') end if if (.not. energy%initialized) then call fatal_error('chemistry', 'energy must be initialized before chemistry') end if if (params%nspecies <= 0) then call fatal_error('chemistry', 'Cantera chemistry requires nspecies > 0 after mechanism discovery') end if if (size(species_Y,1) < params%nspecies .or. size(species_Y,2) < mesh%ncells) then call fatal_error('chemistry', 'species_Y has incompatible shape') end if call resolve_chemistry_active_species(params, nactive_screen_species, active_screen_species) use_named_species_screen = params%chemistry_active_species_threshold > zero .and. & nactive_screen_species > 0 if (present(chemistry_dt)) then dt_local = chemistry_dt else if (step <= 1) then dt_local = params%dt else dt_local = params%dt * real(params%chemistry_update_interval, rk) end if if (dt_local <= zero) return if (allocated(energy%chem_dT_dt)) then allocate(prev_chem_dT_dt(size(energy%chem_dT_dt))) prev_chem_dT_dt = energy%chem_dT_dt end if if (allocated(energy%chem_heat_release_rate)) energy%chem_heat_release_rate = zero if (allocated(energy%chem_dT_dt)) energy%chem_dT_dt = zero if (allocated(energy%chem_dY_max_dt)) energy%chem_dY_max_dt = zero if (allocated(energy%chem_rel_rho_change_dt)) energy%chem_rel_rho_change_dt = zero if (allocated(energy%chem_source_alpha)) energy%chem_source_alpha = one if (allocated(energy%chem_active)) energy%chem_active = zero if (params%enable_chemistry_load_balancing) then nloc = mesh%ncells allocate(active_cell(nloc)) nactive = 0 skipped_temperature_local = 0 skipped_named_species_local = 0 skipped_legacy_mass_fraction_local = 0 global_active = 0 do c = 1, mesh%ncells if (params%chemistry_temperature_cutoff > zero) then if (energy%T(c) < params%chemistry_temperature_cutoff) then skipped_temperature_local = skipped_temperature_local + 1 cycle end if end if if (use_named_species_screen) then reactive_indicator = zero do k = 1, nactive_screen_species reactive_indicator = max(reactive_indicator, & max(zero, species_Y(active_screen_species(k), c))) end do if (reactive_indicator < params%chemistry_active_species_threshold) then skipped_named_species_local = skipped_named_species_local + 1 cycle end if end if if (params%chemistry_min_reactive_mass_fraction > zero) then reactive_indicator = maxval(max(zero, species_Y(1:params%nspecies, c))) if (reactive_indicator < params%chemistry_min_reactive_mass_fraction) then skipped_legacy_mass_fraction_local = skipped_legacy_mass_fraction_local + 1 cycle end if end if global_active = global_active + 1 if (mod(global_active - 1, flow%nprocs) == flow%rank) then nactive = nactive + 1 active_cell(nactive) = c end if end do else nloc = max(0, flow%last_cell - flow%first_cell + 1) allocate(active_cell(max(1, nloc))) nactive = 0 skipped_temperature_local = 0 skipped_named_species_local = 0 skipped_legacy_mass_fraction_local = 0 if (nloc > 0) then do i = 1, nloc c = flow%first_cell + i - 1 if (params%chemistry_temperature_cutoff > zero) then if (energy%T(c) < params%chemistry_temperature_cutoff) then skipped_temperature_local = skipped_temperature_local + 1 cycle end if end if if (use_named_species_screen) then reactive_indicator = zero do k = 1, nactive_screen_species reactive_indicator = max(reactive_indicator, & max(zero, species_Y(active_screen_species(k), c))) end do if (reactive_indicator < params%chemistry_active_species_threshold) then skipped_named_species_local = skipped_named_species_local + 1 cycle end if end if if (params%chemistry_min_reactive_mass_fraction > zero) then reactive_indicator = maxval(max(zero, species_Y(1:params%nspecies, c))) if (reactive_indicator < params%chemistry_min_reactive_mass_fraction) then skipped_legacy_mass_fraction_local = skipped_legacy_mass_fraction_local + 1 cycle end if end if nactive = nactive + 1 active_cell(nactive) = c end do end if call MPI_Allreduce(nactive, global_active, 1, MPI_INTEGER, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allreduce failed for active chemistry count') end if if (global_active <= 0) then total_time_local = real(MPI_Wtime(), rk) - t_total_start call accumulate_chemistry_screening_stats(nloc, nactive, skipped_temperature_local, & skipped_named_species_local, skipped_legacy_mass_fraction_local, & total_time_local, reactor_time_local) call write_chemistry_diagnostics_row(flow, params, step, dt_local, nloc, nactive, nloc - nactive, & total_time_local, reactor_time_local, zero, zero, zero, 1.0_rk, & zero, zero, zero, zero, zero) deallocate(active_cell) return end if if (nactive > 0) then allocate(T_local(nactive), P_local(nactive), Y_local(params%nspecies, nactive)) allocate(T_before(nactive), rho_before(nactive), Y_before(params%nspecies, nactive)) allocate(h_local(nactive), cp_local(nactive), lambda_local(nactive), rho_local(nactive), alpha_local(nactive), qdot_local(nactive)) alpha_local = 1.0_rk P_local = params%background_press do i = 1, nactive c = active_cell(i) T_local(i) = energy%T(c) sum_Y = zero do k = 1, params%nspecies Y_local(k, i) = max(zero, species_Y(k, c)) sum_Y = sum_Y + Y_local(k, i) end do if (sum_Y > tiny_safe) then Y_local(:, i) = Y_local(:, i) / sum_Y else Y_local(:, i) = zero Y_local(1, i) = 1.0_rk end if T_before(i) = T_local(i) rho_before(i) = max(energy%rho_thermo(c), tiny_safe) Y_before(:, i) = Y_local(:, i) end do call build_c_species_names(params, c_names_flat, n_len) call profiler_start('Chemistry_Cantera_ReactorNet') t_reactor_start = real(MPI_Wtime(), rk) if (params%enable_chemistry_subcycling) then dT_max_est = zero if (allocated(prev_chem_dT_dt)) then do i = 1, nactive c = active_cell(i) dT_max_est = max(dT_max_est, abs(prev_chem_dT_dt(c)) * dt_local) end do end if n_sub = max(1, min(params%max_chemistry_subcycles, & ceiling(dT_max_est / max(params%subcycling_dT_threshold, tiny_safe)))) dt_sub = dt_local / real(n_sub, rk) do j = 1, n_sub call cantera_integrate_const_p_chemistry_c(nactive, dt_sub, T_local, P_local, params%nspecies, & Y_local, h_local, cp_local, lambda_local, rho_local, & params%energy_reference_T, c_names_flat, n_len, & params%chemistry_rtol, params%chemistry_atol, & params%chemistry_max_steps, merge(1, 0, params%chemistry_energy_enabled)) end do else call cantera_integrate_const_p_chemistry_c(nactive, dt_local, T_local, P_local, params%nspecies, & Y_local, h_local, cp_local, lambda_local, rho_local, & params%energy_reference_T, c_names_flat, n_len, & params%chemistry_rtol, params%chemistry_atol, & params%chemistry_max_steps, merge(1, 0, params%chemistry_energy_enabled)) end if reactor_time_local = real(MPI_Wtime(), rk) - t_reactor_start call profiler_stop('Chemistry_Cantera_ReactorNet') raw_max_dT_local = zero raw_max_dY_local = zero raw_max_rel_drho_local = zero do i = 1, nactive raw_max_dT_local = max(raw_max_dT_local, abs(T_local(i) - T_before(i))) raw_max_dY_local = max(raw_max_dY_local, maxval(abs(Y_local(:, i) - Y_before(:, i)))) raw_max_rel_drho_local = max(raw_max_rel_drho_local, & abs(rho_local(i) - rho_before(i)) / max(rho_before(i), tiny_safe)) end do use_source_limiter = params%chemistry_limit_source_update alpha_base = min(1.0_rk, max(zero, params%chemistry_source_relaxation)) min_alpha_local = 1.0_rk qdot_max_local = zero qdot_integral_local = zero if (use_source_limiter) then do i = 1, nactive alpha = alpha_base dT_abs = abs(T_local(i) - T_before(i)) if (params%chemistry_max_dT_per_step > zero .and. dT_abs > params%chemistry_max_dT_per_step) then alpha = min(alpha, params%chemistry_max_dT_per_step / max(dT_abs, tiny_safe)) end if dY_abs = maxval(abs(Y_local(:, i) - Y_before(:, i))) if (params%chemistry_max_dY_per_step > zero .and. dY_abs > params%chemistry_max_dY_per_step) then alpha = min(alpha, params%chemistry_max_dY_per_step / max(dY_abs, tiny_safe)) end if rel_drho = abs(rho_local(i) - rho_before(i)) / max(rho_before(i), tiny_safe) if (params%chemistry_max_rel_rho_change_per_step > zero .and. & rel_drho > params%chemistry_max_rel_rho_change_per_step) then alpha = min(alpha, params%chemistry_max_rel_rho_change_per_step / max(rel_drho, tiny_safe)) end if alpha = max(zero, min(1.0_rk, alpha)) T_local(i) = T_before(i) + alpha * (T_local(i) - T_before(i)) Y_local(:, i) = Y_before(:, i) + alpha * (Y_local(:, i) - Y_before(:, i)) Y_local(:, i) = max(zero, Y_local(:, i)) sum_Y = sum(Y_local(:, i)) if (sum_Y > tiny_safe) then Y_local(:, i) = Y_local(:, i) / sum_Y else Y_local(:, i) = zero Y_local(1, i) = 1.0_rk end if alpha_local(i) = alpha min_alpha_local = min(min_alpha_local, alpha) end do call cantera_update_thermo_c(nactive, T_local, P_local, params%nspecies, Y_local, & h_local, cp_local, lambda_local, rho_local, & params%energy_reference_T, c_names_flat, n_len) end if call cantera_heat_release_rate_c(nactive, T_local, P_local, params%nspecies, Y_local, qdot_local, & c_names_flat, n_len) qdot_max_local = zero qdot_integral_local = zero do i = 1, nactive c = active_cell(i) qdot_max_local = max(qdot_max_local, qdot_local(i)) qdot_integral_local = qdot_integral_local + qdot_local(i) * mesh%cells(c)%volume end do max_dT_local = zero max_dY_local = zero max_rel_drho_local = zero do i = 1, nactive max_dT_local = max(max_dT_local, abs(T_local(i) - T_before(i))) max_dY_local = max(max_dY_local, maxval(abs(Y_local(:, i) - Y_before(:, i)))) max_rel_drho_local = max(max_rel_drho_local, abs(rho_local(i) - rho_before(i)) / max(rho_before(i), tiny_safe)) end do if (params%enable_chemistry_load_balancing) then nfields = 10 + params%nspecies allocate(pack_cell(max(1, nactive))) allocate(pack_data(nfields, max(1, nactive))) do i = 1, nactive c = active_cell(i) pack_cell(i) = c pack_data(1, i) = T_local(i) pack_data(2, i) = h_local(i) pack_data(3, i) = cp_local(i) pack_data(4, i) = lambda_local(i) pack_data(5, i) = rho_local(i) pack_data(6, i) = alpha_local(i) pack_data(7, i) = qdot_local(i) pack_data(8, i) = T_before(i) pack_data(9, i) = rho_before(i) pack_data(10, i) = maxval(abs(Y_local(:, i) - Y_before(:, i))) pack_data(11:10+params%nspecies, i) = Y_local(1:params%nspecies, i) end do allocate(active_counts(flow%nprocs)) allocate(active_displs(flow%nprocs)) allocate(active_matrix_counts(flow%nprocs)) allocate(active_matrix_displs(flow%nprocs)) call MPI_Allgather(nactive, 1, MPI_INTEGER, & active_counts, 1, MPI_INTEGER, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allgather active counts failed') active_displs(1) = 0 active_matrix_counts(1) = active_counts(1) * nfields active_matrix_displs(1) = 0 do r = 2, flow%nprocs active_displs(r) = active_displs(r - 1) + active_counts(r - 1) active_matrix_counts(r) = active_counts(r) * nfields active_matrix_displs(r) = active_matrix_displs(r - 1) + active_matrix_counts(r - 1) end do global_active = sum(active_counts) allocate(all_active_cells(max(1, global_active))) allocate(all_data_new(nfields, max(1, global_active))) ! Gather active cell indices call MPI_Allgatherv(pack_cell, nactive, MPI_INTEGER, & all_active_cells, active_counts, active_displs, MPI_INTEGER, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allgatherv active cells failed') ! Gather packed data call MPI_Allgatherv(pack_data, nactive * nfields, MPI_DOUBLE_PRECISION, & all_data_new, active_matrix_counts, active_matrix_displs, MPI_DOUBLE_PRECISION, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allgatherv active data failed') ! Unpack gathered data directly into global replicated arrays do i = 1, global_active c = all_active_cells(i) energy%T(c) = all_data_new(1, i) energy%h(c) = all_data_new(2, i) energy%cp(c) = all_data_new(3, i) energy%lambda(c) = all_data_new(4, i) energy%rho_thermo(c) = all_data_new(5, i) alpha = all_data_new(6, i) qdot_val = all_data_new(7, i) T_prev = all_data_new(8, i) rho_prev = all_data_new(9, i) dY_max_val = all_data_new(10, i) species_Y(1:params%nspecies, c) = all_data_new(11:10+params%nspecies, i) if (allocated(energy%chem_active)) energy%chem_active(c) = one if (allocated(energy%chem_source_alpha)) energy%chem_source_alpha(c) = alpha if (allocated(energy%chem_dT_dt)) energy%chem_dT_dt(c) = (energy%T(c) - T_prev) / max(dt_local, tiny_safe) if (allocated(energy%chem_dY_max_dt)) energy%chem_dY_max_dt(c) = dY_max_val / max(dt_local, tiny_safe) if (allocated(energy%chem_rel_rho_change_dt)) energy%chem_rel_rho_change_dt(c) = & (abs(energy%rho_thermo(c) - rho_prev) / max(rho_prev, tiny_safe)) / max(dt_local, tiny_safe) if (allocated(energy%chem_heat_release_rate)) energy%chem_heat_release_rate(c) = qdot_val end do ! Enforce mass fraction normalization do i = 1, global_active c = all_active_cells(i) species_Y(1:params%nspecies, c) = max(zero, species_Y(1:params%nspecies, c)) sum_Y = sum(species_Y(1:params%nspecies, c)) if (sum_Y > tiny_safe) then species_Y(1:params%nspecies, c) = species_Y(1:params%nspecies, c) / sum_Y else species_Y(:, c) = zero species_Y(1, c) = 1.0_rk end if end do deallocate(pack_cell, pack_data) deallocate(active_counts, active_displs, active_matrix_counts, active_matrix_displs) deallocate(all_active_cells, all_data_new) else do i = 1, nactive c = active_cell(i) energy%T(c) = T_local(i) energy%h(c) = h_local(i) energy%cp(c) = cp_local(i) energy%lambda(c) = lambda_local(i) energy%rho_thermo(c) = rho_local(i) do k = 1, params%nspecies species_Y(k, c) = max(zero, Y_local(k, i)) end do sum_Y = sum(species_Y(1:params%nspecies, c)) if (sum_Y > tiny_safe) then species_Y(1:params%nspecies, c) = species_Y(1:params%nspecies, c) / sum_Y else species_Y(:, c) = zero species_Y(1, c) = 1.0_rk end if if (allocated(energy%chem_active)) energy%chem_active(c) = one if (allocated(energy%chem_source_alpha)) energy%chem_source_alpha(c) = alpha_local(i) if (allocated(energy%chem_dT_dt)) energy%chem_dT_dt(c) = (T_local(i) - T_before(i)) / max(dt_local, tiny_safe) if (allocated(energy%chem_dY_max_dt)) energy%chem_dY_max_dt(c) = & maxval(abs(Y_local(:, i) - Y_before(:, i))) / max(dt_local, tiny_safe) if (allocated(energy%chem_rel_rho_change_dt)) energy%chem_rel_rho_change_dt(c) = & (abs(rho_local(i) - rho_before(i)) / max(rho_before(i), tiny_safe)) / max(dt_local, tiny_safe) if (allocated(energy%chem_heat_release_rate)) energy%chem_heat_release_rate(c) = qdot_local(i) end do end if end if if (.not. allocated(T_local)) then max_dT_local = zero max_dY_local = zero max_rel_drho_local = zero raw_max_dT_local = zero raw_max_dY_local = zero raw_max_rel_drho_local = zero min_alpha_local = 1.0_rk qdot_max_local = zero qdot_integral_local = zero end if reduce_in = [max_dT_local, max_dY_local, max_rel_drho_local] raw_reduce_in = [raw_max_dT_local, raw_max_dY_local, raw_max_rel_drho_local] call MPI_Allreduce(raw_reduce_in, raw_reduce_out, 3, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allreduce failed for raw chemistry source diagnostics') call MPI_Allreduce(reduce_in, reduce_out, 3, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allreduce failed for chemistry source diagnostics') call MPI_Allreduce(min_alpha_local, alpha_reduce, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allreduce failed for chemistry source limiter alpha') if (present(max_dT_chem)) max_dT_chem = reduce_out(1) if (present(max_dY_chem)) max_dY_chem = reduce_out(2) if (present(max_rel_drho_chem)) max_rel_drho_chem = reduce_out(3) if (present(min_source_alpha)) min_source_alpha = alpha_reduce if (present(raw_max_dT_chem)) raw_max_dT_chem = raw_reduce_out(1) if (present(raw_max_dY_chem)) raw_max_dY_chem = raw_reduce_out(2) if (present(raw_max_rel_drho_chem)) raw_max_rel_drho_chem = raw_reduce_out(3) if (.not. params%enable_chemistry_load_balancing) then call flow_exchange_cell_matrix(flow, species_Y) call flow_exchange_cell_scalars(flow, energy%T, energy%h, energy%cp, energy%lambda) call flow_exchange_cell_scalar(flow, energy%rho_thermo) end if total_time_local = real(MPI_Wtime(), rk) - t_total_start call accumulate_chemistry_screening_stats(nloc, nactive, skipped_temperature_local, & skipped_named_species_local, skipped_legacy_mass_fraction_local, & total_time_local, reactor_time_local, nactive * n_sub) call write_chemistry_diagnostics_row(flow, params, step, dt_local, nloc, nactive, nloc - nactive, & total_time_local, reactor_time_local, & max_dT_local, max_dY_local, max_rel_drho_local, min_alpha_local, & raw_max_dT_local, raw_max_dY_local, raw_max_rel_drho_local, & qdot_max_local, qdot_integral_local) if (allocated(T_local)) deallocate(T_local, P_local, Y_local, T_before, Y_before, rho_before, & h_local, cp_local, lambda_local, rho_local, alpha_local, qdot_local) if (allocated(c_names_flat)) deallocate(c_names_flat) deallocate(active_cell) end subroutine advance_cantera_chemistry !> Advance transported sensible enthalpy one explicit finite-volume step. !! !! The explicit finite-volume update is: !! !! V dh/dt = - sum_f F_f h_f !! + (1/rho) sum_f lambda A_f (T_nb - T_c)/d_f !! + (qrad/rho) V !! !! `fields%face_flux` is volumetric flux in \( \mathrm{m^3/s} \). !! `fields%mass_flux` is density-weighted flux in \( \mathrm{kg/s} \). !! Both are stored owner-to-neighbor and re-oriented outward from the !! currently updated cell before upwind advection. !! !! Constant-density mode updates `h` with `params%rho`. Variable-density !! mode uses a conservative `rho*h` update with `fields%mass_flux`, !! `transport%rho_old`, and current `transport%rho`, then stores the !! transported specific enthalpy back in `energy%h`. subroutine advance_energy_transport(mesh, flow, bc, params, fields, energy, species_Y, transport, step) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(inout) :: flow type(bc_set_t), intent(in) :: bc type(case_params_t), intent(in) :: params type(flow_fields_t), intent(in) :: fields type(energy_fields_t), intent(inout) :: energy real(rk), intent(in), optional :: species_Y(:,:) type(transport_properties_t), intent(in), optional :: transport integer, intent(in), optional :: step integer :: c, lf, fid, neighbor, k, ierr real(rk) :: flux_out real(rk) :: face_area, dist, vol real(rk) :: h_cell, h_other, h_face real(rk) :: T_cell, T_other real(rk) :: rhs_h, rhs_h_current, rhs_h_update, diff_term, lambda_face real(rk) :: c_curr, c_old real(rk) :: conductive_boundary_flux_out_local real(rk) :: advective_boundary_rho_h_flux_out_local real(rk) :: qrad_integral_step_local real(rk) :: rho_species_hdiff_integral_step_local real(rk) :: rho_diag real(rk) :: energy_update_delta_rate_local, energy_update_rhs_integral_local real(rk) :: energy_update_delta_rate_global, energy_update_rhs_integral_global real(rk) :: energy_update_balance_global, energy_update_denom real(rk) :: update_send(3), update_recv(3) real(rk) :: operator_consistent_rho_h_integral_local real(rk) :: operator_consistent_rho_h_integral_global real(rk) :: rho_cell, rho_old_cell, rho_face real(rk) :: sum_diff_flux, diff_face, species_hdiff_rhs real(rk) :: Y_cell, Y_other logical :: is_dirichlet, do_diffusion, use_species_hdiff, has_species_diffusion_face logical :: use_variable_density_energy, use_ab2 real(rk), allocatable :: Y_point(:) ! Explicitly null-initialized pointer scratch avoids false gfortran ! release-build descriptor warnings on the conditional h_k diffusion path. real(rk), pointer :: hk_species(:,:) => null(), hk_boundary(:) => null() real(rk), allocatable :: sp_diff_flux(:), sp_Y_face_lin(:) real(rk), allocatable :: grad_h(:,:) if (.not. params%enable_energy) return if (.not. energy%initialized) then call fatal_error('energy', 'advance requested before energy initialization') end if if (params%rho <= tiny_safe) call fatal_error('energy', 'rho must be positive') if (params%energy_cp <= tiny_safe) call fatal_error('energy', 'energy_cp must be positive') if (params%energy_lambda < zero) call fatal_error('energy', 'energy_lambda must be non-negative') use_ab2 = trim(params%energy_time_scheme) == 'ab2' .or. & trim(params%energy_time_scheme) == 'adams_bashforth2' .or. & trim(params%energy_time_scheme) == 'adams-bashforth2' c_curr = 1.5_rk c_old = -0.5_rk if (use_ab2 .and. energy%rhs_old_valid .and. params%dt_old > tiny_safe) then c_curr = 1.0_rk + 0.5_rk * (params%dt / params%dt_old) c_old = -0.5_rk * (params%dt / params%dt_old) end if use_variable_density_energy = params%enable_variable_density if (use_variable_density_energy) then if (.not. present(transport)) then call fatal_error('energy', 'variable-density energy transport requires transport properties') end if if (.not. allocated(transport%rho)) then call fatal_error('energy', 'variable-density energy transport requires transport%rho') end if if (.not. allocated(transport%rho_old)) then call fatal_error('energy', 'variable-density energy transport requires transport%rho_old') end if if (.not. allocated(fields%mass_flux)) then call fatal_error('energy', 'variable-density energy transport requires fields%mass_flux') end if end if use_species_hdiff = params%enable_species_enthalpy_diffusion if (use_species_hdiff) then if (.not. present(species_Y)) then call fatal_error('energy', 'species enthalpy diffusion requires species_Y') end if if (.not. present(transport)) then call fatal_error('energy', 'species enthalpy diffusion requires transport properties') end if if (.not. allocated(transport%diffusivity)) then call fatal_error('energy', 'species enthalpy diffusion requires transport diffusivity') end if if (.not. params%enable_cantera_thermo) then call fatal_error('energy', 'species enthalpy diffusion requires Cantera thermo') end if if (params%nspecies <= 0) then call fatal_error('energy', 'species enthalpy diffusion requires nspecies > 0') end if if (size(species_Y, 1) < params%nspecies .or. size(species_Y, 2) < mesh%ncells) then call fatal_error('energy', 'species enthalpy diffusion species_Y has incompatible shape') end if if (size(transport%diffusivity, 1) < params%nspecies .or. & size(transport%diffusivity, 2) < mesh%ncells) then call fatal_error('energy', 'species enthalpy diffusion diffusivity has incompatible shape') end if end if ! Option A: preserve transported enthalpy across species updates. ! Species transport may have changed Y before this energy step. Do not ! recompute h from the old temperature and the new composition. Instead, ! keep h as the transported state, recover T(h,Y,p0), then refresh the ! thermo/transport properties used by the enthalpy equation. call profiler_start('Energy_Exchange_H') call flow_exchange_cell_scalar(flow, energy%h) call profiler_stop('Energy_Exchange_H') energy%h_old = energy%h if (params%enable_cantera_thermo .and. params%enable_species .and. & present(species_Y) .and. params%nspecies > 0) then allocate(Y_point(params%nspecies)) end if if (use_species_hdiff) then if (.not. allocated(Y_point)) allocate(Y_point(params%nspecies)) allocate(hk_species(params%nspecies, mesh%ncells)) allocate(hk_boundary(params%nspecies)) allocate(sp_diff_flux(params%nspecies)) allocate(sp_Y_face_lin(params%nspecies)) hk_species = zero hk_boundary = zero sp_diff_flux = zero sp_Y_face_lin = zero end if if (params%enable_cantera_thermo) then ! A pre-flux thermo sync is required only when species transport may ! have changed the composition since the previous post-flux sync. ! Without species, current T/cp/lambda/rho_thermo are already valid ! from initialization or the previous energy step's post-sync. if (params%enable_species .and. present(species_Y) .and. params%nspecies > 0) then call write_variable_density_energy_debug(mesh, flow, params, fields, energy, & 'pre_cantera_presync', species_Y, transport, step) if (use_species_hdiff) then call profiler_start('Energy_Cantera_PreSyncSpeciesH') call recover_temperature_update_thermo_and_species_h_cantera_owned(mesh, flow, params, energy, species_Y, hk_species) call profiler_stop('Energy_Cantera_PreSyncSpeciesH') else call profiler_start('Energy_Cantera_PreSync') call recover_temperature_and_update_thermo_cantera_owned(mesh, flow, params, energy, species_Y) call profiler_stop('Energy_Cantera_PreSync') end if energy%h = energy%h_old end if end if ! Ensure off-rank temperature/property values are current before fluxes. call profiler_start('Energy_PreFlux_Exchange') if (allocated(energy%lambda)) then call flow_exchange_cell_scalars(flow, energy%T, energy%lambda) else call flow_exchange_cell_scalar(flow, energy%T) end if call profiler_stop('Energy_PreFlux_Exchange') if (allocated(energy%T_old)) energy%T_old = energy%T allocate(grad_h(3, mesh%ncells)) call compute_scalar_gradients(mesh, energy%h_old, grad_h) if (use_species_hdiff .and. .not. associated(hk_species)) then call fatal_error('energy', 'species enthalpy diffusion scratch was not initialized') end if if (allocated(energy%species_enthalpy_diffusion)) energy%species_enthalpy_diffusion = zero conductive_boundary_flux_out_local = zero energy_update_delta_rate_local = zero energy_update_rhs_integral_local = zero operator_consistent_rho_h_integral_local = zero if (allocated(energy%operator_consistent_rho_h)) energy%operator_consistent_rho_h = zero energy%operator_consistent_rho_h_available = 0 advective_boundary_rho_h_flux_out_local = zero qrad_integral_step_local = zero rho_species_hdiff_integral_step_local = zero call profiler_start('Energy_Flux_Update') do c = flow%first_cell, flow%last_cell rhs_h = zero h_cell = energy%h_old(c) T_cell = energy%T(c) do lf = 1, mesh%ncell_faces(c) fid = mesh%cell_faces(lf, c) ! `fields%face_flux` and `fields%mass_flux` are oriented ! owner -> neighbor. Reorient outward from the current cell. ! Variable-density energy uses mass flux for conservative rho*h ! advection; constant-density energy uses volumetric flux. if (mesh%faces(fid)%owner == c) then if (use_variable_density_energy) then flux_out = fields%mass_flux(fid) else flux_out = fields%face_flux(fid) end if else if (use_variable_density_energy) then flux_out = -fields%mass_flux(fid) else flux_out = -fields%face_flux(fid) end if end if neighbor = face_effective_neighbor(mesh, bc, fid, c) face_area = mesh%faces(fid)%area is_dirichlet = .false. if (neighbor > 0) then h_other = energy%h_old(neighbor) T_other = energy%T(neighbor) do_diffusion = .true. else call boundary_temperature(mesh, bc, fid, T_cell, T_other, is_dirichlet) if (is_dirichlet) then if (allocated(Y_point)) then call build_boundary_thermo_Y(mesh, bc, params, fid, c, species_Y, Y_point) h_other = boundary_enthalpy_from_temperature(mesh, params, T_other, Y_point) else h_other = boundary_enthalpy_from_temperature(mesh, params, T_other) end if do_diffusion = .true. else h_other = h_cell T_other = T_cell do_diffusion = .false. end if end if ! Advective h flux using the selected scalar scheme. The flux ! is outward from the current cell and is volumetric in the ! constant-density path, mass-weighted in variable-density mode. h_face = scalar_face_value(mesh, energy%h_old, grad_h, & params%energy_convection_scheme, params%scalar_limiter, & c, neighbor, fid, flux_out, h_other, is_dirichlet) rhs_h = rhs_h - flux_out * h_face if (neighbor <= 0) then if (use_variable_density_energy) then advective_boundary_rho_h_flux_out_local = advective_boundary_rho_h_flux_out_local + flux_out * h_face else advective_boundary_rho_h_flux_out_local = advective_boundary_rho_h_flux_out_local + params%rho * flux_out * h_face end if end if ! Fourier heat conduction uses grad(T), not grad(h). if (do_diffusion .and. (params%enable_cantera_thermo .or. params%energy_lambda > zero)) then dist = energy_face_normal_distance(mesh, bc, fid, c, neighbor) if (allocated(energy%lambda)) then if (neighbor > 0) then lambda_face = 0.5_rk * (energy%lambda(c) + energy%lambda(neighbor)) else lambda_face = energy%lambda(c) end if else lambda_face = params%energy_lambda end if if (lambda_face > zero) then if (use_variable_density_energy) then diff_term = lambda_face * (T_other - T_cell) / dist * face_area else diff_term = (lambda_face / params%rho) * & (T_other - T_cell) / dist * face_area end if rhs_h = rhs_h + diff_term if (neighbor <= 0) then if (use_variable_density_energy) then conductive_boundary_flux_out_local = conductive_boundary_flux_out_local - diff_term else conductive_boundary_flux_out_local = conductive_boundary_flux_out_local - params%rho * diff_term end if end if end if end if ! Species-enthalpy diffusion correction: ! corrected_diff_flux(k) matches mod_species and equals -J_k,out*A/rho. ! Therefore sum_k h_k,face * corrected_diff_flux(k) is added to the ! sensible-enthalpy RHS for this cell. if (use_species_hdiff) then dist = energy_face_normal_distance(mesh, bc, fid, c, neighbor) sum_diff_flux = zero has_species_diffusion_face = .false. do k = 1, params%nspecies Y_cell = species_Y(k, c) if (neighbor > 0) then Y_other = species_Y(k, neighbor) is_dirichlet = .true. else call boundary_species(mesh, bc, fid, k, Y_cell, Y_other, is_dirichlet) end if sp_diff_flux(k) = zero sp_Y_face_lin(k) = 0.5_rk * (Y_cell + Y_other) if (neighbor /= 0 .or. is_dirichlet) then has_species_diffusion_face = .true. if (neighbor == 0) then diff_face = transport%diffusivity(k, c) else diff_face = 0.5_rk * (transport%diffusivity(k, c) + transport%diffusivity(k, neighbor)) end if if (use_variable_density_energy) then if (neighbor == 0) then rho_face = transport%rho(c) else rho_face = 0.5_rk * (transport%rho(c) + transport%rho(neighbor)) end if sp_diff_flux(k) = rho_face * diff_face * (Y_other - Y_cell) / dist * face_area else sp_diff_flux(k) = diff_face * (Y_other - Y_cell) / dist * face_area end if end if sum_diff_flux = sum_diff_flux + sp_diff_flux(k) end do species_hdiff_rhs = zero if (has_species_diffusion_face) then if (neighbor > 0) then do k = 1, params%nspecies species_hdiff_rhs = species_hdiff_rhs + & 0.5_rk * (hk_species(k, c) + hk_species(k, neighbor)) * & (sp_diff_flux(k) - sp_Y_face_lin(k) * sum_diff_flux) end do else call build_boundary_thermo_Y(mesh, bc, params, fid, c, species_Y, Y_point) call species_sensible_enthalpies_from_temperature_cantera_point(params, T_other, hk_boundary, Y_point) do k = 1, params%nspecies species_hdiff_rhs = species_hdiff_rhs + & 0.5_rk * (hk_species(k, c) + hk_boundary(k)) * & (sp_diff_flux(k) - sp_Y_face_lin(k) * sum_diff_flux) end do end if end if rhs_h = rhs_h + species_hdiff_rhs if (allocated(energy%species_enthalpy_diffusion)) then if (use_variable_density_energy) then rho_cell = max(transport%rho(c), tiny_safe) energy%species_enthalpy_diffusion(c) = energy%species_enthalpy_diffusion(c) + & species_hdiff_rhs / & (rho_cell * mesh%cells(c)%volume) else energy%species_enthalpy_diffusion(c) = energy%species_enthalpy_diffusion(c) + & species_hdiff_rhs / mesh%cells(c)%volume end if end if end if end do if (use_variable_density_energy) then rho_cell = max(transport%rho(c), tiny_safe) rho_old_cell = max(transport%rho_old(c), tiny_safe) rhs_h_current = rhs_h / mesh%cells(c)%volume + energy%qrad(c) if (use_ab2 .and. energy%rhs_old_valid) then rhs_h_update = c_curr * rhs_h_current + c_old * energy%rhs_old(c) else rhs_h_update = rhs_h_current end if energy%h(c) = (rho_old_cell * energy%h_old(c) + params%dt * rhs_h_update) / rho_cell energy_update_delta_rate_local = energy_update_delta_rate_local + & ((rho_cell * energy%h(c) - rho_old_cell * energy%h_old(c)) * mesh%cells(c)%volume) / params%dt energy_update_rhs_integral_local = energy_update_rhs_integral_local + & rhs_h_update * mesh%cells(c)%volume operator_consistent_rho_h_integral_local = operator_consistent_rho_h_integral_local + & rho_cell * energy%h(c) * mesh%cells(c)%volume if (allocated(energy%operator_consistent_rho_h)) then energy%operator_consistent_rho_h(c) = rho_cell * energy%h(c) end if else rhs_h_current = rhs_h / mesh%cells(c)%volume + energy%qrad(c) / params%rho if (use_ab2 .and. energy%rhs_old_valid) then rhs_h_update = c_curr * rhs_h_current + c_old * energy%rhs_old(c) else rhs_h_update = rhs_h_current end if energy%h(c) = energy%h_old(c) + params%dt * rhs_h_update energy_update_delta_rate_local = energy_update_delta_rate_local + & params%rho * (energy%h(c) - energy%h_old(c)) * mesh%cells(c)%volume / params%dt energy_update_rhs_integral_local = energy_update_rhs_integral_local + & params%rho * rhs_h_update * mesh%cells(c)%volume operator_consistent_rho_h_integral_local = operator_consistent_rho_h_integral_local + & params%rho * energy%h(c) * mesh%cells(c)%volume if (allocated(energy%operator_consistent_rho_h)) then energy%operator_consistent_rho_h(c) = params%rho * energy%h(c) end if end if energy%rhs_old(c) = rhs_h_current end do call profiler_stop('Energy_Flux_Update') energy%rhs_old_valid = .true. update_send(1) = energy_update_delta_rate_local update_send(2) = energy_update_rhs_integral_local update_send(3) = operator_consistent_rho_h_integral_local call MPI_Allreduce(update_send, update_recv, 3, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing direct energy update closure') energy_update_delta_rate_global = update_recv(1) energy_update_rhs_integral_global = update_recv(2) operator_consistent_rho_h_integral_global = update_recv(3) energy_update_balance_global = energy_update_delta_rate_global - energy_update_rhs_integral_global energy_update_denom = max(abs(energy_update_delta_rate_global), abs(energy_update_rhs_integral_global), tiny_safe) energy%last_energy_update_delta_rate_integral = energy_update_delta_rate_global energy%last_energy_update_rhs_integral = energy_update_rhs_integral_global energy%last_energy_update_balance_defect = energy_update_balance_global energy%relative_last_energy_update_balance_defect = abs(energy_update_balance_global) / energy_update_denom energy%cumulative_energy_update_delta_integral = & energy%cumulative_energy_update_delta_integral + params%dt * energy_update_delta_rate_global energy%cumulative_energy_update_rhs_integral = & energy%cumulative_energy_update_rhs_integral + params%dt * energy_update_rhs_integral_global energy%last_operator_consistent_rho_h_integral = operator_consistent_rho_h_integral_global energy%operator_consistent_rho_h_available = 1 energy%last_conductive_boundary_flux_out = conductive_boundary_flux_out_local energy%last_conductive_boundary_flux_available = 1 qrad_integral_step_local = zero rho_species_hdiff_integral_step_local = zero do c = flow%first_cell, flow%last_cell vol = mesh%cells(c)%volume if (allocated(energy%qrad)) qrad_integral_step_local = qrad_integral_step_local + energy%qrad(c) * vol if (allocated(energy%species_enthalpy_diffusion)) then if (use_variable_density_energy) then rho_diag = max(transport%rho(c), tiny_safe) else rho_diag = params%rho end if rho_species_hdiff_integral_step_local = rho_species_hdiff_integral_step_local + & rho_diag * energy%species_enthalpy_diffusion(c) * vol end if end do energy%cumulative_boundary_rho_h_advective_flux_out = & energy%cumulative_boundary_rho_h_advective_flux_out + params%dt * advective_boundary_rho_h_flux_out_local energy%cumulative_boundary_rho_h_conductive_flux_out = & energy%cumulative_boundary_rho_h_conductive_flux_out + params%dt * conductive_boundary_flux_out_local energy%cumulative_qrad_integral = energy%cumulative_qrad_integral + params%dt * qrad_integral_step_local energy%cumulative_rho_species_hdiff_integral = & energy%cumulative_rho_species_hdiff_integral + params%dt * rho_species_hdiff_integral_step_local energy%cumulative_energy_budget_available = 1 call write_variable_density_energy_debug(mesh, flow, params, fields, energy, & 'post_flux_pre_cantera_postsync', species_Y, transport, step) if (allocated(grad_h)) deallocate(grad_h) if (associated(hk_species)) deallocate(hk_species) if (associated(hk_boundary)) deallocate(hk_boundary) if (allocated(sp_diff_flux)) deallocate(sp_diff_flux) if (allocated(sp_Y_face_lin)) deallocate(sp_Y_face_lin) if (allocated(Y_point)) deallocate(Y_point) if (params%enable_cantera_thermo) then ! Protect transported h from property-refresh roundoff: recover T from ! the newly transported h, refresh cp/lambda/rho_thermo, then restore h. energy%h_old = energy%h call profiler_start('Energy_Cantera_PostSync') call recover_temperature_and_update_thermo_cantera_owned(mesh, flow, params, energy, species_Y) call profiler_stop('Energy_Cantera_PostSync') energy%h = energy%h_old else call profiler_start('Energy_ConstantCp_RecoverT') call recover_temperature_constant_cp(params, energy) call profiler_stop('Energy_ConstantCp_RecoverT') end if ! Synchronize updated owned cells for output and the next step. ! Pack related scalar fields to reduce neighbor-message latency. call profiler_start('Energy_Final_Exchange') if (allocated(energy%lambda) .and. allocated(energy%species_enthalpy_diffusion)) then call flow_exchange_cell_scalars(flow, energy%h, energy%T, energy%lambda, energy%species_enthalpy_diffusion) else if (allocated(energy%lambda)) then call flow_exchange_cell_scalars(flow, energy%h, energy%T, energy%lambda) else if (allocated(energy%species_enthalpy_diffusion)) then call flow_exchange_cell_scalars(flow, energy%h, energy%T, energy%species_enthalpy_diffusion) else call flow_exchange_cell_scalars(flow, energy%h, energy%T) end if call profiler_stop('Energy_Final_Exchange') end subroutine advance_energy_transport !> Targeted diagnostics for experimental variable-density energy/Cantera HP failures. subroutine write_variable_density_energy_debug(mesh, flow, params, fields, energy, label, species_Y, transport, step) 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(energy_fields_t), intent(in) :: energy character(len=*), intent(in) :: label real(rk), intent(in), optional :: species_Y(:,:) type(transport_properties_t), intent(in), optional :: transport integer, intent(in), optional :: step integer :: c, k, ierr integer :: local_bad_cell real(rk) :: local_min_h, local_max_h, local_min_T, local_max_T real(rk) :: local_min_rho, local_max_rho, local_min_rho_old, local_max_rho_old real(rk) :: local_min_S, local_max_S real(rk) :: global_min_h, global_max_h, global_min_T, global_max_T real(rk) :: global_min_rho, global_max_rho, global_min_rho_old, global_max_rho_old real(rk) :: global_min_S, global_max_S real(rk) :: local_min_vec(5), local_max_vec(5) real(rk) :: global_min_vec(5), global_max_vec(5) logical :: have_transport, suspicious if (.not. params%enable_variable_density) return if (.not. params%variable_density_debug) return if (present(step)) then if (mod(step, max(1, params%variable_density_debug_interval)) /= 0) return end if if (.not. allocated(energy%h) .or. .not. allocated(energy%T)) return associate(dummy_mesh_ncells => mesh%ncells) end associate have_transport = present(transport) local_min_h = huge(zero) local_max_h = -huge(zero) local_min_T = huge(zero) local_max_T = -huge(zero) local_min_rho = huge(zero) local_max_rho = -huge(zero) local_min_rho_old = huge(zero) local_max_rho_old = -huge(zero) local_min_S = huge(zero) local_max_S = -huge(zero) local_bad_cell = -1 do c = flow%first_cell, flow%last_cell local_min_h = min(local_min_h, energy%h(c)) local_max_h = max(local_max_h, energy%h(c)) local_min_T = min(local_min_T, energy%T(c)) local_max_T = max(local_max_T, energy%T(c)) if (have_transport) then if (allocated(transport%rho)) then local_min_rho = min(local_min_rho, transport%rho(c)) local_max_rho = max(local_max_rho, transport%rho(c)) end if if (allocated(transport%rho_old)) then local_min_rho_old = min(local_min_rho_old, transport%rho_old(c)) local_max_rho_old = max(local_max_rho_old, transport%rho_old(c)) end if end if if (allocated(fields%divergence_source)) then local_min_S = min(local_min_S, fields%divergence_source(c)) local_max_S = max(local_max_S, fields%divergence_source(c)) end if if (local_bad_cell < 0) then if (energy%h(c) < -1.0e5_rk .or. energy%h(c) > 1.0e7_rk) then local_bad_cell = c end if end if end do if (.not. have_transport) then local_min_rho = zero local_max_rho = zero local_min_rho_old = zero local_max_rho_old = zero end if if (.not. allocated(fields%divergence_source)) then local_min_S = zero local_max_S = zero end if local_min_vec = [local_min_h, local_min_T, local_min_rho, local_min_rho_old, local_min_S] local_max_vec = [local_max_h, local_max_T, local_max_rho, local_max_rho_old, local_max_S] call MPI_Allreduce(local_min_vec, global_min_vec, 5, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing vd debug minima') call MPI_Allreduce(local_max_vec, global_max_vec, 5, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing vd debug maxima') global_min_h = global_min_vec(1) global_min_T = global_min_vec(2) global_min_rho = global_min_vec(3) global_min_rho_old = global_min_vec(4) global_min_S = global_min_vec(5) global_max_h = global_max_vec(1) global_max_T = global_max_vec(2) global_max_rho = global_max_vec(3) global_max_rho_old = global_max_vec(4) global_max_S = global_max_vec(5) if (flow%rank == 0) then write(output_unit,'(a,a,a,2(1x,es13.5),a,2(1x,es13.5),a,2(1x,es13.5),a,2(1x,es13.5),a,2(1x,es13.5))') & 'vd-energy-debug[', trim(label), '] h=', global_min_h, global_max_h, & ' T=', global_min_T, global_max_T, & ' rho=', global_min_rho, global_max_rho, & ' rho_old=', global_min_rho_old, global_max_rho_old, & ' S=', global_min_S, global_max_S end if suspicious = local_bad_cell > 0 if (suspicious) then c = local_bad_cell write(output_unit,'(a,i0,a,i0,a,a,a,es16.8,a,es16.8)') & 'vd-energy-bad-cell rank=', flow%rank, ' cell=', c, ' label=', trim(label), & ' h=', energy%h(c), ' T=', energy%T(c) if (have_transport) then if (allocated(transport%rho) .and. allocated(transport%rho_old)) then write(output_unit,'(a,es16.8,a,es16.8)') & ' rho=', transport%rho(c), ' rho_old=', transport%rho_old(c) end if end if if (allocated(fields%divergence_source)) then write(output_unit,'(a,es16.8)') ' divergence_source=', fields%divergence_source(c) end if if (present(species_Y)) then do k = 1, min(params%nspecies, size(species_Y, 1)) write(output_unit,'(a,i0,a,es16.8)') ' Y(', k, ')=', species_Y(k, c) end do end if end if end subroutine write_variable_density_energy_debug !> Outward unit normal from cell_id for face_id. function energy_outward_normal(mesh, face_id, cell_id) result(nvec) type(mesh_t), intent(in) :: mesh integer, intent(in) :: face_id integer, intent(in) :: cell_id real(rk) :: nvec(3) if (mesh%faces(face_id)%owner == cell_id) then nvec = mesh%faces(face_id)%normal else nvec = -mesh%faces(face_id)%normal end if end function energy_outward_normal !> Normal distance used for cell-cell and cell-boundary temperature gradients. function energy_face_normal_distance(mesh, bc, face_id, cell_id, nb) result(dist) type(mesh_t), intent(in) :: mesh type(bc_set_t), intent(in) :: bc integer, intent(in) :: face_id integer, intent(in) :: cell_id integer, intent(in) :: nb real(rk) :: dist integer :: pair_face integer :: btype real(rk) :: nvec(3) nvec = energy_outward_normal(mesh, face_id, cell_id) if (nb > 0) then if (mesh%faces(face_id)%neighbor == 0) then btype = patch_type_for_face(mesh, bc, face_id) if (btype == bc_periodic) then pair_face = mesh%faces(face_id)%periodic_face if (pair_face <= 0) then call fatal_error('energy', 'periodic face has no paired face') end if dist = abs(dot_product(mesh%faces(face_id)%center - & mesh%cells(cell_id)%center, nvec)) + & abs(dot_product(mesh%cells(nb)%center - & mesh%faces(pair_face)%center, nvec)) dist = max(dist, tiny_safe) return end if end if dist = abs(dot_product(mesh%cells(nb)%center - & mesh%cells(cell_id)%center, nvec)) else dist = abs(dot_product(mesh%faces(face_id)%center - & mesh%cells(cell_id)%center, nvec)) end if dist = max(dist, tiny_safe) end function energy_face_normal_distance !> Writes the CSV header for energy diagnostics. subroutine write_energy_diagnostics_header(params, flow) type(case_params_t), intent(in) :: params type(flow_mpi_t), intent(in) :: flow integer :: unit_id character(len=256 + 32) :: filename if (flow%rank /= 0 .or. .not. params%write_diagnostics) return if (.not. params%enable_energy) return filename = trim(params%output_dir)//'/diagnostics/energy_diagnostics.csv' open(newunit=unit_id, file=trim(filename), status='replace', action='write') write(unit_id,'(a)') 'step,time,min_T,max_T,mean_T,min_h,max_h,mean_h,min_qrad,max_qrad,integral_qrad,max_delta_T,rel_h_residual' close(unit_id) end subroutine write_energy_diagnostics_header !> Appends one row of global energy diagnostics. !! !! The temperature and enthalpy means are volume weighted over owned cells. !! qrad integral is the domain integral of qrad dV [W]. subroutine write_energy_diagnostics_row(mesh, flow, params, energy, step, time, write_file) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(case_params_t), intent(in) :: params type(energy_fields_t), intent(in) :: energy integer, intent(in) :: step real(rk), intent(in) :: time logical, intent(in), optional :: write_file integer :: c, unit_id, ierr character(len=256 + 32) :: filename real(rk) :: vol, old_T real(rk) :: local_min_T, local_max_T, local_sum_T real(rk) :: local_min_h, local_max_h, local_sum_h real(rk) :: local_min_qrad, local_max_qrad, local_integral_qrad real(rk) :: local_volume real(rk) :: local_max_delta_T, local_delta_h2, local_h_old2 real(rk) :: global_min_T, global_max_T, global_sum_T real(rk) :: global_min_h, global_max_h, global_sum_h real(rk) :: global_min_qrad, global_max_qrad, global_integral_qrad real(rk) :: global_volume real(rk) :: global_max_delta_T, global_delta_h2, global_h_old2 real(rk) :: mean_T, mean_h, rel_h_residual logical :: do_write_file if (.not. params%enable_energy) return ! File diagnostics only. Terminal health output is handled in main.f90. do_write_file = params%write_diagnostics if (present(write_file)) do_write_file = write_file .and. params%write_diagnostics if (.not. allocated(energy%T) .or. .not. allocated(energy%h) .or. & .not. allocated(energy%h_old) .or. .not. allocated(energy%qrad)) then call fatal_error('energy', 'energy diagnostics requested but energy arrays are not allocated') end if local_min_T = huge(0.0_rk) local_max_T = -huge(0.0_rk) local_sum_T = zero local_min_h = huge(0.0_rk) local_max_h = -huge(0.0_rk) local_sum_h = zero local_min_qrad = huge(0.0_rk) local_max_qrad = -huge(0.0_rk) local_integral_qrad = zero local_volume = zero local_max_delta_T = zero local_delta_h2 = zero local_h_old2 = zero do c = 1, mesh%ncells if (.not. flow%owned(c)) cycle vol = mesh%cells(c)%volume local_volume = local_volume + vol local_min_T = min(local_min_T, energy%T(c)) local_max_T = max(local_max_T, energy%T(c)) local_sum_T = local_sum_T + energy%T(c) * vol local_min_h = min(local_min_h, energy%h(c)) local_max_h = max(local_max_h, energy%h(c)) local_sum_h = local_sum_h + energy%h(c) * vol local_min_qrad = min(local_min_qrad, energy%qrad(c)) local_max_qrad = max(local_max_qrad, energy%qrad(c)) local_integral_qrad = local_integral_qrad + energy%qrad(c) * vol if (allocated(energy%T_old)) then old_T = energy%T_old(c) else old_T = params%energy_reference_T + & (energy%h_old(c) - params%energy_reference_h) / params%energy_cp end if local_max_delta_T = max(local_max_delta_T, abs(energy%T(c) - old_T)) local_delta_h2 = local_delta_h2 + (energy%h(c) - energy%h_old(c))**2 * vol local_h_old2 = local_h_old2 + energy%h_old(c)**2 * vol end do call MPI_Allreduce(local_min_T, global_min_T, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing min_T') call MPI_Allreduce(local_max_T, global_max_T, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing max_T') call MPI_Allreduce(local_sum_T, global_sum_T, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing sum_T') call MPI_Allreduce(local_min_h, global_min_h, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing min_h') call MPI_Allreduce(local_max_h, global_max_h, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing max_h') call MPI_Allreduce(local_sum_h, global_sum_h, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing sum_h') call MPI_Allreduce(local_min_qrad, global_min_qrad, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing min_qrad') call MPI_Allreduce(local_max_qrad, global_max_qrad, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing max_qrad') call MPI_Allreduce(local_integral_qrad, global_integral_qrad, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing integral_qrad') call MPI_Allreduce(local_volume, global_volume, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing energy volume') call MPI_Allreduce(local_max_delta_T, global_max_delta_T, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing max_delta_T') call MPI_Allreduce(local_delta_h2, global_delta_h2, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing delta_h2') call MPI_Allreduce(local_h_old2, global_h_old2, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing h_old2') if (global_volume > tiny_safe) then mean_T = global_sum_T / global_volume mean_h = global_sum_h / global_volume else mean_T = zero mean_h = zero end if rel_h_residual = sqrt(global_delta_h2) / max(sqrt(global_h_old2), tiny_safe) if (flow%rank /= 0) return if (do_write_file) then filename = trim(params%output_dir)//'/diagnostics/energy_diagnostics.csv' open(newunit=unit_id, file=trim(filename), status='old', position='append', action='write') write(unit_id,'(i0,12(",",es16.8))') step, time, global_min_T, global_max_T, mean_T, & global_min_h, global_max_h, mean_h, global_min_qrad, global_max_qrad, & global_integral_qrad, global_max_delta_T, rel_h_residual close(unit_id) end if end subroutine write_energy_diagnostics_row !> Print bridge-level Cantera cache statistics, summed over flow ranks. !! !! The counters are diagnostic only. They help determine whether Cantera !! transport, thermo-sync, and species-enthalpy caches are doing useful work !! before attempting further performance changes. subroutine write_cantera_cache_stats(flow) type(flow_mpi_t), intent(in) :: flow integer, parameter :: nstats = 16 real(c_double) :: local_stats(nstats), global_stats(nstats) integer :: ierr local_stats = 0.0_c_double global_stats = 0.0_c_double call cantera_get_cache_stats_c(local_stats, nstats) call MPI_Reduce(local_stats, global_stats, nstats, MPI_DOUBLE_PRECISION, MPI_SUM, 0, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI_Reduce failed for Cantera cache statistics') if (flow%rank /= 0) return if (sum(global_stats) <= 0.0_c_double) return write(output_unit,'(a)') ' ======================================================================' write(output_unit,'(a)') ' CANTERA CACHE STATISTICS' write(output_unit,'(a)') ' Counts are summed over flow ranks; hits/misses are cell/point states.' write(output_unit,'(a)') ' ======================================================================' write(output_unit,'(a)') ' Cache Calls Points Hits Misses Hit%' write(output_unit,'(a)') ' ----------------------------------------------------------------------' call print_cache_row('Transport_mu_Dk', 1) call print_cache_row('Energy_ThermoSync', 5) call print_cache_row('SpeciesH_Bulk', 9) call print_cache_row('SpeciesH_Point', 13) write(output_unit,'(a)') ' ======================================================================' contains subroutine print_cache_row(name, offset) character(len=*), intent(in) :: name integer, intent(in) :: offset real(c_double) :: calls, points, hits, misses, hit_rate calls = global_stats(offset) points = global_stats(offset + 1) hits = global_stats(offset + 2) misses = global_stats(offset + 3) if (calls <= 0.0_c_double .and. points <= 0.0_c_double) return if (hits + misses > 0.0_c_double) then hit_rate = 100.0_c_double * hits / (hits + misses) else hit_rate = 0.0_c_double end if write(output_unit,'(2x,a24,1x,f12.0,1x,f14.0,1x,f14.0,1x,f14.0,1x,f8.2)') & trim(name), calls, points, hits, misses, hit_rate end subroutine print_cache_row end subroutine write_cantera_cache_stats end module mod_energy