!> Fractional-step projection solver for constant- and variable-density modes. !! !! This module implements the core hydrodynamic solver for the !! LowMachReact-Hex framework. Constant-density mode uses `params%rho` and !! enforces \( \nabla \cdot u = 0 \). Guarded variable-density low-Mach mode !! uses active density `transport%rho` and enforces !! \( \nabla \cdot u = S \), where \(S\) is stored in !! `fields%projection_divergence_source` for the projection time level. !! !! 1. **Predictor Step**: An intermediate velocity \(\mathbf{u}^*\) is !! calculated by advancing the momentum equation explicitly, excluding !! the new pressure gradient. !! $$\frac{\mathbf{u}^* - \mathbf{u}^n}{\Delta t} = -(\mathbf{u} \cdot \nabla) \mathbf{u} + \nu \nabla^2 \mathbf{u} - \frac{1}{\rho} \nabla p^n + \mathbf{f}$$ !! !! 2. **Poisson Solve**: A pressure correction potential \(\phi = p^{n+1} - p^n\) !! is found by solving the Poisson equation derived from the active !! continuity constraint. !! $$\nabla^2 \phi = \frac{\rho}{\Delta t} \nabla \cdot \mathbf{u}^*$$ !! !! 3. **Corrector Step**: The final velocity \(\mathbf{u}^{n+1}\) and !! pressure \(p^{n+1}\) are updated using the potential gradient. !! $$\mathbf{u}^{n+1} = \mathbf{u}^* - \frac{\Delta t}{\rho} \nabla \phi$$ !! $$p^{n+1} = p^n + \phi$$ !! !! The linear system for the Poisson equation is solved using a !! **Preconditioned Conjugate Gradient (PCG)** method with a !! diagonal (Jacobi) preconditioner. module mod_flow_projection use mpi_f08 use mod_precision, only : rk, zero, one, half, tiny_safe, fatal_error use mod_input, only : case_params_t, max_species use mod_mesh, only : mesh_t use mod_mpi_flow, only : flow_mpi_t, flow_exchange_cell_scalar, & flow_exchange_cell_matrix, flow_exchange_face_scalar, & flow_global_dot_owned, flow_global_two_dots_owned, & flow_global_max_owned use mod_bc, only : bc_set_t, bc_periodic, bc_neumann, bc_dirichlet, & patch_type_for_face, boundary_velocity, & face_effective_neighbor, boundary_pressure_type, & boundary_pressure, boundary_species, boundary_temperature use mod_profiling, only : profiler_start, profiler_stop use mod_flow_fields, only : flow_fields_t use mod_cantera, only : transport_properties_t, evaluate_cantera_boundary_state use mod_linear_solver, only : solver_stats_t, & projection_work, & solve_pressure_correction, & subtract_owned_mean_scalar, & finalize_flow_projection_workspace, & ensure_projection_workspace, & ensure_pressure_operator_cache, & face_normal_distance, & outward_normal, & pressure_cache_has_dirichlet_pressure, & pressure_cache_has_neumann_outlet, & pressure_face_inverse_density implicit none private public :: advance_projection_step, compute_flow_diagnostics public :: compute_and_update_cfl, update_lowmach_divergence_source_from_density public :: correct_projection_to_lowmach_source public :: finalize_flow_projection_workspace, face_normal_distance, solver_stats_t contains !> Orchestrates one full fractional-step iteration. !! !! This routine implements the high-level logic of the projection method. !! It coordinates the explicit momentum update, the elliptic pressure !! solve, and the final divergence-free correction. !! !! @param mesh The computational mesh. !! @param flow MPI decomposition context. !! @param bc Boundary condition set. !! @param params Case parameters (dt, rho, etc). !! @param transport Physical property fields. !! @param fields Flow field containers to be updated. !! @param stats Diagnostic stats to be populated. !> Advances the low-Mach pressure projection system by one fractional step. !! !! Solves the fractional-step projection flow equations: !! 1. Predictor: Evaluates the momentum advection and diffusion terms using !! explicit second-order Adams-Bashforth time integration to compute !! the intermediate velocity field \( u^* \). !! 2. Volumetric Flux Predictor: Interpolates \( u^* \) to cell faces to obtain !! face-centered predicted volumetric fluxes. !! 3. Neumann Boundary Balancing: Adjusts Neumann outlet fluxes globally to !! enforce strict numerical conservation constraints and prevent Poisson drift. !! 4. Pressure Poisson Solve: Resolves the pressure potential \( \phi \) that !! enforces \( \nabla \cdot u = S_u \) where \( S_u \) is the low-Mach divergence source. !! 5. Corrector Step: Updates cell velocities, face fluxes, and thermodynamic pressure. !! !! @param mesh Computational mesh context. !! @param flow MPI parallelization data structure. !! @param bc Boundary condition patches. !! @param params Configuration input parameters. !! @param transport Physical and thermodynamic properties (e.g. density). !! @param fields Flow field arrays (u, u_star, p, phi, face_flux, mass_flux) updated in-place. !! @param stats Solver diagnostics statistics container. subroutine advance_projection_step(mesh, flow, bc, params, transport, fields, stats) 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(transport_properties_t), intent(inout) :: transport type(flow_fields_t), intent(inout) :: fields type(solver_stats_t), intent(inout) :: stats integer :: c !< Cell index counter. integer, save :: projection_step_counter = 0 !< Persistent timestep counter to schedule diagnostics. logical :: do_projection_diagnostics !< True if diagnostics/VTK output is scheduled this step. real(rk) :: source_value !< Resolved divergence source for the current cell. projection_step_counter = projection_step_counter + 1 do_projection_diagnostics = (mod(projection_step_counter, params%output_interval) == 0) .or. & (projection_step_counter == params%nsteps) call ensure_projection_workspace(mesh) call ensure_pressure_operator_cache(mesh, flow, bc) associate( & local_vec => projection_work%local_vec, & local_scalar => projection_work%local_scalar, & rhs_poisson => projection_work%rhs_poisson, & local_face_flux => projection_work%local_face_flux, & local_vec_star => projection_work%local_vec_star, & predicted_face_flux => projection_work%predicted_face_flux) ! 1. Predictor: Compute Explicit momentum RHS and Advance intermediate velocity u* ! Constant-density path has S = 0. Variable-density mode preserves the ! source prepared after the previous thermo-density sync. if (.not. params%enable_variable_density) then if (allocated(fields%divergence_source)) fields%divergence_source = zero end if ! Preserve the source time level actually consumed by this projection. ! Later energy/thermo updates may advance fields%divergence_source before ! diagnostics are written, so projection residuals compare against ! fields%projection_divergence_source. if (allocated(fields%projection_divergence_source)) then fields%projection_divergence_source = zero if (allocated(fields%divergence_source)) then fields%projection_divergence_source = fields%divergence_source end if end if call profiler_start('Projection_Momentum_RHS') call compute_momentum_rhs(mesh, flow, bc, params, transport, fields%u, fields%p, local_vec) call profiler_stop('Projection_Momentum_RHS') ! 2. Advance intermediate velocity u* locally ! This uses the PREVIOUS fields%rhs_old and the CURRENT local_vec call profiler_start('Projection_AB2') call advance_ab2(mesh, flow, params, fields, local_vec, local_vec_star) call profiler_stop('Projection_AB2') ! Store current RHS for next step (AB2) ! rhs_old is effectively partitioned; each rank only needs its owned cell values. fields%rhs_old = local_vec fields%rhs_old_valid = .true. fields%u_star = local_vec_star call profiler_start('Projection_Predict_Flux') call profiler_start('PredictFlux_Exchange_UStar') call flow_exchange_cell_matrix(flow, fields%u_star) call profiler_stop('PredictFlux_Exchange_UStar') ! 3. Interpolate predicted cell velocity to intermediate face fluxes call profiler_start('PredictFlux_Compute') call compute_predicted_face_flux(mesh, flow, bc, params, transport, fields%u_star, predicted_face_flux) call profiler_stop('PredictFlux_Compute') ! 3b. Balance mass at open boundaries to prevent drift in closed-loop systems call profiler_start('PredictFlux_Balance') call balance_neumann_outlet_flux(mesh, flow, bc, params, fields, predicted_face_flux) call profiler_stop('PredictFlux_Balance') call profiler_start('PredictFlux_Exchange_Face') call flow_exchange_face_scalar(flow, predicted_face_flux) call profiler_stop('PredictFlux_Exchange_Face') call profiler_stop('Projection_Predict_Flux') ! 4. Construct the Poisson RHS for div(u)=S. Constant-density mode ! scales by rho/dt; variable-density mode solves the variable- ! coefficient projection in flux form with the source volume integral. call profiler_start('Projection_Poisson_RHS') call compute_flux_divergence(mesh, flow, predicted_face_flux, local_scalar) fields%div = local_scalar rhs_poisson = zero do c = flow%first_cell, flow%last_cell source_value = zero if (allocated(fields%divergence_source)) then source_value = fields%divergence_source(c) end if if (params%enable_variable_density) then rhs_poisson(c) = -(fields%div(c) - source_value) * & mesh%cells(c)%volume / params%dt else rhs_poisson(c) = -params%rho / params%dt * & (fields%div(c) - source_value) * mesh%cells(c)%volume end if end do ! Floating pressure handle. ! Constant-density mode keeps the historical cell-1 pin. Experimental ! variable-density mode keeps every continuity equation active and ! removes the pure-Neumann nullspace with a zero-mean RHS gauge. if (.not. pressure_cache_has_dirichlet_pressure()) then if (params%enable_variable_density) then call subtract_owned_mean_scalar(flow, rhs_poisson) else rhs_poisson(1) = zero end if end if fields%phi = zero call profiler_stop('Projection_Poisson_RHS') ! Solve Laplacian system for potential phi call profiler_start('Projection_PCG') call solve_pressure_correction(mesh, flow, bc, params, transport, rhs_poisson, fields%phi, stats) call profiler_stop('Projection_PCG') ! 5. Update Pressure: p(n+1) = p(n) + phi call profiler_start('Projection_Pressure_Update') do c = flow%first_cell, flow%last_cell fields%p(c) = fields%p(c) + fields%phi(c) end do if (.not. pressure_cache_has_dirichlet_pressure()) then if (params%enable_variable_density) then call subtract_owned_mean_scalar(flow, fields%p) else fields%p(1) = zero end if end if call flow_exchange_cell_scalar(flow, fields%p) call profiler_stop('Projection_Pressure_Update') ! 6. Correct face fluxes and cell-centered velocity locally call profiler_start('Projection_Correction') call correct_face_flux(mesh, flow, bc, params, transport, predicted_face_flux, fields%phi, local_face_flux) fields%face_flux = local_face_flux call flow_exchange_face_scalar(flow, fields%face_flux) ! Keep density-weighted face flux available for scalar/enthalpy ! transport. In constant-density mode this is exactly ! params%rho * face_flux; in variable-density mode rho_f is ! interpolated from active density. call compute_face_mass_flux(mesh, flow, bc, transport, fields%face_flux, fields%mass_flux) call flow_exchange_face_scalar(flow, fields%mass_flux) call correct_cell_velocity(mesh, flow, bc, params, transport, fields%u_star, fields%phi, local_vec) ! 7. Keep velocity ghosts current for the next predictor/species step. fields%u = local_vec call flow_exchange_cell_matrix(flow, fields%u) ! Divergence is only needed for diagnostics/output, so avoid computing ! and exchanging it on non-output steps. if (do_projection_diagnostics) then call compute_flux_divergence(mesh, flow, fields%face_flux, local_scalar) fields%div = local_scalar call flow_exchange_cell_scalar(flow, fields%div) end if call profiler_stop('Projection_Correction') if (do_projection_diagnostics) then call profiler_start('Projection_Diagnostics') call compute_flow_diagnostics(mesh, flow, bc, params, fields, stats) call profiler_stop('Projection_Diagnostics') end if end associate end subroutine advance_projection_step !> Advances velocity using the 2nd-order Adams-Bashforth scheme. !! !! Falls back to 1st-order Forward Euler on the very first timestep !! when previous RHS data is unavailable. subroutine advance_ab2(mesh, flow, params, fields, rhs, local_ustar) 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 real(rk), intent(in) :: rhs(:,:) real(rk), intent(out) :: local_ustar(:,:) integer :: c real(rk) :: c_curr, c_old local_ustar = zero if (fields%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) do c = flow%first_cell, flow%last_cell local_ustar(:, c) = fields%u(:, c) + params%dt * & (c_curr * rhs(:, c) + c_old * fields%rhs_old(:, c)) end do else do c = flow%first_cell, flow%last_cell local_ustar(:, c) = fields%u(:, c) + params%dt * rhs(:, c) end do end if associate(dummy => mesh%ncells) end associate end subroutine advance_ab2 !> Evaluates the advective, diffusive, and pressure terms of the momentum equation. !! !! Supports both **Upwind** (stable) and **Central** (high-accuracy) !! convection schemes. Advection is computed using a flux-form !! discretization. subroutine compute_momentum_rhs(mesh, flow, bc, params, transport, u, p, local_rhs) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(bc_set_t), intent(in) :: bc type(case_params_t), intent(in) :: params type(transport_properties_t), intent(in) :: transport real(rk), intent(in) :: u(:,:) real(rk), intent(in) :: p(:) real(rk), intent(out) :: local_rhs(:,:) integer :: c, lf, f, nb real(rk) :: nvec(3), nhat(3) real(rk) :: uf(3), ub(3), advected(3) real(rk) :: un_area real(rk) :: conv(3), diff(3), gradp(3) real(rk) :: dist, diff_face, rho_cell logical :: use_central use_central = .false. select case (trim(params%momentum_convection_scheme)) case ('central', 'central_difference', 'central-difference') use_central = .true. case ('upwind', 'first_order_upwind', 'first-order-upwind') use_central = .false. case default call fatal_error('flow', 'unknown momentum_convection_scheme '//trim(params%momentum_convection_scheme)) end select local_rhs = zero do c = flow%first_cell, flow%last_cell conv = zero diff = zero ! Compute local pressure gradient at cell center call pressure_gradient_cell(mesh, bc, p, c, gradp) do lf = 1, mesh%ncell_faces(c) f = mesh%cell_faces(lf, c) nvec = outward_normal(mesh, f, c) nhat = nvec nb = face_effective_neighbor(mesh, bc, f, c) if (nb > 0) then uf = face_linear_vector(mesh, bc, f, c, nb, u(:, c), u(:, nb)) ub = u(:, nb) dist = face_normal_distance(mesh, bc, f, c, nb) else call boundary_velocity(mesh, bc, f, u(:, c), ub, & boundary_density=boundary_density_for_cell(params, transport, c)) uf = ub dist = face_normal_distance(mesh, bc, f, c, 0) end if un_area = dot_product(uf, nhat) * mesh%faces(f)%area ! Scheme selection for advection if (use_central) then advected = uf else if (un_area >= zero) then advected = u(:, c) else advected = ub end if end if conv = conv + un_area * advected ! Viscous diffusion Term if (nb > 0) then diff_face = half * (transport%nu(c) + transport%nu(nb)) else diff_face = transport%nu(c) end if diff = diff + diff_face * mesh%faces(f)%area * (ub - u(:, c)) / dist end do ! Assemble total RHS. Constant-density mode uses params%rho; guarded ! variable-density mode uses the active cell density. rho_cell = params%rho if (params%enable_variable_density) rho_cell = max(transport%rho(c), tiny_safe) local_rhs(:, c) = -conv / mesh%cells(c)%volume + & diff / mesh%cells(c)%volume + & params%body_force - gradp / rho_cell end do end subroutine compute_momentum_rhs real(rk) function boundary_density_for_cell(params, transport, cell_id) result(rho_b) type(case_params_t), intent(in) :: params type(transport_properties_t), intent(in) :: transport integer, intent(in) :: cell_id rho_b = params%rho if (params%enable_variable_density) then if (allocated(transport%rho)) then if (cell_id >= 1 .and. cell_id <= size(transport%rho)) then rho_b = max(transport%rho(cell_id), tiny_safe) end if end if end if end function boundary_density_for_cell !> Calculates the pressure gradient at a cell center using Gauss's Theorem. !! !! $$\nabla p = \frac{1}{V_c} \int_S p \mathbf{n} dS$$ subroutine pressure_gradient_cell(mesh, bc, p, cell_id, gradp) type(mesh_t), intent(in) :: mesh type(bc_set_t), intent(in) :: bc real(rk), intent(in) :: p(:) integer, intent(in) :: cell_id real(rk), intent(out) :: gradp(3) integer :: lf, f, nb, patch_id real(rk) :: nvec(3) real(rk) :: pf, dist logical :: is_dirichlet gradp = zero do lf = 1, mesh%ncell_faces(cell_id) f = mesh%cell_faces(lf, cell_id) nvec = outward_normal(mesh, f, cell_id) nb = face_effective_neighbor(mesh, bc, f, cell_id) if (nb > 0) then pf = face_linear_scalar(mesh, bc, f, cell_id, nb, p(cell_id), p(nb)) else if (boundary_pressure_type(mesh, bc, f) == bc_dirichlet) then call boundary_pressure(mesh, bc, f, p(cell_id), pf, is_dirichlet) else patch_id = mesh%faces(f)%patch if (patch_id > 0) then dist = face_normal_distance(mesh, bc, f, cell_id, 0) pf = p(cell_id) + bc%patches(patch_id)%dpdn * dist else pf = p(cell_id) end if end if end if gradp = gradp + pf * nvec * mesh%faces(f)%area end do gradp = gradp / mesh%cells(cell_id)%volume end subroutine pressure_gradient_cell !> Linearly interpolates cell-centered intermediate velocity to mesh faces. subroutine compute_predicted_face_flux(mesh, flow, bc, params, transport, ustar, face_flux) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(bc_set_t), intent(in) :: bc type(case_params_t), intent(in) :: params type(transport_properties_t), intent(in) :: transport real(rk), intent(in) :: ustar(:,:) real(rk), intent(out) :: face_flux(:) integer :: i, f, owner, nb real(rk) :: uf(3), ub(3) face_flux = zero do i = 1, size(flow%owned_faces) f = flow%owned_faces(i) owner = mesh%faces(f)%owner nb = face_effective_neighbor(mesh, bc, f, owner) if (nb > 0) then uf = face_linear_vector(mesh, bc, f, owner, nb, ustar(:, owner), ustar(:, nb)) else call boundary_velocity(mesh, bc, f, ustar(:, owner), ub, & boundary_density=boundary_density_for_cell(params, transport, owner)) uf = ub end if face_flux(f) = dot_product(uf, mesh%faces(f)%normal) * mesh%faces(f)%area end do end subroutine compute_predicted_face_flux !> Corrects face fluxes using the pressure potential gradient. subroutine correct_face_flux(mesh, flow, bc, params, transport, predicted_flux, phi, face_flux) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(bc_set_t), intent(in) :: bc type(case_params_t), intent(in) :: params type(transport_properties_t), intent(in) :: transport real(rk), intent(in) :: predicted_flux(:) real(rk), intent(in) :: phi(:) real(rk), intent(out) :: face_flux(:) integer :: i, f, owner, nb real(rk) :: dist, beta_face face_flux = zero do i = 1, size(flow%owned_faces) f = flow%owned_faces(i) owner = mesh%faces(f)%owner nb = face_effective_neighbor(mesh, bc, f, owner) face_flux(f) = predicted_flux(f) if (nb > 0) then dist = face_normal_distance(mesh, bc, f, owner, nb) if (params%enable_variable_density) then beta_face = pressure_face_inverse_density(transport, owner, nb) face_flux(f) = predicted_flux(f) - params%dt * beta_face * & mesh%faces(f)%area * (phi(nb) - phi(owner)) / dist else face_flux(f) = predicted_flux(f) - params%dt / params%rho * & mesh%faces(f)%area * (phi(nb) - phi(owner)) / dist end if else if (boundary_pressure_type(mesh, bc, f) == bc_dirichlet) then dist = face_normal_distance(mesh, bc, f, owner, 0) if (params%enable_variable_density) then beta_face = pressure_face_inverse_density(transport, owner, 0) face_flux(f) = predicted_flux(f) - params%dt * beta_face * & mesh%faces(f)%area * (zero - phi(owner)) / dist else face_flux(f) = predicted_flux(f) - params%dt / params%rho * & mesh%faces(f)%area * (zero - phi(owner)) / dist end if end if end if end do end subroutine correct_face_flux !> Computes density-weighted face mass flux from volumetric flux. !! !! This keeps a globally valid face-centered mass flux for scalar and !! enthalpy transport. `face_flux` is volumetric flux in \( \mathrm{m^3/s} \); !! `mass_flux` is \( \rho_f \) times that flux in \( \mathrm{kg/s} \). Both !! use owner-to-neighbor orientation. !! !! In constant-density mode \(\rho_f = \rho_{\text{params}}\). In !! variable-density mode \(\rho_f\) is the arithmetic face density from !! `transport%rho` for interior faces and the owner density on physical !! boundaries. subroutine compute_face_mass_flux(mesh, flow, bc, transport, face_flux, mass_flux) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(bc_set_t), intent(in) :: bc type(transport_properties_t), intent(in) :: transport real(rk), intent(in) :: face_flux(:) real(rk), intent(out) :: mass_flux(:) integer :: i, f, owner, nb real(rk) :: rho_face if (.not. allocated(transport%rho)) then call fatal_error('flow', 'transport rho must be allocated before mass-flux update') end if mass_flux = zero do i = 1, size(flow%owned_faces) f = flow%owned_faces(i) owner = mesh%faces(f)%owner nb = face_effective_neighbor(mesh, bc, f, owner) if (nb > 0) then rho_face = 0.5_rk * (transport%rho(owner) + transport%rho(nb)) else rho_face = transport%rho(owner) end if mass_flux(f) = rho_face * face_flux(f) end do end subroutine compute_face_mass_flux !> Computes the discrete divergence of face fluxes for each cell. subroutine compute_flux_divergence(mesh, flow, face_flux, local_div) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow real(rk), intent(in) :: face_flux(:) real(rk), intent(out) :: local_div(:) integer :: c, lf, f real(rk) :: sgn local_div = zero do c = flow%first_cell, flow%last_cell do lf = 1, mesh%ncell_faces(c) f = mesh%cell_faces(lf, c) if (mesh%faces(f)%owner == c) then sgn = one else sgn = -one end if local_div(c) = local_div(c) + sgn * face_flux(f) / mesh%cells(c)%volume end do end do end subroutine compute_flux_divergence !> Updates cell-centered velocity using the pressure potential gradient. subroutine correct_cell_velocity(mesh, flow, bc, params, transport, ustar, phi, local_u) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(bc_set_t), intent(in) :: bc type(case_params_t), intent(in) :: params type(transport_properties_t), intent(in) :: transport real(rk), intent(in) :: ustar(:,:) real(rk), intent(in) :: phi(:) real(rk), intent(out) :: local_u(:,:) integer :: c, lf, f, nb real(rk) :: nvec(3), grad(3), phi_face, rho_cell local_u = zero do c = flow%first_cell, flow%last_cell grad = zero do lf = 1, mesh%ncell_faces(c) f = mesh%cell_faces(lf, c) nvec = outward_normal(mesh, f, c) nb = face_effective_neighbor(mesh, bc, f, c) if (nb > 0) then phi_face = face_linear_scalar(mesh, bc, f, c, nb, phi(c), phi(nb)) else if (boundary_pressure_type(mesh, bc, f) == bc_dirichlet) then phi_face = zero else phi_face = phi(c) end if end if grad = grad + phi_face * nvec * mesh%faces(f)%area end do grad = grad / mesh%cells(c)%volume rho_cell = params%rho if (params%enable_variable_density) rho_cell = max(transport%rho(c), tiny_safe) local_u(:, c) = ustar(:, c) - params%dt / rho_cell * grad end do end subroutine correct_cell_velocity !> Aggregates global raw-divergence, kinetic-energy, and boundary-flux data. !! !! `stats%max_divergence` and `stats%rms_divergence` are raw !! \( \nabla \cdot u \) diagnostics. They are primary projection metrics for !! constant-density runs, but variable-density validation should use !! `div(u)-S_projection` diagnostics because the target is \( \nabla \cdot u=S \). subroutine compute_flow_diagnostics(mesh, flow, bc, params, fields, stats) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(bc_set_t), intent(in) :: bc type(case_params_t), intent(in) :: params type(flow_fields_t), intent(in) :: fields type(solver_stats_t), intent(inout) :: stats real(rk) :: local_ke, global_ke real(rk) :: local_rms, global_rms real(rk) :: local_flux, global_flux real(rk) :: local_sum_vec(3), global_sum_vec(3) integer :: c, ierr stats%max_divergence = flow_global_max_owned(flow, fields%div) local_rms = zero local_ke = zero do c = flow%first_cell, flow%last_cell local_rms = local_rms + fields%div(c) * fields%div(c) local_ke = local_ke + half * params%rho * mesh%cells(c)%volume * & dot_product(fields%u(:, c), fields%u(:, c)) end do call compute_boundary_flux(mesh, flow, bc, fields%face_flux, local_flux) local_sum_vec = [local_rms, local_ke, local_flux] call MPI_Allreduce(local_sum_vec, global_sum_vec, 3, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) call check_mpi(ierr, 'diagnostic flow sums') global_rms = global_sum_vec(1) global_ke = global_sum_vec(2) global_flux = global_sum_vec(3) stats%rms_divergence = sqrt(global_rms / max(real(mesh%ncells, rk), one)) stats%kinetic_energy = global_ke stats%net_boundary_flux = global_flux end subroutine compute_flow_diagnostics !> Adjusts flux at Neumann outlets to ensure strict global mass balance. !! !! This is critical for solvers without a pressure-pinned cell, !! as numerical drift in the Poisson RHS can lead to a singular system. subroutine balance_neumann_outlet_flux(mesh, flow, bc, params, fields, face_flux) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(bc_set_t), intent(in) :: bc type(case_params_t), intent(in) :: params type(flow_fields_t), intent(in) :: fields real(rk), intent(inout) :: face_flux(:) integer :: c, f, owner, btype, ierr real(rk) :: local_net_flux real(rk) :: global_net_flux real(rk) :: local_target_flux real(rk) :: global_target_flux real(rk) :: local_outlet_area real(rk) :: global_outlet_area real(rk) :: correction_per_area real(rk) :: local_sum_vec(3), global_sum_vec(3) if (pressure_cache_has_dirichlet_pressure()) return ! If the case has no Neumann outlet faces, there is nothing to balance. ! This avoids two unnecessary MPI_Allreduce calls per timestep for closed ! domains such as lid-driven cavity. if (.not. pressure_cache_has_neumann_outlet()) return local_net_flux = zero local_target_flux = zero local_outlet_area = zero ! Low-Mach compatibility target. For incompressible/constant-density ! flow the target net volume flux is zero. For variable-density ! low-Mach flow, the boundary volume flux must equal integral(S)dV. if (params%enable_variable_density .and. allocated(fields%divergence_source)) then do c = flow%first_cell, flow%last_cell local_target_flux = local_target_flux + fields%divergence_source(c) * mesh%cells(c)%volume end do end if do f = 1, mesh%nfaces if (mesh%faces(f)%neighbor /= 0) cycle btype = patch_type_for_face(mesh, bc, f) if (btype == bc_periodic) cycle owner = mesh%faces(f)%owner if (.not. flow%owned(owner)) cycle local_net_flux = local_net_flux + face_flux(f) if (btype == bc_neumann) then local_outlet_area = local_outlet_area + mesh%faces(f)%area end if end do local_sum_vec = [local_net_flux, local_outlet_area, local_target_flux] call profiler_start('MPI_Communication') call MPI_Allreduce(local_sum_vec, global_sum_vec, 3, & MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) call check_mpi(ierr, 'outlet balance sums') call profiler_stop('MPI_Communication') global_net_flux = global_sum_vec(1) global_outlet_area = global_sum_vec(2) global_target_flux = global_sum_vec(3) if (global_outlet_area <= tiny_safe) then if (abs(global_net_flux - global_target_flux) > tiny_safe) then call fatal_error('flow', 'boundary flux incompatible with low-Mach target and no neumann outlet faces found') end if return end if correction_per_area = (global_target_flux - global_net_flux) / global_outlet_area do f = 1, mesh%nfaces if (mesh%faces(f)%neighbor /= 0) cycle btype = patch_type_for_face(mesh, bc, f) if (btype /= bc_neumann) cycle owner = mesh%faces(f)%owner if (.not. flow%owned(owner)) cycle face_flux(f) = face_flux(f) + correction_per_area * mesh%faces(f)%area end do end subroutine balance_neumann_outlet_flux !> Helper to integrated boundary flux on MPI-owned partitions. subroutine compute_boundary_flux(mesh, flow, bc, face_flux, local_flux) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(bc_set_t), intent(in) :: bc real(rk), intent(in) :: face_flux(:) real(rk), intent(out) :: local_flux integer :: f, owner, btype local_flux = zero do f = 1, mesh%nfaces if (mesh%faces(f)%neighbor /= 0) cycle btype = patch_type_for_face(mesh, bc, f) if (btype == bc_periodic) cycle owner = mesh%faces(f)%owner if (.not. flow%owned(owner)) cycle local_flux = local_flux + face_flux(f) end do end subroutine compute_boundary_flux !> Allocates and resets temporary solver vectors. !> Calculates the normal distance between cell centers or cell-to-face. !! !! Handles periodic boundary logic by accounting for the face pair offset. !! !! @param mesh The mesh. !! @param bc Boundary conditions. !! @param face_id Face index. !! @param cell_id Source cell index. !! @param nb Neighbor cell index (or 0 for boundaries). !> Computes the linear interpolation weight for a neighbor cell. function face_neighbor_weight(mesh, bc, face_id, cell_id, nb) result(w_nb) 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) :: w_nb integer :: pair_face integer :: btype real(rk) :: nvec(3) real(rk) :: d_owner real(rk) :: d_nb real(rk) :: d_total if (nb <= 0) then w_nb = zero return end if nvec = outward_normal(mesh, face_id, cell_id) d_owner = abs(dot_product(mesh%faces(face_id)%center - & mesh%cells(cell_id)%center, nvec)) 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('flow', 'periodic face has no paired face') end if d_nb = abs(dot_product(mesh%cells(nb)%center - & mesh%faces(pair_face)%center, nvec)) else d_nb = abs(dot_product(mesh%cells(nb)%center - & mesh%faces(face_id)%center, nvec)) end if else d_nb = abs(dot_product(mesh%cells(nb)%center - & mesh%faces(face_id)%center, nvec)) end if d_total = max(d_owner + d_nb, tiny_safe) w_nb = d_owner / d_total end function face_neighbor_weight !> Linearly interpolates a scalar field to a face. function face_linear_scalar(mesh, bc, face_id, cell_id, nb, owner_value, neighbor_value) result(face_value) 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), intent(in) :: owner_value, neighbor_value real(rk) :: face_value real(rk) :: w_nb w_nb = face_neighbor_weight(mesh, bc, face_id, cell_id, nb) face_value = (one - w_nb) * owner_value + w_nb * neighbor_value end function face_linear_scalar !> Linearly interpolates a vector field to a face. function face_linear_vector(mesh, bc, face_id, cell_id, nb, owner_value, neighbor_value) result(face_value) 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), intent(in) :: owner_value(3), neighbor_value(3) real(rk) :: face_value(3) real(rk) :: w_nb w_nb = face_neighbor_weight(mesh, bc, face_id, cell_id, nb) face_value = (one - w_nb) * owner_value + w_nb * neighbor_value end function face_linear_vector !> Internal helper for MPI error checking. subroutine check_mpi(ierr, where) integer, intent(in) :: ierr character(len=*), intent(in) :: where if (ierr /= MPI_SUCCESS) call fatal_error('flow', trim(where)//' MPI failure') end subroutine check_mpi !> Update the guarded variable-density low-Mach divergence source. !! !! Constant-density mode returns without changing the zero source. In !! variable-density mode the next projection targets \( \nabla \cdot u=S \) !! with the conservative source !! $$ !! S = !! \frac{\rho^{old} - \rho}{\rho \Delta t} !! - !! \frac{u \cdot \nabla \rho}{\rho}. !! $$ !! The advective term is evaluated from corrected volumetric face fluxes as !! a cell-local outward sum, !! \( u\cdot\nabla\rho \approx V_c^{-1}\sum_f F_{f,out}(\rho_f-\rho_c) \). !! Interior \(\rho_f\) values are centered; physical boundary faces use a !! zero-normal-density-gradient assumption because no explicit density !! boundary state exists yet. The routine updates `fields%divergence_source` !! for the next projection and advances `transport%rho_old` after doing so. !! !! @param mesh Computational mesh context. !! @param flow MPI parallelization communicator data. !! @param bc Boundary condition set. !! @param params Global simulation configuration parameters. !! @param transport Thermodynamic property array containing density. !! @param fields Flow fields where the divergence source will be populated. !! @param species_Y Optional species mass fraction array for boundary density resolution. !! @param T_state Optional temperature field for boundary density resolution. subroutine update_lowmach_divergence_source_from_density(mesh, flow, bc, params, transport, fields, species_Y, T_state) 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(transport_properties_t), intent(inout) :: transport type(flow_fields_t), intent(inout) :: fields real(rk), intent(in), optional :: species_Y(:,:) real(rk), intent(in), optional :: T_state(:) integer :: c !< Cell index counter. integer :: lf !< Local cell-to-face loop counter. integer :: fid !< Global face index. integer :: nb !< Neighbor cell index across the face. integer :: first_cell !< Starting cell index for the MPI partition. integer :: last_cell !< Ending cell index for the MPI partition. real(rk) :: vol !< Cell volume [m^3]. real(rk) :: rho_c !< Active cell density at current time level [kg/m^3]. real(rk) :: rho_old_c !< Cell density at the previous time level [kg/m^3]. real(rk) :: rho_nb !< Neighbor cell density [kg/m^3]. real(rk) :: rho_f !< Face density resolved by upwind/interpolation [kg/m^3]. real(rk) :: flux_out !< Outward volumetric face flux [m^3/s]. real(rk) :: history_source !< Temporal history source contribution (drho/dt) / rho. real(rk) :: advective_density_term !< Sum of outward convective density flux [kg/s]. real(rk) :: advective_source !< Advective source contribution -(u . grad(rho)) / rho. real(rk) :: rho_floor !< Small positive density safety floor to prevent division-by-zero. logical :: have_face_flux !< Flag indicating whether the face volumetric flux field is allocated. ! Local boundary-evaluation variables real(rk) :: T_b !< Interpolated boundary temperature [K]. real(rk) :: Y_b(max_species) !< Interpolated boundary species mass fraction vector. real(rk) :: h_sens_b !< Boundary enthalpy placeholder. real(rk) :: rho_b !< Resolved thermodynamic density at boundary [kg/m^3]. real(rk) :: cp_b !< Boundary specific heat placeholder. real(rk) :: lambda_b !< Boundary conductivity placeholder. real(rk) :: mu_b !< Boundary dynamic viscosity placeholder. logical :: is_dirichlet_T !< True if temperature boundary is Dirichlet. logical :: is_dirichlet_Y !< True if species boundaries are Dirichlet. integer :: k !< Species index loop counter. if (.not. params%enable_variable_density) return if (.not. allocated(transport%rho)) return if (.not. allocated(fields%divergence_source)) then allocate(fields%divergence_source(mesh%ncells)) fields%divergence_source = 0.0_rk else if (size(fields%divergence_source) /= mesh%ncells) then deallocate(fields%divergence_source) allocate(fields%divergence_source(mesh%ncells)) fields%divergence_source = 0.0_rk end if if (.not. allocated(transport%rho_old)) then allocate(transport%rho_old(mesh%ncells)) transport%rho_old = transport%rho else if (size(transport%rho_old) /= mesh%ncells) then deallocate(transport%rho_old) allocate(transport%rho_old(mesh%ncells)) transport%rho_old = transport%rho end if if (size(transport%rho) < mesh%ncells) return rho_floor = tiny(1.0_rk) have_face_flux = allocated(fields%face_flux) if (have_face_flux) have_face_flux = size(fields%face_flux) >= mesh%nfaces first_cell = max(1, flow%first_cell) last_cell = min(mesh%ncells, flow%last_cell) if (last_cell < first_cell) return do c = first_cell, last_cell vol = mesh%cells(c)%volume rho_c = max(transport%rho(c), rho_floor) rho_old_c = transport%rho_old(c) if (params%dt > rho_floor) then history_source = (rho_old_c - rho_c) / (rho_c * params%dt) else history_source = 0.0_rk end if advective_density_term = 0.0_rk if (have_face_flux .and. vol > rho_floor) then do lf = 1, mesh%ncell_faces(c) fid = mesh%cell_faces(lf, c) if (mesh%faces(fid)%owner == c) then flux_out = fields%face_flux(fid) nb = mesh%faces(fid)%neighbor if (nb <= 0) nb = mesh%faces(fid)%periodic_neighbor else flux_out = -fields%face_flux(fid) nb = mesh%faces(fid)%owner end if if (nb > 0 .and. nb <= mesh%ncells) then rho_nb = max(transport%rho(nb), rho_floor) if (params%upwind_dilatation_density) then if (flux_out >= 0.0_rk) then rho_f = rho_c else rho_f = rho_nb end if else rho_f = 0.5_rk * (rho_c + rho_nb) end if else ! Physical boundary face. If variable density and Cantera are enabled, ! and species/energy are present, evaluate the boundary density using Cantera. rho_f = rho_c if (params%enable_variable_density .and. present(species_Y) .and. present(T_state)) then if (allocated(transport%rho) .and. size(T_state) >= c) then ! 1. Retrieve the boundary temperature call boundary_temperature(mesh, bc, fid, T_state(c), T_b, is_dirichlet_T) ! 2. Retrieve the boundary species mass fractions Y_b = zero do k = 1, params%nspecies call boundary_species(mesh, bc, fid, k, species_Y(k, c), Y_b(k), is_dirichlet_Y) end do ! 3. Evaluate the boundary state density using Cantera call evaluate_cantera_boundary_state(params, T_b, Y_b(1:params%nspecies), & h_sens_b, rho_b, cp_b, lambda_b, mu_b) if (params%upwind_dilatation_density) then if (flux_out >= 0.0_rk) then rho_f = rho_c else rho_f = rho_b end if else rho_f = 0.5_rk * (rho_c + rho_b) end if end if end if end if advective_density_term = advective_density_term + flux_out * (rho_f - rho_c) end do advective_density_term = advective_density_term / vol end if advective_source = -advective_density_term / rho_c fields%divergence_source(c) = history_source + advective_source ! Advance density history only after storing the source that will be ! used by the next projection. transport%rho_old(c) = transport%rho(c) end do call flow_exchange_cell_scalar(flow, fields%divergence_source) call flow_exchange_cell_scalar(flow, transport%rho_old) end subroutine update_lowmach_divergence_source_from_density !> Calculates domain-wide CFL and optionally scales the timestep size. !! !! For stability in reacting flows, the timestep growth is capped !! to 2% per step when `use_dynamic_dt` is active. subroutine compute_and_update_cfl(mesh, flow, params, fields, stats) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(case_params_t), intent(inout) :: params type(flow_fields_t), intent(in) :: fields type(solver_stats_t), intent(inout) :: stats integer :: c, lf, f real(rk) :: local_cfl_rate, max_cfl_rate real(rk) :: cell_outward_flux integer :: ierr local_cfl_rate = zero do c = flow%first_cell, flow%last_cell cell_outward_flux = zero do lf = 1, mesh%ncell_faces(c) f = mesh%cell_faces(lf, c) cell_outward_flux = cell_outward_flux + abs(fields%face_flux(f)) end do cell_outward_flux = half * cell_outward_flux / mesh%cells(c)%volume if (cell_outward_flux > local_cfl_rate) then local_cfl_rate = cell_outward_flux end if end do call MPI_Allreduce(local_cfl_rate, max_cfl_rate, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) call check_mpi(ierr, 'cfl max rate') if (params%use_dynamic_dt) then if (max_cfl_rate > tiny_safe) then params%dt = min(params%max_cfl / max_cfl_rate, & params%dt * params%dt_growth_limit, & params%max_dt) else params%dt = min(params%dt * params%dt_growth_limit, params%max_dt) end if params%dt = max(params%min_dt, min(params%dt, params%max_dt)) end if stats%cfl = max_cfl_rate * params%dt end subroutine compute_and_update_cfl !> Performs a post-chemistry corrective pressure projection step. !! !! Re-projects the velocity field onto the divergence manifold after stiff chemical !! integration has altered the local density, and thus the divergence source terms. !! This step resolves the new pressure correction potential \( \phi \) via Poisson !! solver and updates cell velocities, face volumetric fluxes, and mass fluxes to !! guarantee strict low-Mach mass conservation. !! !! @param mesh Computational mesh context. !! @param flow MPI parallelization data structure. !! @param bc Boundary condition patches. !! @param params Configuration input parameters. !! @param transport Physical and thermodynamic properties (e.g. density). !! @param fields Flow field arrays corrected in-place. !! @param stats Solver diagnostics statistics container. subroutine correct_projection_to_lowmach_source(mesh, flow, bc, params, transport, fields, stats) 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(transport_properties_t), intent(inout) :: transport type(flow_fields_t), intent(inout) :: fields type(solver_stats_t), intent(inout) :: stats integer :: c !< Cell index loop counter. real(rk) :: source_value !< Local divergence source term. if (.not. params%enable_density_corrector_projection) return if (.not. params%enable_variable_density) return call ensure_projection_workspace(mesh) call ensure_pressure_operator_cache(mesh, flow, bc) associate(rhs_poisson => projection_work%rhs_poisson, & local_scalar => projection_work%local_scalar, & local_face_flux => projection_work%local_face_flux, & local_vec => projection_work%local_vec) ! Compute current face-flux divergence using the live face_flux. call compute_flux_divergence(mesh, flow, fields%face_flux, local_scalar) ! Build Poisson RHS for the residual div(u) - S. rhs_poisson = zero do c = flow%first_cell, flow%last_cell source_value = zero if (allocated(fields%divergence_source)) then source_value = fields%divergence_source(c) end if if (params%enable_variable_density) then rhs_poisson(c) = -(local_scalar(c) - source_value) * & mesh%cells(c)%volume / params%dt else rhs_poisson(c) = -params%rho / params%dt * & (local_scalar(c) - source_value) * mesh%cells(c)%volume end if end do ! Remove pure-Neumann nullspace (zero-mean gauge for VD mode). if (.not. pressure_cache_has_dirichlet_pressure()) then if (params%enable_variable_density) then call subtract_owned_mean_scalar(flow, rhs_poisson) else rhs_poisson(1) = zero end if end if fields%phi = zero call solve_pressure_correction(mesh, flow, bc, params, transport, rhs_poisson, fields%phi, stats) ! Pressure update. do c = flow%first_cell, flow%last_cell fields%p(c) = fields%p(c) + fields%phi(c) end do if (.not. pressure_cache_has_dirichlet_pressure()) then if (params%enable_variable_density) then call subtract_owned_mean_scalar(flow, fields%p) else fields%p(1) = zero end if end if call flow_exchange_cell_scalar(flow, fields%p) ! Correct face fluxes. call correct_face_flux(mesh, flow, bc, params, transport, fields%face_flux, fields%phi, local_face_flux) fields%face_flux = local_face_flux call flow_exchange_face_scalar(flow, fields%face_flux) call compute_face_mass_flux(mesh, flow, bc, transport, fields%face_flux, fields%mass_flux) call flow_exchange_face_scalar(flow, fields%mass_flux) ! Correct cell-centred velocity. call correct_cell_velocity(mesh, flow, bc, params, transport, fields%u, fields%phi, local_vec) fields%u = local_vec call flow_exchange_cell_matrix(flow, fields%u) end associate end subroutine correct_projection_to_lowmach_source end module mod_flow_projection