Pressure-Poisson projection: the multigrid pressure-correction solve, the pressure update, and the velocity correction that makes Ucont divergence-free. More...

#include "variables.h"
#include "Boundaries.h"

Include dependency graph for poisson.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
PetscErrorCode	PoissonSolver_Multigrid (UserMG *usermg)
	Solves the pressure-correction equation for every block with geometric multigrid.

PetscErrorCode	AssemblePoissonOperator (UserCtx *user)
	Assembles the pressure-correction operator on one multigrid level.

PetscErrorCode	ComputePoissonRHS (UserCtx *user, Vec B)
	Forms the right-hand side of the pressure-correction equation.

PetscErrorCode	UpdatePressure (UserCtx *user)
	Adds the pressure correction to the pressure, `P += Phi`, and refreshes both fields' periodic images and ghosts.

PetscErrorCode	ProjectVelocity (UserCtx *user)
	Corrects the contravariant flux with the gradient of `Phi`.

Detailed Description

Pressure-Poisson projection: the multigrid pressure-correction solve, the pressure update, and the velocity correction that makes Ucont divergence-free.

The discrete operator, the right-hand side, and the projection gradient are built from one face-gradient definition, so the Laplacian the solver inverts is exactly the divergence of the gradient the projection applies.

Every masked operation reads the solid field from the UserCtx of the level it acts on (lNvert). Solid cells are excluded from the operator, the right-hand side, the projection, and the multigrid transfers; cells next to a solid face use one-sided transverse differences.

Definition in file poisson.h.

Function Documentation

◆ PoissonSolver_Multigrid()

PetscErrorCode PoissonSolver_Multigrid ( UserMG * usermg )

extern

Solves the pressure-correction equation for every block with geometric multigrid.

On the first call for a block, assembles the operator on every multigrid level and builds the outer Krylov solver with its PCMG preconditioner, grid transfers, level smoothers, coarse solve, and null space. The solver is kept in the finest level's UserCtx::ksp and reused on every later call; it is destroyed with the context.

Each call refreshes the ghosted contravariant flux, forms the right-hand side, solves for Phi on the finest level, and appends the iteration history to Poisson_Solver_Convergence_History_Block_<bi>.log in the run's log directory.

Solver controls are read from PETSc options under the ps_ prefix (-ps_ksp_*, -ps_mg_levels_N_*, -ps_mg_coarse_*) when the solver is built.

Parameters

[in,out] usermg Multigrid hierarchy; the finest level's Phi receives the solution.

Returns: PETSc error code. A non-finite residual or a preconditioner that could not be built stops the run with PETSC_ERR_NOT_CONVERGED.

Solves the pressure-correction equation for every block with geometric multigrid.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also: PoissonSolver_Multigrid()

Definition at line 952 of file poisson.c.

