foamlib provides a simple, modern, ergonomic and fast Python interface for interacting with OpenFOAM.

Parsing a volVectorField with 200k cells.¹

🚀 Introduction

foamlib is a Python package designed to simplify the manipulation of OpenFOAM cases and files. Its standalone parser makes it easy to work with OpenFOAM’s input/output files from Python, while its case-handling capabilities facilitate various execution workflows—reducing boilerplate code and enabling efficient Python-based pre- and post-processing, as well as simulation management.

Compared to PyFoam and other similar tools like fluidfoam, fluidsimfoam, and Ofpp, foamlib offers advantages such as modern Python compatibility, support for binary-formatted fields, a fully type-hinted API, and asynchronous operations; making OpenFOAM workflows more accessible and streamlined.

👋 Basics

foamlib offers the following Python classes:

• FoamFile (and FoamFieldFile): read-write access to OpenFOAM configuration and field files as if they were Python dicts, using foamlib's own parser and in-place editor. Supports ASCII and binary field formats (with or without compression).
• FoamCase: a class for configuring, running, and accessing the results of OpenFOAM cases.
• AsyncFoamCase: variant of FoamCase with asynchronous methods for running multiple cases at once.
• AsyncSlurmFoamCase: subclass of AsyncFoamCase used for running cases on a Slurm cluster.

☑️ Get started

📦 Install

• With pip:

  pip install foamlib

• With conda:

  conda install -c conda-forge foamlib

🐑 Clone a case

import os
from pathlib import Path
from foamlib import FoamCase

pitz_tutorial = FoamCase(Path(os.environ["FOAM_TUTORIALS"]) / "incompressible/simpleFoam/pitzDaily")

my_pitz = pitz_tutorial.clone("myPitz")

🏃 Run the case

my_pitz.run()

🔎 Access the results

latest_time = my_pitz[-1]

p = latest_time["p"]
U = latest_time["U"]

print(p.internal_field)
print(U.internal_field)

🧹 Clean the case

my_pitz.clean()

⚙️ Edit the controlDict file

my_pitz.control_dict["writeInterval"] = 10

📝 Make multiple file reads and writes in a single go

with my_pitz.fv_schemes as f:
    f["gradSchemes"]["default"] = f["divSchemes"]["default"]
    f["snGradSchemes"]["default"] = "uncorrected"

⏳ Run a case asynchronously

import asyncio
from foamlib import AsyncFoamCase

async def run_case():
    my_pitz_async = AsyncFoamCase(my_pitz)
    await my_pitz_async.run()

asyncio.run(run_case())

🔢 Parse a field using the FoamFieldFile class directly

from foamlib import FoamFieldFile

U = FoamFieldFile(Path(my_pitz) / "0/U")

print(U.internal_field)

🔁 Run an optimization loop on a Slurm-based cluster

import os
from pathlib import Path
from foamlib import AsyncSlurmFoamCase
from scipy.optimize import differential_evolution

base = AsyncSlurmFoamCase(Path(os.environ["FOAM_TUTORIALS"]) / "incompressible/simpleFoam/pitzDaily")

async def cost(x):
    async with base.clone() as clone:
        clone[0]["U"].boundary_field["inlet"].value = [x[0], 0, 0]
        await clone.run(fallback=True) # Run locally if Slurm is not available
        return abs(clone[-1]["U"].internal_field[0][0])

result = differential_evolution(cost, bounds=[(-1, 1)], workers=AsyncSlurmFoamCase.map, polish=False)

📄 Use it to create a run (or clean) script

#!/usr/bin/env python3
from pathlib import Path
from foamlib import FoamCase

case = FoamCase(Path(__file__).parent)
# Any additional configuration here
case.run()

▶️ A complete example

The following is a fully self-contained example that demonstrates how to create an OpenFOAM case from scratch, run it, and analyze the results.

Example

#!/usr/bin/env python3
"""Check the diffusion of a scalar field in a scalarTransportFoam case."""

import shutil
from pathlib import Path

import numpy as np
from scipy.special import erfc
from foamlib import FoamCase

path = Path(__file__).parent / "diffusionCheck"
shutil.rmtree(path, ignore_errors=True)
path.mkdir(parents=True)
(path / "system").mkdir()
(path / "constant").mkdir()
(path / "0").mkdir()

case = FoamCase(path)

with case.control_dict as f:
    f["application"] = "scalarTransportFoam"
    f["startFrom"] = "latestTime"
    f["stopAt"] = "endTime"
    f["endTime"] = 5
    f["deltaT"] = 1e-3
    f["writeControl"] = "adjustableRunTime"
    f["writeInterval"] = 1
    f["purgeWrite"] = 0
    f["writeFormat"] = "ascii"
    f["writePrecision"] = 6
    f["writeCompression"] = False
    f["timeFormat"] = "general"
    f["timePrecision"] = 6
    f["adjustTimeStep"] = False
    f["runTimeModifiable"] = False

