Coupling OpenFOAM and python

Hello preCICE-rs!

I am trying to run an OpenFOAM simulation using boundary values provided by python. To get an understanding of how to perform this using preCICE adapter, I tried to modify the partitioned-pipe tutorial as I wanted to avoid FEniCS installation just for the test.

This test is composed of 2 steps:

  1. Coupled simulation using the two pimpleFoam solvers:
    During this step, the solution of the Fluid2 pimpleFoam solver is stored for each iteration of the simulation using the export:csv option.
    <participant name="Fluid2">
      ...
      <export:csv every-iteration="true" directory="./fluid2-output" />
    </participant>
  1. Coupling the Fluid1 pimpleFoam solver with Fluid2 python pseudo-solver:
    A python pseudo-solver (“pseudo” because no actual calculation takes place in the python script) reads the pressure data stored in the fluid2-openfoam-pimplefoam/fluid2-output directory and supply the data to the Fluid1 pimpleFOAM solver.

Here is a graphics of the configuration for the steps 1 (top) and 2 (bottom):

Here is the python script that serves as the pseudo-solver in step 2:

from __future__ import division, print_function

import numpy as np
import os
import sys
import argparse
import glob
import pandas as po
import precice
from precice import action_write_initial_data, \
    action_read_iteration_checkpoint, action_write_iteration_checkpoint

# ---------------------------------------------------------------------------
# UTILITY FUNCTION(S)
# ---------------------------------------------------------------------------
def loadParallelCSV(fileID):
  return po.concat([po.read_csv(name, sep=";") for name in glob.glob(f"{fileID}_*.csv")],
                   ignore_index=True)
def loadSerialCSV(fileID):
  return po.read_csv(f"{fileID}.csv", sep=";")

# ---------------------------------------------------------------------------
# INPUT
# ---------------------------------------------------------------------------
print("Starting Fluid2 Solver...")

parser = argparse.ArgumentParser(description='Fluid2 pseudo-solver')
parser.add_argument('-c', '--configurationFileName', 
                    help="Name of the xml config file.", 
                    nargs='?', 
                    type=str,
                    default="precice-config.xml")
parser.add_argument('-p', '--parentDIR', 
                    help="Location where the code is launched.", 
                    nargs='?', 
                    type=str)

try:
    args = parser.parse_args()
except SystemExit:
    print("preCICE configuration file not specified")
    quit()

inputPATH = os.path.join(args.parentDIR,'..',
                         'fluid2-openfoam-pimplefoam','fluid2-output')
inputStr = os.path.join(inputPATH, 'Fluid2-Mesh-Fluid2')

# ---------------------------------------------------------------------------
# SOLVER
# ---------------------------------------------------------------------------
print("Configure preCICE...")
solver_process_index = 0
solver_process_size = 1
interface = precice.Interface("Fluid2", args.configurationFileName, 
                              solver_process_index, solver_process_size)
print("preCICE configured...")

#- Read init
pointData = loadSerialCSV(inputStr+'.init')

N = len(pointData.index)
dimensions = interface.get_dimensions()

pressure = pointData['Pressure'].to_numpy()

meshID = interface.get_mesh_id("Fluid2-Mesh")
velocityID = interface.get_data_id("Velocity", meshID)
pressureGradientID = interface.get_data_id("PressureGradient", meshID)
pressureID = interface.get_data_id("Pressure", meshID)

vertexIDs = np.zeros(N + 1)
grid = pointData[['PosX', 'PosY', 'PosZ']].to_numpy()

vertexIDs = interface.set_mesh_vertices(meshID, grid)

t = 0
it = 1

print("FLuid2: initialize preCICE...")

#- preCICE defines timestep size of solver via precice-config.xml
precice_dt = interface.initialize()

if interface.is_action_required(action_write_initial_data()):
    interface.write_block_scalar_data(pressureID, vertexIDs, pressure)
    interface.mark_action_fulfilled(action_write_initial_data())

interface.initialize_data()

if interface.is_read_data_available():
    #- NOTE: velocity and pressureGradient are imported but not used
    velocity = interface.read_block_vector_data(velocityID, vertexIDs)
    pressureGradient = interface.read_block_scalar_data(pressureGradientID, vertexIDs)

while interface.is_coupling_ongoing():
    #- When an implicit coupling scheme is used, checkpointing is required
    if interface.is_action_required(action_write_iteration_checkpoint()):
        interface.mark_action_fulfilled(action_write_iteration_checkpoint())

    pointData = loadSerialCSV(inputStr+'.it'+str(it))
    pressure = pointData['Pressure'].to_numpy()

    interface.write_block_scalar_data(pressureID, vertexIDs, pressure)
    precice_dt = interface.advance(precice_dt)
    #- NOTE: velocity and pressureGradient are imported but not used
    velocity = interface.read_block_vector_data(velocityID, vertexIDs)
    pressureGradient = interface.read_block_scalar_data(pressureGradientID, vertexIDs)

    #- Not yet converged
    if interface.is_action_required(action_read_iteration_checkpoint()):
        interface.mark_action_fulfilled(action_read_iteration_checkpoint())
    #- Converged, timestep complete
    else: 
        t += precice_dt
    it += 1