{
    SimCtx       *simCtx = usermg->mgctx[0].user[0].simCtx;
    const FieldId staggered_fields[] = {FIELD_ID_UCONT};
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting Multigrid Poisson Solve...\n");
 
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        UserCtx           *user = &usermg->mgctx[usermg->mglevels - 1].user[bi];
        KSPConvergedReason reason;
 
        if (!user->ksp) PetscCall(PoissonMultigrid_Build(usermg, bi));
 
        PetscCall(SynchronizePeriodicStaggeredFields(user, 1, staggered_fields));
        PetscCall(ComputePoissonRHS(user, user->B));
 
        PetscCall(PoissonMultigrid_OpenConvergenceLog(user->ksp, simCtx, bi));
        PetscCall(KSPSolve(user->ksp, user->B, user->Phi));
        PetscCall(PoissonMultigrid_CloseConvergenceLog(user->ksp));
 
        /* A non-finite residual or a preconditioner that could not be built leaves Phi
           meaningless, and the projection would carry it into the next momentum step,
           where it surfaces as a failure of the wrong solver. Stopping at max_it is this
           solve's normal mode and stays silent; any other divergence is reported, because
           the projection proceeds on that Phi. */
        PetscCall(KSPGetConvergedReason(user->ksp, &reason));
        PetscCheck(reason != KSP_DIVERGED_NANORINF && reason != KSP_DIVERGED_PC_FAILED,
                   PETSC_COMM_WORLD, PETSC_ERR_NOT_CONVERGED,
                   "Pressure Poisson solve on block %" PetscInt_FMT " failed at step %" PetscInt_FMT
                   " (KSP reason %s). Known causes: a multigrid hierarchy coarsened too far "
                   "(reduce poisson_solver.multigrid.levels or refine the grid), or a momentum "
                   "field that has already diverged, such as an explicit time step beyond its "
                   "stability limit.",
                   bi, simCtx->step, KSPConvergedReasons[reason]);
        if (reason < 0 && reason != KSP_DIVERGED_ITS) {
            LOG(GLOBAL, LOG_WARNING, "Pressure Poisson solve on block %" PetscInt_FMT
                " diverged at step %" PetscInt_FMT " (KSP reason %s); the projection uses the last iterate.\n",
                bi, simCtx->step, KSPConvergedReasons[reason]);
        }
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Multigrid Poisson Solve complete.\n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ AssemblePoissonOperator()

PetscErrorCode AssemblePoissonOperator ( UserCtx * user )

extern

Assembles the pressure-correction operator on one multigrid level.

Allocates user->A on the first call and reassembles its entries on later calls into the same nonzero structure. Fluid rows carry the 19-point curvilinear Laplacian; dummy rows on the domain boundary and solid rows are identities. Non-periodic faces are homogeneous Neumann; periodic faces wrap to the opposite interior layer.

Parameters

[in,out] user Level context supplying metrics, lNvert, and boundary types.

Returns: PETSc error code.

Assembles the pressure-correction operator on one multigrid level.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also: AssemblePoissonOperator()

Definition at line 251 of file poisson.c.

{
    const DMDALocalInfo info = user->info;
    const PetscInt      m[3] = {info.mx, info.my, info.mz};
    PetscBool           periodic[3];
    PoissonFaceMetrics  metrics;
    const PetscReal  ***nvert, ***aj;
    AO                  ao;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    if (!user->A) {
        PetscInt local_rows;
        PetscCall(VecGetLocalSize(user->Phi, &local_rows));
        PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, local_rows, local_rows, m[0] * m[1] * m[2], m[0] * m[1] * m[2],
                               19, NULL, 19, NULL, &user->A));
    }
    PetscCall(MatZeroEntries(user->A));
 
    PoissonOperator_PeriodicAxes(user, periodic);
    PetscCall(DMDAGetAO(user->da, &ao));
    PetscCall(PoissonOperator_GetFaceMetrics(user, &metrics));
    PetscCall(DMDAVecGetArrayRead(user->da, user->lNvert, (void *)&nvert));
    PetscCall(DMDAVecGetArrayRead(user->da, user->lAj, (void *)&aj));
 
    for (PetscInt k = info.zs; k < info.zs + info.zm; k++) {
        for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
            for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
                const PetscInt c[3] = {i, j, k};
                PetscInt       row = i + j * m[0] + k * m[0] * m[1];
                PetscInt       columns[19];
                PetscScalar    coefficients[19] = {0.0};
 
                PetscCall(AOApplicationToPetsc(ao, 1, &row));
                if (i == 0 || i == m[0] - 1 || j == 0 || j == m[1] - 1 || k == 0 || k == m[2] - 1) {
                    const PetscScalar one = 1.0;
                    PetscCall(MatSetValues(user->A, 1, &row, 1, &row, &one, INSERT_VALUES));
                    continue;
                }
 
                for (PetscInt s = 0; s < 19; s++) {
                    const PetscInt *d = POISSON_STENCIL_OFFSETS[s];
                    columns[s] = PoissonOperator_NeighborIndex(i, d[0], m[0], periodic[0]) +
                                 PoissonOperator_NeighborIndex(j, d[1], m[1], periodic[1]) * m[0] +
                                 PoissonOperator_NeighborIndex(k, d[2], m[2], periodic[2]) * m[0] * m[1];
                }
                PetscCall(AOApplicationToPetsc(ao, 19, columns));
 
                if (nvert[k][j][i] > POISSON_SOLID_THRESHOLD) {
                    /* Solid rows keep the fluid row's structure with zero couplings, so a
                       solid field that changes during a run reassembles in place. */
                    coefficients[0] = 1.0;
                    PetscCall(MatSetValues(user->A, 1, &row, 19, columns, coefficients, INSERT_VALUES));
                    continue;
                }
 
                /* Faces in the order east, west, north, south, top, bottom. Non-periodic
                   boundary faces carry no flux (homogeneous Neumann). */
                for (PetscInt n = 0; n < 3; n++) {
                    const PetscInt first = periodic[n] ? 0 : 1;
                    const PetscInt last  = periodic[n] ? m[n] - 1 : m[n] - 2;
                    for (PetscInt side = 0; side < 2; side++) {
                        PetscInt across[3] = {0, 0, 0}, own[3] = {0, 0, 0}, face_cell[3] = {i, j, k};
                        const PetscBool upper = (PetscBool)(side == 0);
 
                        across[n] = upper ? 1 : -1;
                        if (PoissonOperator_At(nvert, c, across) >= POISSON_SOLID_THRESHOLD) continue;
                        if (c[n] == (upper ? last : first)) continue;
                        if (!upper) { own[n] = -1; face_cell[n] -= 1; }
 
                        const PoissonFaceGradient face =
                            PoissonOperator_FaceGradientStencil(&metrics, nvert, face_cell, n, m, periodic);
                        PoissonOperator_AddFaceFlux(&face, own, n, upper ? 1.0 : -1.0, coefficients);
                    }
                }
 
                for (PetscInt s = 0; s < 19; s++) coefficients[s] *= -aj[k][j][i];
                PetscCall(MatSetValues(user->A, 1, &row, 19, columns, coefficients, INSERT_VALUES));
            }
        }
    }
 
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->lAj, (void *)&aj));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->lNvert, (void *)&nvert));
    PetscCall(PoissonOperator_RestoreFaceMetrics(user, &metrics));
    PetscCall(MatAssemblyBegin(user->A, MAT_FINAL_ASSEMBLY));
    PetscCall(MatAssemblyEnd(user->A, MAT_FINAL_ASSEMBLY));
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Poisson operator assembled on level %d.\n", user->thislevel);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ ComputePoissonRHS()

