Technical Summary: A General Approach for Running Python Codes in OpenFOAM Using an Embedded pybind11 Python Interpreter

Problem Statement

The intersection of traditional computational mechanics and machine learning (ML) has created a demand for seamless interfaces between Python-based procedures and C++-based solvers like OpenFOAM. While Python offers rapid prototyping, readability, and access to a vast ecosystem of scientific and ML tools (e.g., TensorFlow, Keras, scikit-learn), OpenFOAM is built on C++ for computational efficiency.

Existing methods to bridge this gap often require translating Python code or models into C++, which negates the primary advantage of Python: coding simplicity. Current approaches include:

• Pre/Post-processing tools: Libraries like PyFoam and fluidfoam handle case manipulation and plotting but do not integrate Python into the solver loop.
• Code Translation: Converting neural networks to C++ via tools like frugally-deep or linking directly to TensorFlow C++ APIs. These methods require complex compilation settings and source-level integration, creating a barrier to entry.
• Expression Parsers: Tools like swak4Foam allow some Python execution but often rely on lower-level interfaces or specific parsers.

There is a lack of a standard, general methodology to execute "pure" Python code directly within an OpenFOAM solver without code translation, particularly for complex tasks like boundary conditions, solver prototyping, and ML-based constitutive laws.

Methodology

The authors propose a general approach using the pybind11 library to embed a Python interpreter directly within an OpenFOAM application. This method allows OpenFOAM to call Python functions and scripts, transfer data, and retrieve results without converting the Python code to C++.

Core Implementation Steps

1. Integration: The pybind11 header-only library is installed (e.g., via Conda) and linked to the OpenFOAM build system by modifying the Make/options file to include necessary Python headers and libraries. A specific modification to setRootCase.H is required to prevent conflicts with Python modules like TensorFlow regarding command-line arguments.
2. Interpreter Initialization: The Python interpreter is initialized once within the OpenFOAM code (e.g., in a solver constructor or createFields.H) using py::initialize_interpreter().
3. Data Transfer Mechanisms:
   • Scope Object: A py::object scope is created to manage the Python namespace, allowing C++ variables to be injected into Python and vice versa.
   • Pass-by-Copy: Data is transferred by converting OpenFOAM containers (e.g., Field, List) to std::vector, then casting them to py::array (NumPy arrays) for Python. Results are cast back to C++ types. This supports both primitive types and multi-dimensional arrays.
   • Pass-by-Reference: To avoid the overhead of copying large datasets, the authors utilize ctypes in Python. The memory address and size of an OpenFOAM field are passed to Python, which creates a NumPy array view directly onto the C++ memory. Modifications in Python are immediately reflected in OpenFOAM without data duplication.
4. Execution: Python scripts are loaded at runtime using py::eval_file, and specific functions are executed via py::exec.

Key Contributions and Test Cases

The paper validates this methodology through three complementary test cases:

1. Python-based Velocity Profile Boundary Condition:

   • A custom boundary condition (pythonVelocity) was created where the velocity profile is defined by a runtime Python script.
   • The script receives face-centre coordinates and time, calculates velocities, and returns them to OpenFOAM.
   • Result: Successfully applied to the classic cavity tutorial, reproducing analytical velocity profiles.

2. Python-based Heat Transfer Prototyping Solver:

   • A wrapper for the laplacianFoam solver was developed to demonstrate rapid prototyping.
   • The governing equations (Laplace equation) were solved using a simple explicit finite difference method implemented entirely in Python.
   • Result: The Python solver produced temperature fields consistent with the native C++ laplacianFoam, demonstrating the feasibility of prototyping solution strategies in Python.

3. Field Calculations with Machine Learning (Solids4foam):

   • A mechanical constitutive law for solids4foam was implemented to calculate stress tensors from strain tensors.
   • Comparison: Seven approaches were tested:
     1. Native C++ (Reference).
     2. Analytical expression in Python.
     3. Neural Network (NN) via Keras Sequential API.
     4. NN via Keras Functional API.
     5. NN with TensorFlow eager execution disabled.
     6. NN translated to pure NumPy.
     7. NN translated to NumPy with Numba JIT compilation.
   • Data Transfer: Each approach was tested with both "cell-by-cell" transfer and "entire field" transfer, as well as "pass-by-copy" vs. "pass-by-reference."

Results and Performance Analysis

The study yielded critical insights regarding efficiency and implementation strategies:

• Accuracy: All Python-based approaches (analytical and ML) produced results indistinguishable from the native C++ reference case, with negligible errors in von Mises stress and displacement magnitude.
• Data Transfer Strategy:
  • Cell-by-Cell vs. Entire Field: Passing data as an entire field is significantly faster than cell-by-cell transfer. For the analytical Python method, the entire field approach was ~3x faster; for neural networks, the difference was orders of magnitude.
  • Copy vs. Reference: Pass-by-reference (using memory addresses) provided significant speedups over pass-by-copy. The analytical Python method using pass-by-reference achieved execution times nearly identical to the native C++ solver.
• Overhead:
  • ML Overhead: Neural networks introduced a fixed overhead (loading the interpreter and model). On a coarse mesh (1,000 cells), this overhead was substantial (>600%). On a fine mesh (400,000 cells), the overhead dropped to ~18% for Keras-based models.
  • Optimization: Pass-by-reference reduced the overhead of the best-performing neural network (Base Keras) to less than 4% on the fine mesh.
• Implementation Efficiency:
  • The Analytical-Python approach was the most efficient, adding only ~8% overhead (pass-by-copy) or ~1% (pass-by-reference) compared to C++.
  • Pure NumPy/Numba: Contrary to expectations, translating Keras models to pure NumPy or using Numba did not provide a performance advantage over the native TensorFlow/Keras implementations in this context; in some cases, they were slower.

Significance and Claims

The paper claims that this approach offers a straight-forward method to integrate Python into OpenFOAM without code translation, preserving the simplicity of Python development while leveraging OpenFOAM's robust framework.

• Feasibility: The authors demonstrate that complex tasks, including boundary conditions, solver algorithms, and ML-based constitutive laws, can be implemented in Python with minimal modification to the OpenFOAM source.
• Efficiency: The study concludes that the most computationally inexpensive approach is passing data as entire fields by reference. With this method, the overhead for feed-forward neural networks can be reduced to a few percent, making it viable for production use.
• Limitations of Current Alternatives: The paper notes that while cell-by-cell transfer is possible, it is prohibitively expensive for large fields. Furthermore, the effort to convert Python models to C++ (e.g., via frugally-deep) conflicts with the goal of rapid prototyping.
• Future Outlook: The authors state that the approach runs in parallel via standard OpenFOAM MPI (with each processor holding a local interpreter) and anticipate its utility in expanding the application of ML in computational mechanics, such as in elasto-plasticity problems.

The work positions itself not as a replacement for C++ solvers, but as a flexible bridge that allows researchers to utilize the Python ecosystem for specific solver components without sacrificing the ability to run large-scale simulations.