From OpenFOAMWiki

icoFoam solves the incompressible laminar Navier-Stokes equations using the PISO algorithm. The code is inherently transient, requiring an initial condition (such as zero velocity) and boundary conditions. The icoFOAM code can take mesh non-orthogonality into account with successive non-orthogonality iterations. The number of PISO corrections and non-orthogonality corrections are controlled through user input.

A highly annotated version of the code is shown below. Ferziger and Peric, below, give a fairly standard overview of PISO, but OpenFOAM uses a notation and form closer to Jasak and Rusche's theses. References to PISO documentation are:

  1. J.H. Ferziger and M. Peric Computational Methods for Fluid Dynamics 3rd ed. Springer 2002.
  2. Hrvoje Jasak, PhD 1996, PDF of thesis posted at: Error analysis and estimation in the Finite Volume method with applications to fluid flows.
  3. Henrik Rusche, PhD 2002 Computational fluid dynamics of dispersed two-phase flows at high phase fractions.

int main(int argc, char *argv[])
// these header files contain source code for common tasks that are used in numerous applications.
#   include "setRootCase.H"
#   include "createTime.H"
#   include "createMesh.H"
#   include "createFields.H"
#   include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    Info<< "\nStarting time loop\n" << endl;
// use the runTime object to control time stepping
    for (runTime++; !runTime.end(); runTime++)
        Info<< "Time = " << runTime.timeName() << nl << endl;
#   include "readPISOControls.H"
#   include "CourantNo.H"
//set up the linear algebra for the momentum equation.  The flux of U, phi, is treated explicity
//using the last known value of U.  
        fvVectorMatrix UEqn
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
// solve using the last known value of p on the RHS.  This gives us a velocity field that is
// not divergence free, but approximately satisfies momentum.  See Eqn. 7.31 of Ferziger & Peric
        solve(UEqn == -fvc::grad(p));
        // --- PISO loop---- take nCorr corrector steps
        for (int corr=0; corr<nCorr; corr++)
// from the last solution of velocity, extract the diag. term from the matrix and store the reciprocal
// note that the matrix coefficients are functions of U due to the non-linearity of convection.
            volScalarField rUA = 1.0/UEqn.A();
// take a Jacobi pass and update U.  See Hrv Jasak's thesis eqn. 3.137 and Henrik Rusche's thesis, eqn. 2.43
// UEqn.H is the right-hand side of the UEqn minus the product of (the off-diagonal terms and U).
// Note that since the pressure gradient is not included in the UEqn. above, this gives us U without
// the pressure gradient.  Also note that UEqn.H() is a function of U.
            U = rUA*UEqn.H();
// calculate the fluxes by dotting the interpolated velocity (to cell faces) with face normals
// The ddtPhiCorr term accounts for the divergence of the face velocity field by taking out the 
//  difference between the interpolated velocity and the flux.
            phi = (fvc::interpolate(U) & mesh.Sf()) 
                + fvc::ddtPhiCorr(rUA, U, phi);  
// adjusts the inlet and outlet fluxes to obey continuity, which is necessary for creating a well-posed
// problem where a solution for pressure exists.
            adjustPhi(phi, U, p);  
// iteratively correct for non-orthogonality.  The non-orthogonal part of the Laplacian is calculated from the most recent
// solution for pressure, using a deferred-correction approach.
            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
// set up the pressure equation
                fvScalarMatrix pEqn
                    fvm::laplacian(rUA, p) == fvc::div(phi)
// in incompressible flow, only relative pressure matters.  Unless there is a pressure BC present,
// one cell's pressure can be set arbitrarily to produce a unique pressure solution
                pEqn.setReference(pRefCell, pRefValue);
// on the last non-orthogonality correction, correct the flux using the most up-to-date pressure
                if (nonOrth == nNonOrthCorr)
//The .flux method includes contributions from all implicit terms of the pEqn (the Laplacian)
                    phi -= pEqn.flux();  
            }  // end of non-orthogonality looping
#           include "continuityErrs.H"
// add pressure gradient to interior velocity and BC's.  Note that this pressure is not just a small
// correction to a previous pressure, but is the entire pressure field.  Contrast this to the use of p'
// in Ferziger & Peric, Eqn. 7.37.
            U -= rUA*fvc::grad(p);
        } // end of the PISO loop
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    } // end of the time step loop
    Info<< "End\n" << endl;
// ************************************************************************* //