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.
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 . 2. Volumetric Flux Predictor: Interpolates 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 that enforces where is the low-Mach divergence source. 5. Corrector Step: Updates cell velocities, face fluxes, and thermodynamic pressure.
| Type | Intent | Optional | 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 (u, u_star, p, phi, face_flux, mass_flux) updated in-place. |
||
| type(solver_stats_t), | intent(inout) | :: | 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