PetscErrorCode ComputePoissonRHS	(	UserCtx *	user,
		Vec	B
	)

extern

Forms the right-hand side of the pressure-correction equation.

Writes the scaled divergence of the ghosted contravariant flux lUcont into B, zero on dummy and solid cells, and stores its domain integral in SimCtx::poissonSourceImbalance.

Parameters

[in]	user	Finest-level context.
[out]	B	Right-hand-side vector on `user->da`.

Returns: PETSc error code.

Forms the right-hand side of the pressure-correction equation.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also: ComputePoissonRHS()

Definition at line 353 of file poisson.c.

{
    SimCtx            *simCtx = user->simCtx;
    const DMDALocalInfo info = user->info;
    const PetscInt     mx = info.mx, my = info.my, mz = info.mz;
    const PetscReal    dt = simCtx->dt;
    const Cmpnts    ***ucont;
    const PetscReal ***nvert, ***aj;
    PetscReal       ***rhs;
    PetscReal          local_sum = 0.0, global_sum = 0.0;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCall(DMDAVecGetArray(user->da, B, &rhs));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lUcont, (void *)&ucont));
    PetscCall(DMDAVecGetArrayRead(user->da, user->lNvert, (void *)&nvert));
    PetscCall(DMDAVecGetArrayRead(user->da, user->lAj, (void *)&aj));
 
    for (PetscInt k = info.zs; k < info.zs + info.zm; k++) {
        for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
            for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
                if (i == 0 || i == mx - 1 || j == 0 || j == my - 1 || k == 0 || k == mz - 1 ||
                    nvert[k][j][i] > POISSON_SOLID_THRESHOLD) {
                    rhs[k][j][i] = 0.0;
                } else {
                    rhs[k][j][i] = -(ucont[k][j][i].x - ucont[k][j][i-1].x +
                                     ucont[k][j][i].y - ucont[k][j-1][i].y +
                                     ucont[k][j][i].z - ucont[k-1][j][i].z) / dt * aj[k][j][i] * COEF_TIME_ACCURACY;
                }
            }
        }
    }
 
    /* The integral of the right-hand side is the net volume flux into the domain carried
       by the uncorrected velocity. With Neumann pressure boundaries it must vanish for
       the equation to have a solution. */
    for (PetscInt k = info.zs; k < info.zs + info.zm; k++) {
        for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
            for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
                local_sum += rhs[k][j][i] / aj[k][j][i] * dt / COEF_TIME_ACCURACY;
            }
        }
    }
    PetscCallMPI(MPI_Allreduce(&local_sum, &global_sum, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)B)));
    simCtx->poissonSourceImbalance = global_sum;
    LOG_ALLOW(GLOBAL, LOG_INFO, "Poisson source imbalance: %le\n", (double)global_sum);
 
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->lAj, (void *)&aj));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->lNvert, (void *)&nvert));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lUcont, (void *)&ucont));
    PetscCall(DMDAVecRestoreArray(user->da, B, &rhs));
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