with case.fv_schemes as f:
    f["ddtSchemes"] = {"default": "Euler"}
    f["gradSchemes"] = {"default": "Gauss linear"}
    f["divSchemes"] = {"default": "none", "div(phi,U)": "Gauss linear", "div(phi,T)": "Gauss linear"}
    f["laplacianSchemes"] = {"default": "Gauss linear corrected"}

with case.fv_solution as f:
    f["solvers"] = {"T": {"solver": "PBiCG", "preconditioner": "DILU", "tolerance": 1e-6, "relTol": 0}}

with case.block_mesh_dict as f:
    f["scale"] = 1
    f["vertices"] = [
        [0, 0, 0],
        [1, 0, 0],
        [1, 0.5, 0],
        [1, 1, 0],
        [0, 1, 0],
        [0, 0.5, 0],
        [0, 0, 0.1],
        [1, 0, 0.1],
        [1, 0.5, 0.1],
        [1, 1, 0.1],
        [0, 1, 0.1],
        [0, 0.5, 0.1],
    ]
    f["blocks"] = [
        "hex", [0, 1, 2, 5, 6, 7, 8, 11], [400, 20, 1], "simpleGrading", [1, 1, 1],
        "hex", [5, 2, 3, 4, 11, 8, 9, 10], [400, 20, 1], "simpleGrading", [1, 1, 1],
    ]
    f["edges"] = []
    f["boundary"] = [
        ("inletUp", {"type": "patch", "faces": [[5, 4, 10, 11]]}),
        ("inletDown", {"type": "patch", "faces": [[0, 5, 11, 6]]}),
        ("outletUp", {"type": "patch", "faces": [[2, 3, 9, 8]]}),
        ("outletDown", {"type": "patch", "faces": [[1, 2, 8, 7]]}),
        ("walls", {"type": "wall", "faces": [[4, 3, 9, 10], [0, 1, 7, 6]]}),
        ("frontAndBack", {"type": "empty", "faces": [[0, 1, 2, 5], [5, 2, 3, 4], [6, 7, 8, 11], [11, 8, 9, 10]]}),
    ]
    f["mergePatchPairs"] = []

with case.transport_properties as f:
    f["DT"] = f.Dimensioned(1e-3, f.DimensionSet(length=2, time=-1), "DT")

with case[0]["U"] as f:
    f.dimensions = f.DimensionSet(length=1, time=-1)
    f.internal_field = [1, 0, 0]
    f.boundary_field = {
        "inletUp": {"type": "fixedValue", "value": [1, 0, 0]},
        "inletDown": {"type": "fixedValue", "value": [1, 0, 0]},
        "outletUp": {"type": "zeroGradient"},
        "outletDown": {"type": "zeroGradient"},
        "walls": {"type": "zeroGradient"},
        "frontAndBack": {"type": "empty"},
    }

with case[0]["T"] as f:
    f.dimensions = f.DimensionSet(temperature=1)
    f.internal_field = 0
    f.boundary_field = {
        "inletUp": {"type": "fixedValue", "value": 0},
        "inletDown": {"type": "fixedValue", "value": 1},
        "outletUp": {"type": "zeroGradient"},
        "outletDown": {"type": "zeroGradient"},
        "walls": {"type": "zeroGradient"},
        "frontAndBack": {"type": "empty"},
    }

case.run()

x, y, z = case[0].cell_centers().internal_field.T

end = x == x.max()
x = x[end]
y = y[end]
z = z[end]

DT = case.transport_properties["DT"].value
U = case[0]["U"].internal_field[0]

for time in case[1:]:
    if U*time.time < 2*x.max():
        continue

    T = time["T"].internal_field[end]
    analytical = 0.5 * erfc((y - 0.5) / np.sqrt(4 * DT * x/U))
    if np.allclose(T, analytical, atol=0.1):
        print(f"Time {time.time}: OK")
    else:
        raise RuntimeError(f"Time {time.time}: {T} != {analytical}")

📘 API documentation

For more information on how to use foamlibs's classes and methods, check out the documentation.

🙋 Support

If you have any questions or need help, feel free to open a discussion.

If you believe you have found a bug in foamlib, please open an issue.

🧑‍💻 Contributing

You're welcome to contribute to foamlib! Check out the contributing guidelines for more information.

Footnotes

[1] foamlib 0.8.1 vs PyFoam 2023.7 on a MacBook Air (2020, M1) with 8 GB of RAM. Benchmark script.