icoFlux

Here is a pythonFlu version of referenced icoFoam C++ solver

from Foam.OpenFOAM.include import setRootCase
args = setRootCase( argc, argv )
 
from Foam.OpenFOAM.include import createTime
runTime = createTime( args )
 
from Foam.OpenFOAM.include import createMesh
mesh = createMesh( runTime )
 
transportProperties, nu, p, U, phi, pRefCell, pRefValue = createFields( runTime, mesh )
 
from Foam.finiteVolume.cfdTools.general.include import initContinuityErrs
cumulativeContErr = initContinuityErrs()
 
from Foam.OpenFOAM import ext_Info, nl
ext_Info() << "\nStarting time loop\n"
 
runTime += runTime.deltaT()
while not runTime.end() :
    ext_Info() << "Time = " <<  runTime.timeName() << nl << nl
 
    from Foam.finiteVolume.cfdTools.general.include import readPISOControls
    piso, nCorr, nNonOrthCorr, momentumPredictor, transonic, nOuterCorr = readPISOControls( mesh )
 
    from Foam.finiteVolume.cfdTools.incompressible import CourantNo
    CoNum, meanCoNum = CourantNo( mesh, phi, runTime )
 
    from Foam import fvm
    UEqn = ( fvm.ddt( U ) + fvm.div( phi, U ) - fvm.laplacian( nu, U ) )
 
    from Foam import fvc
    from Foam.finiteVolume import solve
    solve( UEqn == -fvc.grad( p ) )
 
    # --- PISO loop
 
    for corr in range( nCorr ) :
        rUA = 1.0 / UEqn.A()
 
        U.ext_assign( rUA * UEqn.H() )
        phi.ext_assign( ( fvc.interpolate( U ) & mesh.Sf() ) + fvc.ddtPhiCorr( rUA, U, phi ) )
 
        from Foam.finiteVolume import adjustPhi
        adjustPhi( phi, U, p )
 
        for nonOrth in range( nNonOrthCorr + 1 ) :
            pEqn = ( fvm.laplacian( rUA, p ) == fvc.div( phi ) )
 
            pEqn.setReference( pRefCell, pRefValue )
            pEqn.solve()
 
            if nonOrth == nNonOrthCorr:
                phi.ext_assign( phi - pEqn.flux() )
                pass
 
            pass
 
        from Foam.finiteVolume.cfdTools.incompressible import continuityErrs
        cumulativeContErr = continuityErrs( mesh, phi, runTime, cumulativeContErr )
 
        U.ext_assign( U - rUA * fvc.grad( p ) )
        U.correctBoundaryConditions();
 
        pass
 
    runTime.write()
 
    ext_Info() << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << \
                  "  ClockTime = " << runTime.elapsedClockTime() << " s" << nl << nl
 
    runTime += runTime.deltaT()
    pass
 
ext_Info() << "End\n"
 
import os
return os.EX_OK