◆ UpdatePressure()

PetscErrorCode UpdatePressure ( UserCtx * user )

extern

Adds the pressure correction to the pressure, P += Phi, and refreshes both fields' periodic images and ghosts.

Parameters

[in,out] user Block context holding P and Phi.

Returns: PETSc error code.

Adds the pressure correction to the pressure, P += Phi, and refreshes both fields' periodic images and ghosts.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also: UpdatePressure()

Definition at line 416 of file poisson.c.

{
    const FieldId periodic_fields[] = {FIELD_ID_P, FIELD_ID_PHI};
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCall(VecAXPY(user->P, 1.0, user->Phi));
    PetscCall(SynchronizePeriodicCellFields(user, 2, periodic_fields));
    PetscCall(UpdateLocalGhosts(user, FIELD_ID_P));
    PetscCall(UpdateLocalGhosts(user, FIELD_ID_PHI));
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ ProjectVelocity()

PetscErrorCode ProjectVelocity ( UserCtx * user )

extern

Corrects the contravariant flux with the gradient of Phi.

Subtracts dt / COEF_TIME_ACCURACY times the face pressure-gradient flux from every fluid face of Ucont, then refreshes the periodic images, reconstructs the Cartesian velocity, and finalizes the cell fields that depend on it.

Parameters

[in,out] user Block context; reads the ghosted lPhi.

Returns: PETSc error code.

Corrects the contravariant flux with the gradient of Phi.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also: ProjectVelocity()

Definition at line 438 of file poisson.c.

{
    SimCtx             *simCtx = user->simCtx;
    const DMDALocalInfo info = user->info;
    const PetscInt      m[3] = {info.mx, info.my, info.mz};
    const PetscInt      start[3] = {info.xs, info.ys, info.zs};
    const PetscInt      end[3] = {info.xs + info.xm, info.ys + info.ym, info.zs + info.zm};
    const PetscReal     scale = simCtx->dt / COEF_TIME_ACCURACY;
    const FieldId       staggered_fields[] = {FIELD_ID_UCONT};
    PetscBool           periodic[3];
    PetscInt            interior_start[3], interior_end[3];
    PoissonFaceMetrics  metrics;
    const PetscReal  ***nvert, ***phi;
    Cmpnts           ***ucont;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PoissonOperator_PeriodicAxes(user, periodic);
    for (PetscInt a = 0; a < 3; a++) {
        interior_start[a] = (start[a] == 0) ? 1 : start[a];
        interior_end[a]   = (end[a] == m[a]) ? m[a] - 1 : end[a];
    }
 
    PetscCall(PoissonOperator_GetFaceMetrics(user, &metrics));
    PetscCall(DMDAVecGetArrayRead(user->da, user->lNvert, (void *)&nvert));
    PetscCall(DMDAVecGetArrayRead(user->da, user->lPhi, (void *)&phi));
    PetscCall(DMDAVecGetArray(user->fda, user->Ucont, &ucont));
 
    /* One pass per face orientation. Faces between two interior cells are corrected on
       every axis; a periodic axis also corrects its seam face at index 0. */
    for (PetscInt n = 0; n < 3; n++) {
        PetscInt lo[3], hi[3];
        for (PetscInt a = 0; a < 3; a++) { lo[a] = interior_start[a]; hi[a] = interior_end[a]; }
        if (periodic[n] && start[n] == 0) lo[n] = 0;
        hi[n] = PetscMin(hi[n], periodic[n] ? m[n] - 1 : m[n] - 2);
 
        for (PetscInt k = lo[2]; k < hi[2]; k++) {
            for (PetscInt j = lo[1]; j < hi[1]; j++) {
                for (PetscInt i = lo[0]; i < hi[0]; i++) {
                    const PetscInt c[3] = {i, j, k};
                    PetscInt       across[3] = {0, 0, 0};
                    PetscReal      difference[3];
 
                    across[n] = 1;
                    if (nvert[k][j][i] > POISSON_SOLID_THRESHOLD ||
                        PoissonOperator_At(nvert, c, across) > POISSON_SOLID_THRESHOLD) continue;
 
                    const PoissonFaceGradient face =
                        PoissonOperator_FaceGradientStencil(&metrics, nvert, c, n, m, periodic);
                    for (PetscInt b = 0; b < 3; b++) {
                        if (b == n) {
                            difference[b] = PoissonOperator_At(phi, c, across) - phi[k][j][i];
                            continue;
                        }
                        const PoissonTransverseDifference diff = face.diff[b];
                        PetscInt hi_lower[3] = {0, 0, 0}, hi_upper[3] = {0, 0, 0};
                        PetscInt lo_lower[3] = {0, 0, 0}, lo_upper[3] = {0, 0, 0};
                        hi_lower[b] = diff.hi; hi_upper[b] = diff.hi; hi_upper[n] = 1;
                        lo_lower[b] = diff.lo; lo_upper[b] = diff.lo; lo_upper[n] = 1;
                        difference[b] = (diff.weight == 0.0) ? 0.0 :
                            (PoissonOperator_At(phi, c, hi_lower) + PoissonOperator_At(phi, c, hi_upper) -
                             PoissonOperator_At(phi, c, lo_lower) - PoissonOperator_At(phi, c, lo_upper)) * diff.weight;
                    }
 
                    const PetscReal flux = difference[0] * face.dot[0] * face.aj +
                                           difference[1] * face.dot[1] * face.aj +
                                           difference[2] * face.dot[2] * face.aj;
                    const PetscReal correction = flux * scale;
                    if (n == 0)      ucont[k][j][i].x -= correction;
                    else if (n == 1) ucont[k][j][i].y -= correction;
                    else             ucont[k][j][i].z -= correction;
                }
            }
        }
    }
 
    PetscCall(DMDAVecRestoreArray(user->fda, user->Ucont, &ucont));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->lPhi, (void *)&phi));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->lNvert, (void *)&nvert));
    PetscCall(PoissonOperator_RestoreFaceMetrics(user, &metrics));
 
    PetscCall(SynchronizePeriodicStaggeredFields(user, 1, staggered_fields));
    PetscCall(Contra2Cart(user));
    PetscCall(FinalizePostProjectionCellFields(user));
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Functions

Detailed Description

Function Documentation

◆ PoissonSolver_Multigrid()

◆ AssemblePoissonOperator()

◆ ComputePoissonRHS()

◆ UpdatePressure()

◆ ProjectVelocity()