23_pisoSol

Important

Visit https://aerosand.cc for the latest updates.

0. Preface

In the previous discussion, we examined the PISO algorithm in detail and briefly looked at some code snippets of the PISO algorithm framework in OpenFOAM.

Here, we will implement the code based on our discussion of the algorithm and understand the usage of some code.

This section primarily discusses:

Implementation of the PISO algorithm
Computing the cavity test case

1. Governing Equations

The governing equations are as follows:

Continuity equation (mass equation):

\nabla\cdot U = 0

Momentum equation:

\frac{\partial U}{\partial t} + \nabla \cdot (UU) = \nabla\cdot(\nu\nabla U)-\nabla p

The following assumptions are still applied:

Viscous term is simplified
Gravity is neglected
Density has been accounted for

Note that the momentum equation includes a transient term.

2. Project Preparation

Run the following commands in the terminal to create the project:

terminal

ofsp
foamNewApp ofsp_23_pisoSol
cd ofsp_21_pisoSol
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity debug_case
code .

2.1. Documentation File

Provide a documentation file for the project:

README.md

## About

A solver that simply reproduces the PISO algorithm.

## Bio

aerosand @ aerosand

## Caution

Pay attention to the OpenFOAM version.

## Deploy

Ensure the solver files are complete and include the debug_case folder

Enter the following commands in the terminal:

./Allclean
./Allrun


## Event

@20260313

- Added scripts #done

.2. Script Files

We continue to use the previous scripts. Unless modified in the future, they will not be elaborated further.

Run the following command in the terminal to create the scripts:

terminal

code {Allclean,Allrun}

The cleaning script Allclean is as follows:

Allclean

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#!/bin/sh
cd "${0%/*}" || exit                                # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
#------------------------------------------------------------------------------

cd debug_case

# Function: initialize the 0 folder
init_0() {
    if [ -d "0.orig" ] && [ -d "0" ]; then
        echo "Both 0 and 0.orig exist. Removing 0 and restoring from 0.orig..."
        rm -rf 0
        cp -a 0.orig 0
    elif [ -d "0" ] && [ ! -d "0.orig" ]; then
        echo "Backing up 0 to 0.orig..."
        cp -a 0 0.orig
    elif [ -d "0.orig" ] && [ ! -d "0" ]; then
        echo "Restoring 0 from 0.orig...".
        cp -a 0.orig 0
    else
        echo "Warning: Neither 0 nor 0.orig exists. Skipping."
    fi
    echo "Aerosand >>> Initial condition done."
}

# Ensure initial conditions
init_0

# Clean the case
. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
cleanCase0

The run script Allrun is as follows:

Allrun

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
#!/bin/sh
cd "${0%/*}" || exit                                # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
#------------------------------------------------------------------------------

cd debug_case

restore0Dir

runApplication blockMesh

runApplication $(getApplication)

Do not worry about the scripts; as the discussion deepens, we will continuously expand the script content and try different writing styles.

3. Project Implementation

We will implement the PISO algorithm with the simplest code possible.

3.1. Main Source Code

Considering the discussion in 22_piso, implement the main framework of the PISO algorithm in the main source code:

ofsp_23_pisoSol.C

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include "fvCFD.H"

#include "pisoControl.H" // Header file containing the PISO algorithm control
// This header file belongs to the finiteVolume library; no additional linking is needed in Make


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Aerosand: a test solver for PISO algorithm."
    );

    #include "postProcess.H" // Refer to 13_commandLine, will not be repeated

    #include "addCheckCaseOptions.H" // Refer to 13_commandLine, will not be repeated
    #include "setRootCaseLists.H" // Refer to 13_commandLine, will not be repeated
    #include "createTime.H"
    #include "createMesh.H"
    #include "createControl.H"
    #include "createFields.H" // User-provided field inclusion


    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop()) // Time advancement
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // Pressure-velocity PISO corrector
        {
            #include "UEqn.H" // Solve the momentum predictor equation

            // --- PISO loop
            while (piso.correct()) // PISO loop, inner loop
            {
                #include "pEqn.H" // Perform multiple pressure-momentum corrections
            }
        }

        runTime.write();

        runTime.printExecutionTime(Info);
    }

    Info<< "End\n" << endl;

    return 0;
}

It can be seen that the main algorithm framework is the same as discussed in the previous section.

3.2. Field Inclusion

The field inclusion createFields.H is the same as in 21_simpleSol, as follows:

createFields.H

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
/*
    # Basic fields
    # transportProperties
    # Pressure reference
*/


// # Basic fields