print("Exiting Fluid2-pseudoSolver")

interface.finalize()

With this setup, the number of iterations is much higher than the original pimpleFoam-pimpleFoam coupling. The python pseudo-solver terminates because it runs out of pressure data to read from the output of the previous simulation in step 1:

...
---[precice] ^[[0m iteration: 26 of 100, time-window: 5, time: 0.04 of 1, time-window-size: 0.01, max-timestep-length: 0.01, ongoing: yes, time-window-complete: no, read-iteration-checkpoint
---[precice] ^[[0m relative convergence measure: relative two-norm diff of data "Pressure" = 6.22e-05, limit = 1.00e-06, normalization = 1.45e+02, conv = false
---[precice] ^[[0m relative convergence measure: relative two-norm diff of data "Velocity" = 2.53e-08, limit = 1.00e-06, normalization = 2.61e+00, conv = true
---[precice] ^[[0m iteration: 27 of 100, time-window: 5, time: 0.04 of 1, time-window-size: 0.01, max-timestep-length: 0.01, ongoing: yes, time-window-complete: no, read-iteration-checkpoint
---[precice] ^[[0m relative convergence measure: relative two-norm diff of data "Pressure" = 6.24e-05, limit = 1.00e-06, normalization = 1.45e+02, conv = false
---[precice] ^[[0m relative convergence measure: relative two-norm diff of data "Velocity" = 1.14e-08, limit = 1.00e-06, normalization = 2.61e+00, conv = true
---[precice] ^[[0m iteration: 28 of 100, time-window: 5, time: 0.04 of 1, time-window-size: 0.01, max-timestep-length: 0.01, ongoing: yes, time-window-complete: no, read-iteration-checkpoint
Traceback (most recent call last):
...
FileNotFoundError: [Errno 2] No such file or directory: '/scratch/nkumar001/OpenFOAM/nkumar001-6/run/precice.partitioned-pipe.of-py.run20230215/fluid2-python/../fluid2-openfoam-pimplefoam/fluid2-output/Fluid2-Mesh-Fluid2.it330.csv'
---[precice] ^[[0m Implicitly finalizing in destructor

The error above basically states that the data at iteration number 330 was not found as the previous iteration saved data only until iteration number 329, at which the simulation terminated. Note that the preCICE configuration file remains the same as that provided in the tutorial for this coupling.

I have two questions in this regard:

  1. Is this a good test to check the coupling between OpenFOAM and python? I am not sure if the precision of the exported data affects the convergence of the OpenFOAM solver?
  2. Is there something I need to modify in the interface.is_coupling_ongoing loop in the python script to make this coupling work?

I understand that the coupling of OpenFOAM and python has previously been attempted, as hinted here: Using preCICE to couple python and OpenFOAM - #3 by Makis. Any input in this regard will be very helpful.

Thanks in advance.

It sounds like you are trying to recreate something like the (new) replay mode of ASTE: Artificial Solver Testing Environment (ASTE) | preCICE - The Coupling Library

With this, you can mock a preCICE participant and give predefined data per time. Of course, this only works in uni-directional coupling.

It looks like in the config you have defined a bi-directional coupling (both participants read and write). But in the Python code, do you really read the data that OpenFOAM sends you? Since you don’t adapt the data you write based on what you read, it sounds reasonable that the coupling will take many iterations. In fact, I would expect that it always takes the maximum number of iterations, until it reaches a steady state.

I was wondering if there is a working example/tutorial available for this kind of problem i.e., Running an OpenFoam simulation with a boundary condition calculated by a Python script. I see plenty of tutorials using different solvers but couldn’t find an OpenFoam + Python example.

Cheers,

Depending on what you are looking for:

  • The solvers of the elastic-tube-1d tutorial are direct implementations in cpp, python and rust. They serve as raw example codes without external adapters.
  • Some tutorials have nutils participants, which are very thin python wrappers around the nutils solver.

Hello, I have created a repo with this specific use case. The case might not run out-of-the-box but you can find the relevant information on how to set up the coupling.

Mostly, you might be interested in the following scripts:
run.precice.overlap/hf-openfoam/system/preciceDict
run.precice.overlap/lf-galfree/lf-solver.py
run.precice.overlap/precice-config.xml

Hope it helps.

Ey, Nish-ant, this is great… very kind of you. Thanks, will check it out!

And fsimons, thanks as well for the suggestions,

Cheers,

Juan