correct_projection_to_lowmach_source Subroutine

public subroutine correct_projection_to_lowmach_source(mesh, flow, bc, params, transport, fields, stats)

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 via Poisson solver and updates cell velocities, face volumetric fluxes, and mass fluxes to guarantee strict low-Mach mass conservation.

Arguments

Type IntentOptional Attributes Name
type(mesh_t), intent(in) :: mesh

Computational mesh context.
type(flow_mpi_t), intent(inout) :: flow

MPI parallelization data structure.
type(bc_set_t), intent(in) :: bc

Boundary condition patches.
type(case_params_t), intent(in) :: params

Configuration input parameters.
type(transport_properties_t), intent(inout) :: transport

Physical and thermodynamic properties (e.g. density).
type(flow_fields_t), intent(inout) :: fields

Flow field arrays corrected in-place.
type(solver_stats_t), intent(inout) :: stats

Solver diagnostics statistics container.

Source Code

   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