Info<< "Reading field p\n" << endl;
volScalarField p
(
    IOobject
    (
        "p",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Reading field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

#include "createPhi.H"


// # transportProperties

Info<< "Reading transportProperties\n" << endl;

IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

dimensionedScalar nu
(
    "nu",
    dimViscosity,
    transportProperties
);


// # Pressure reference

label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, piso.dict(), pRefCell, pRefValue);

3.3. Momentum Predictor

The momentum predictor is in UEqn.H, with the following code:

UEqn.H

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
// Momentum predictor

  fvVectorMatrix UEqn
  (
      fvm::ddt(U) // Transient term
    + fvm::div(phi, U) // Convection term
    - fvm::laplacian(nu, U) // Diffusion term
  );
  // Pressure gradient is not included

  UEqn.relax(); // Equation relaxation

  solve(UEqn == -fvc::grad(p));

Regarding the relaxation of the momentum equation, this will be discussed later. Readers can comment out this relaxation code and compare the differences in calculation results.

3.4. Pressure-Momentum Correction

The pressure-momentum correction is in pEqn.H, with the following code:

pEqn.H

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
{
    volScalarField rAU(1.0/UEqn.A()); // A^-1 
    volVectorField HbyA(rAU*UEqn.H()); // HbyA
    surfaceScalarField phiHbyA // phiHbyA
	(
	    "phiHbyA",
	    fvc::flux(HbyA)
	);

    fvScalarMatrix pEqn
    (
        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
    );
	// Pressure correction equation
	// Similarly, fvm discretization returns a matrix, with the discretization operation targeting the parameter p, which is the unknown quantity to be solved
	// fvc discretization returns a field, which serves as a known quantity in the equation


    pEqn.setReference(pRefCell, pRefValue); // Set the pressure reference

    pEqn.solve(); // Solve the pressure correction equation

	// No pressure relaxation needed

    // Momentum corrector
    U = HbyA - rAU*fvc::grad(p); // Solve for the corrected velocity
    
    // The corrected pressure and velocity participate in the next PISO loop iteration until the specified number of iterations is reached, then proceed to the next time loop
}

Tip

Note that for ease of understanding, non-orthogonal correction is not performed here. Since the cavity test case used later has a simple mesh, this has no impact.

3.5. Project Make

As discussed above, this solver does not use additional libraries, so no extra linking specifications are needed in the project Make.

The Make/files content is as follows:

Make/files

1
2
3
ofsp_23_pisoSol.C

EXE = $(FOAM_USER_APPBIN)/ofsp_23_pisoSol

The Make/options content is as follows:

Make/options

1
2
3
4
5
6
7
EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools

3.6. Compilation

Run the following command in the terminal to compile the entire project:

terminal

wclean
wmake

Compilation is successful with no issues.

3.7. Test Case

We adjust the copied test case to test the above solver.

Modify the control dictionary controlDict as follows:

controlDict

1
2
3
4
5
6
...
application     ofsp_23_pisoSol;
...
endTime         1.0; // Extend the time to compare calculation effects with and without relaxation
// If relaxation is not performed, the calculation may diverge and terminate after some time
...

Modify the solution dictionary fvSolution as follows:

fvSolution

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }
    // Because the PISO algorithm has an inner loop, the final pressure correction needs to verify pressure convergence

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
}

PISO // Specify the dictionary parameters for the PISO algorithm
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0; // Related to non-orthogonal correction; we are not covering this part
    pRefCell        0;
    pRefValue       0;
}

Other files remain unchanged.

3.6. Compilation and Execution

Compile the project:

terminal

wclean
wmake

Run the project:

terminal

./Allclean
./Allrun

Post-process visualization:

terminal

paraFoam -case debug_case

The calculation results can be viewed in ParaView.

4. Summary

Through the implementation and discussion of the solver project, I believe we now have a relatively comprehensive understanding of the PISO algorithm and solver implementation.

In the next section, we will discuss the PIMPLE algorithm.

This section has completed the following discussions:

Implementation of the PISO algorithm
Computing the cavity test case

Support us

Tip

Hopefully, the sharing here can be helpful to you.

If you find this content helpful, your comments or donations would be greatly appreciated. Your support helps ensure the ongoing updates, corrections, refinements, and improvements to this and future series, ultimately benefiting new readers as well.

The information and message provided during donation will be displayed as an acknowledgment of your support.

Support

AliPay

Copyright @ 2026 Aerosand
Course (text, images, etc.)：CC BY-NC-SA 4.0
Code derived from OpenFOAM：GPL v3
Other code：MIT License

Last updated on May 6, 2026 • Aerosand

22_piso 24_pimple