advance_projection_step Subroutine

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

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.

Arguments

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

The computational mesh.
type(flow_mpi_t), intent(inout) :: flow

MPI decomposition context.
type(bc_set_t), intent(in) :: bc

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

Case parameters (dt, rho, etc).
type(transport_properties_t), intent(in) :: transport

Physical property fields.
type(flow_fields_t), intent(inout) :: fields

Flow field containers to be updated.
type(solver_stats_t), intent(inout) :: stats

Diagnostic stats to be populated.

Calls

proc~~advance_projection_step~~CallsGraph proc~advance_projection_step mod_flow_projection::advance_projection_step proc~advance_ab2 mod_flow_projection::advance_ab2 proc~advance_projection_step->proc~advance_ab2 proc~balance_neumann_outlet_flux mod_flow_projection::balance_neumann_outlet_flux proc~advance_projection_step->proc~balance_neumann_outlet_flux proc~compute_flow_diagnostics mod_flow_projection::compute_flow_diagnostics proc~advance_projection_step->proc~compute_flow_diagnostics proc~compute_flux_divergence mod_flow_projection::compute_flux_divergence proc~advance_projection_step->proc~compute_flux_divergence proc~compute_momentum_rhs mod_flow_projection::compute_momentum_rhs proc~advance_projection_step->proc~compute_momentum_rhs proc~compute_predicted_face_flux mod_flow_projection::compute_predicted_face_flux proc~advance_projection_step->proc~compute_predicted_face_flux proc~correct_cell_velocity mod_flow_projection::correct_cell_velocity proc~advance_projection_step->proc~correct_cell_velocity proc~correct_face_flux mod_flow_projection::correct_face_flux proc~advance_projection_step->proc~correct_face_flux proc~ensure_pressure_operator_cache mod_flow_projection::ensure_pressure_operator_cache proc~advance_projection_step->proc~ensure_pressure_operator_cache proc~ensure_projection_workspace mod_flow_projection::ensure_projection_workspace proc~advance_projection_step->proc~ensure_projection_workspace proc~flow_exchange_cell_matrix mod_mpi_flow::flow_exchange_cell_matrix proc~advance_projection_step->proc~flow_exchange_cell_matrix proc~flow_exchange_cell_scalar mod_mpi_flow::flow_exchange_cell_scalar proc~advance_projection_step->proc~flow_exchange_cell_scalar proc~flow_exchange_face_scalar mod_mpi_flow::flow_exchange_face_scalar proc~advance_projection_step->proc~flow_exchange_face_scalar proc~profiler_start mod_profiler::profiler_start proc~advance_projection_step->proc~profiler_start proc~profiler_stop mod_profiler::profiler_stop proc~advance_projection_step->proc~profiler_stop proc~solve_pressure_correction mod_flow_projection::solve_pressure_correction proc~advance_projection_step->proc~solve_pressure_correction proc~balance_neumann_outlet_flux->proc~profiler_start proc~balance_neumann_outlet_flux->proc~profiler_stop mpi_allreduce mpi_allreduce proc~balance_neumann_outlet_flux->mpi_allreduce proc~check_mpi mod_flow_projection::check_mpi proc~balance_neumann_outlet_flux->proc~check_mpi proc~fatal_error mod_kinds::fatal_error proc~balance_neumann_outlet_flux->proc~fatal_error proc~patch_type_for_face mod_bc::patch_type_for_face proc~balance_neumann_outlet_flux->proc~patch_type_for_face proc~compute_flow_diagnostics->mpi_allreduce proc~compute_flow_diagnostics->proc~check_mpi proc~compute_boundary_flux mod_flow_projection::compute_boundary_flux proc~compute_flow_diagnostics->proc~compute_boundary_flux proc~flow_global_max_owned mod_mpi_flow::flow_global_max_owned proc~compute_flow_diagnostics->proc~flow_global_max_owned proc~boundary_velocity mod_bc::boundary_velocity proc~compute_momentum_rhs->proc~boundary_velocity proc~face_effective_neighbor mod_bc::face_effective_neighbor proc~compute_momentum_rhs->proc~face_effective_neighbor proc~face_linear_vector mod_flow_projection::face_linear_vector proc~compute_momentum_rhs->proc~face_linear_vector proc~face_normal_distance mod_flow_projection::face_normal_distance proc~compute_momentum_rhs->proc~face_normal_distance proc~compute_momentum_rhs->proc~fatal_error proc~outward_normal mod_flow_projection::outward_normal proc~compute_momentum_rhs->proc~outward_normal proc~pressure_gradient_cell mod_flow_projection::pressure_gradient_cell proc~compute_momentum_rhs->proc~pressure_gradient_cell proc~compute_predicted_face_flux->proc~boundary_velocity proc~compute_predicted_face_flux->proc~face_effective_neighbor proc~compute_predicted_face_flux->proc~face_linear_vector proc~boundary_pressure_type mod_bc::boundary_pressure_type proc~correct_cell_velocity->proc~boundary_pressure_type proc~correct_cell_velocity->proc~face_effective_neighbor proc~face_linear_scalar mod_flow_projection::face_linear_scalar proc~correct_cell_velocity->proc~face_linear_scalar proc~correct_cell_velocity->proc~outward_normal proc~correct_face_flux->proc~boundary_pressure_type proc~correct_face_flux->proc~face_effective_neighbor proc~correct_face_flux->proc~face_normal_distance proc~ensure_pressure_operator_cache->mpi_allreduce proc~ensure_pressure_operator_cache->proc~boundary_pressure_type proc~ensure_pressure_operator_cache->proc~face_effective_neighbor proc~ensure_pressure_operator_cache->proc~face_normal_distance proc~ensure_pressure_operator_cache->proc~fatal_error proc~finalize_flow_projection_workspace mod_flow_projection::finalize_flow_projection_workspace proc~ensure_pressure_operator_cache->proc~finalize_flow_projection_workspace proc~ensure_pressure_operator_cache->proc~patch_type_for_face proc~ensure_projection_workspace->proc~finalize_flow_projection_workspace proc~flow_exchange_cell_matrix->proc~profiler_start proc~flow_exchange_cell_matrix->proc~profiler_stop mpi_irecv mpi_irecv proc~flow_exchange_cell_matrix->mpi_irecv mpi_isend mpi_isend proc~flow_exchange_cell_matrix->mpi_isend mpi_waitall mpi_waitall proc~flow_exchange_cell_matrix->mpi_waitall proc~check_mpi~2 mod_mpi_flow::check_mpi proc~flow_exchange_cell_matrix->proc~check_mpi~2 proc~flow_exchange_cell_matrix->proc~fatal_error proc~flow_exchange_cell_scalar->proc~profiler_start proc~flow_exchange_cell_scalar->proc~profiler_stop proc~flow_exchange_cell_scalar->mpi_irecv proc~flow_exchange_cell_scalar->mpi_isend proc~flow_exchange_cell_scalar->mpi_waitall proc~flow_exchange_cell_scalar->proc~check_mpi~2 proc~flow_exchange_face_scalar->proc~profiler_start proc~flow_exchange_face_scalar->proc~profiler_stop proc~flow_exchange_face_scalar->mpi_irecv proc~flow_exchange_face_scalar->mpi_isend proc~flow_exchange_face_scalar->mpi_waitall proc~flow_exchange_face_scalar->proc~check_mpi~2 mpi_wtime mpi_wtime proc~profiler_start->mpi_wtime proc~find_or_create_timer mod_profiler::find_or_create_timer proc~profiler_start->proc~find_or_create_timer proc~profiler_stop->mpi_wtime proc~profiler_stop->proc~find_or_create_timer proc~record_edge mod_profiler::record_edge proc~profiler_stop->proc~record_edge proc~solve_pressure_correction->proc~ensure_pressure_operator_cache proc~solve_pressure_correction->proc~ensure_projection_workspace proc~solve_pressure_correction->proc~flow_exchange_cell_scalar proc~flow_global_dot_owned mod_mpi_flow::flow_global_dot_owned proc~solve_pressure_correction->proc~flow_global_dot_owned proc~flow_global_two_dots_owned mod_mpi_flow::flow_global_two_dots_owned proc~solve_pressure_correction->proc~flow_global_two_dots_owned proc~pressure_matvec mod_flow_projection::pressure_matvec proc~solve_pressure_correction->proc~pressure_matvec res res proc~solve_pressure_correction->res proc~check_mpi->proc~fatal_error proc~check_mpi~2->proc~fatal_error proc~compute_boundary_flux->proc~patch_type_for_face proc~is_periodic_face mod_bc::is_periodic_face proc~face_effective_neighbor->proc~is_periodic_face proc~face_neighbor_weight mod_flow_projection::face_neighbor_weight proc~face_linear_scalar->proc~face_neighbor_weight proc~face_linear_vector->proc~face_neighbor_weight proc~face_normal_distance->proc~fatal_error proc~face_normal_distance->proc~outward_normal proc~face_normal_distance->proc~patch_type_for_face proc~flow_global_dot_owned->proc~profiler_start proc~flow_global_dot_owned->proc~profiler_stop proc~flow_global_dot_owned->mpi_allreduce proc~flow_global_dot_owned->proc~check_mpi~2 proc~flow_global_max_owned->proc~profiler_start proc~flow_global_max_owned->proc~profiler_stop proc~flow_global_max_owned->mpi_allreduce proc~flow_global_max_owned->proc~check_mpi~2 proc~flow_global_two_dots_owned->proc~profiler_start proc~flow_global_two_dots_owned->proc~profiler_stop proc~flow_global_two_dots_owned->mpi_allreduce proc~flow_global_two_dots_owned->proc~check_mpi~2 proc~pressure_gradient_cell->proc~boundary_pressure_type proc~pressure_gradient_cell->proc~face_effective_neighbor proc~pressure_gradient_cell->proc~face_linear_scalar proc~pressure_gradient_cell->proc~outward_normal proc~boundary_pressure mod_bc::boundary_pressure proc~pressure_gradient_cell->proc~boundary_pressure proc~pressure_matvec->proc~ensure_pressure_operator_cache proc~pressure_matvec->proc~profiler_start proc~pressure_matvec->proc~profiler_stop proc~face_neighbor_weight->proc~fatal_error proc~face_neighbor_weight->proc~outward_normal proc~face_neighbor_weight->proc~patch_type_for_face proc~is_periodic_face->proc~patch_type_for_face

Called by

proc~~advance_projection_step~~CalledByGraph proc~advance_projection_step mod_flow_projection::advance_projection_step program~lowmach_react_hex lowmach_react_hex program~lowmach_react_hex->proc~advance_projection_step

Source Code

   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(in) :: transport
      type(flow_fields_t), intent(inout) :: fields
      type(solver_stats_t), intent(inout) :: stats

      integer :: c
      integer, save :: projection_step_counter = 0
      logical :: do_projection_diagnostics

      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*
         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, 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, 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: b = -rho/dt * div(u*)
         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
            rhs_poisson(c) = -params%rho / params%dt * &
                             fields%div(c) * mesh%cells(c)%volume
         end do

         ! Floating pressure handle: pin first cell if no Dirichlet BC exists
         if (.not. pressure_cache%has_dirichlet_pressure) then
            rhs_poisson(1) = zero
         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, 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
            fields%p(1) = zero
         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, predicted_face_flux, fields%phi, local_face_flux)
         fields%face_flux = local_face_flux
         call flow_exchange_face_scalar(flow, fields%face_flux)
         call correct_cell_velocity(mesh, flow, bc, params, 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