25_pimpleSol

Important

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

0. Preface

In the previous discussion, we examined the PIMPLE algorithm in detail and briefly looked at some code snippets of the PIMPLE 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 PIMPLE 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_25_pimpleSol
cd ofsp_25_pimpleSol
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 PIMPLE 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.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 PIMPLE algorithm with the simplest code possible.

3.1. Main Source Code

Considering the discussion in 24_pimple, implement the main framework of the PIMPLE algorithm in the main source code:

ofsp_25_pimpleSol.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
54
55
56
57
58
#include "fvCFD.H"

#include "pimpleControl.H" // Header file containing the PIMPLE 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 PIMPLE algorithm."
    );

    #include "postProcess.H"

    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"

    #include "createMesh.H"

    #include "createControl.H"

    #include "createFields.H"

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

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

    while (runTime.loop()) // Time advancement
    {

        Info<< "Time = " << runTime.timeName() << nl << endl;

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop()) // PIMPLE outer loop
        {

            #include "UEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct()) // PIMPLE inner loop
            {
                #include "pEqn.H"
            }

        }

        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 the previous solvers, 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

tmp<fvVectorMatrix> tUEqn // Use OpenFOAM's temporary object class tmp to construct the equation matrix for convenient memory management
(
    fvm::ddt(U) 
  + fvm::div(phi, U)
  - fvm::laplacian(nu, U) // Diffusion term
);
fvVectorMatrix& UEqn = tUEqn.ref();

UEqn.relax();

solve(UEqn == -fvc::grad(p)); // Fully solve the momentum equation

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
volScalarField rAU(1.0/UEqn.A()); // A^-1
volVectorField HbyA(rAU*UEqn.H()); // HbyA
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA)); // phiHbyA

tmp<volScalarField> rAtU(rAU); // Temporary object for convenient memory management

fvScalarMatrix pEqn // Pressure correction equation
(
    fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
);

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

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

p.relax(); // Field relaxation, to be discussed later

// Momentum corrector

U = HbyA - rAtU*fvc::grad(p); // Momentum correction

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_25_pimpleSol;
...
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
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
70
71
72
solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
    
	"(U|k|epsilon)Final"
    {
        $U;
        relTol          0;
    }
    // Because PIMPLE has an outer loop, the final momentum predictor needs to verify velocity convergence
}

PIMPLE // Control dictionary for the PIMPLE algorithm
{
    // Basic loop control
    nOuterCorrectors         2;          // Number of outer loop iterations (PIMPLE-specific)
    nCorrectors              2;          // Number of inner loop (PISO corrector) iterations
    nNonOrthogonalCorrectors 0;          // Number of non-orthogonal correction iterations
    
    // Residual control (for early exit from outer loop)
    residualControl
    {
        p
        {
            tolerance       1e-4;        // Absolute residual threshold
            relTol          0;           // Relative residual threshold (0 means disabled)
        }
        U
        {
            tolerance       1e-4;
            relTol          0;
        }
    }
    
    // Pressure reference
    pRefCell                0;           // Mesh cell index where pressure reference is set
    pRefValue               0;           // Pressure reference value (usually 0)
}

// Relaxation factors (usually outside the PIMPLE dictionary)
relaxationFactors
{
    fields
    {
        p                   0.3;         // Relaxation factor for the pressure field
    }
    equations
    {
        U                   0.7;         // Relaxation factor for the momentum equation
        ".*"                1;           // No relaxation for other equations (matched by regular expression)
    }
}

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 PIMPLE algorithm and solver implementation.

Looking back at these algorithm discussions, we have roughly clarified the main ideas of various algorithms and their simplified code implementations. The idea of pressure-velocity coupling runs through almost all types of computational fluid dynamics solvers in OpenFOAM. Many cases are extensions based on these core algorithms.

As readers adjust the code themselves and compare it with native solvers, they will find many code statements that have been streamlined in our discussion.

On one hand, to facilitate understanding of the main algorithm ideas, we avoided involving other details. Do not worry; we will discuss and explain them all in the future. On the other hand, we used the cavity test case as a debugging case, which is relatively simple. Even without “boundary condition constraints,” “consistency checks,” etc., it can still produce usable results.

This section has completed the following discussions:

Implementation of the PIMPLE 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

24_pimple 26_epilogue