Field exchange: Slurm cavity¶

Apply a small seeded velocity perturbation at the start of each time step, observe U and p, and inspect the completed case.

Before running¶

Install FoamNordic, build its integration and check the OpenFOAM environment. Start inside a FoamNordic checkout containing the tutorial cases. Runs are copied below the notebook output directory.

This example uses Slurm and module load openfoam/2512 with bash. Replace the allocation account, partition and resource settings with those supported by the cluster before submitting. Run the notebook client where Slurm commands and the selected OpenFOAM module are available. The random perturbation illustrates a field transform; it does not define an equation-level closure.

Read or download the original notebook. Run the original cells below in order.

The code shows the main workload choices and leaves advanced execution options at their defaults; see the execution reference. Figure styling uses Hedgehogs helpers on Matplotlib axes.

Preamble¶

# Standard Library
from pathlib import Path

# FoamNordic
import foamnordic as fno

# Data and Numerical Computing
import numpy as np

# Visualisation
import matplotlib.pyplot as plt
from IPython.display import display
import hedgehogs as hdg

Directories¶

# Case Type and Name
CASE_TYPE = "laminar"
CASE_NAME = "lidDrivenCavity"
CATEGORY = "incompressible"

# FoamNordic Project Directory
PROJECT_DIR = Path.cwd().resolve()
BASE_DIR = next((path for path in (PROJECT_DIR, *PROJECT_DIR.parents) if (path / "tutorials" / "openfoam_tutorials").is_dir()), None)
if BASE_DIR is None:
    raise FileNotFoundError("Run this notebook from inside the FoamNordic source tree")
MAIN_DIR = BASE_DIR / "tutorials" / "foamnordic_tutorials" / CATEGORY
OF_SCRIPT_DIR = BASE_DIR / "tutorials" / "openfoam_tutorials" / CATEGORY / CASE_TYPE / CASE_NAME

# Output Directory
OUTPUT_DIR = MAIN_DIR / "output"
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

Configuration¶

# Hedgehogs configuration
FIG_X, FIG_Y = hdg.figsize()
FIGURE_SIZE = (FIG_X, FIG_Y)
PALETTE = hdg.get_palette()
hdg.set_style()

# HPC configuration
ACCOUNT = "<allocation-account>"  # e.g. project_xxxxxxx
PARTITION = "test"
TIME = "00:15:00"
N_NODES = 1
N_TASKS = 1
CPUS_PER_TASK = 1
MEM_PER_CPU = "4G"

# FoamNordic Slurm configuration
of_scheduler = fno.Slurm.openfoam(
    nodes=N_NODES,
    ntasks=N_TASKS,
    cpus_per_task=CPUS_PER_TASK,
    mem_per_cpu=MEM_PER_CPU,
)

model_scheduler = fno.Slurm.model(
    cpus_per_task=1,
    mem_per_cpu=MEM_PER_CPU
)

scheduler = fno.Slurm(
    account=ACCOUNT,
    partition=PARTITION,
    time=TIME,
    openfoam=of_scheduler,
    model=model_scheduler
)

# FoamNordic configuration
SEED = 42
key = fno.Random.key(seed=SEED, scope="global")

Example - Field Coupling¶

Case Definition¶

# Define a custom operator to perturb the velocity field
def perturb_velocity(velocity, *, key):
    key, subkey = fno.Random.split(value=key, count=2)
    scale = fno.Random.uniform(value=subkey, low=0.995, high=1.005, shape=None)
    updated_velocity = velocity * scale
    return {"updated_velocity": updated_velocity}

velocity_transform = fno.Transform(
    name="perturbVelocity",
    operator=fno.Operator.function(function=perturb_velocity),
    inputs={"velocity": fno.Field(name="U")},
    outputs={"updated_velocity": fno.Field(name="U")},
    at="time_step_start",
    key=key,
)

# Initialize the OpenFOAM case
case = fno.OpenFOAM.Case(
    name="lidDrivenCavity",
    case_dir=OF_SCRIPT_DIR,
    run_dir=OUTPUT_DIR,
    of_cmd="module load openfoam/2512",
    shell="bash",
    application="pimpleFoam",
    ranks=N_TASKS,
)

case.initialize(ranks=N_TASKS, mesh="blockMesh", validate_mesh=True);

Observation Definition¶

# Define field observations to monitor during the simulation
observations = fno.Observe(
    summaries={
        "U": ("min", "max", "mean", "l2"),
        "p": ("min", "max", "mean"),
    },
    interval=10,
)

Submit Job¶

# Connect the OpenFOAM case and execution resources
longship = fno.Longship(
    case=case,
    transforms=(velocity_transform,),
    observations=(observations,),
    scheduler=scheduler,
    verbose=True,
)

# Launch the workload
run = longship.launch(verbose=None)

# Plot observations while the solver is running
records = []
figure, axes = plt.subplots(1, 2, figsize=hdg.figsize(column="double"), layout="constrained")
hdg.plots.square(ax=axes)
velocity_line, = axes[0].plot([], [], color=PALETTE["black"])
pressure_line, = axes[1].plot([], [], color=PALETTE["black"])

axes[0].set_xlabel("Physical time")
axes[0].set_ylabel(r"$\Vert U \Vert_2$")
axes[1].set_xlabel("Physical time")
axes[1].set_ylabel(r"$\overline{p}$")

for axis in axes:
    hdg.plots.apply_grid(ax=axis)

display_handle = display(figure, display_id=True)

for record in run.observe(poll_interval=0.1, progress=True):
    records.append(record)
    physical_time = [item.time for item in records]
    velocity_l2 = [item.summary["U"].l2 for item in records]
    pressure_mean = [item.summary["p"].mean for item in records]

    velocity_line.set_data(physical_time, velocity_l2)
    pressure_line.set_data(physical_time, pressure_mean)

    for axis in axes:
        axis.relim()
        axis.autoscale_view()

    display_handle.update(figure)

plt.close(figure)

# Wait for the job to complete (polling the job status)
result = run.stop(force=False, timeout=3600, progress=True)

# Summary of the job result
result.summary(style="compact", display=True);

Post-processing¶

# Postprocessing
post = result.postprocess

velocity = post.field(name="U", time_idx=-1, physical_time=None)
pressure = post.field(name="p", time_idx=-1, physical_time=None)

print("U shape:", velocity.shape)
print("p shape:", pressure.shape)

statistics = post.statistics(
    fields=["U", "p"],
    time_idx=-1,
    physical_time=None,
    verbose=True,
)

# Visualisation of the field statistics
physical_time = np.asarray([record.time for record in records])

def extract_stat(records, field, stat):
    return np.asarray([getattr(rec.summary[field], stat) for rec in records])

u_l2 = extract_stat(records, "U", "l2")
p_mean = extract_stat(records, "p", "mean")

fig, axes = plt.subplots(1, 2, figsize=hdg.figsize(column="double"))

axes[0].plot(physical_time, u_l2, c=PALETTE["black"])
axes[0].set_xlabel("Physical time")
axes[0].set_ylabel(r"$\Vert U \Vert_2$")

axes[1].plot(physical_time, p_mean, c=PALETTE["black"])
axes[1].set_xlabel("Physical time")
axes[1].set_ylabel(r"$\overline{p}$")

for axis in axes:
    hdg.plots.apply_grid(ax=axis)

hdg.plots.square(ax=axes)
hdg.plots.show(fig=fig)

# Postprocessing with class
postprocessor = fno.Postprocess.Case(source=result)

# Load fields from the postprocessor (shape: (4096, 3) for U, (4096,) for p)
U = postprocessor.field(name="U", time_idx=-1, physical_time=None)
p = postprocessor.field(name="p", time_idx=-1, physical_time=None)

# Grid parameters for 2D lid-driven cavity (64x64 uniform grid, L = 1.0m)
n_cells = int(np.sqrt(len(U)))
L = 1.0
x = np.linspace(0, L, n_cells)
y = np.linspace(0, L, n_cells)
X, Y = np.meshgrid(x, y)

# Reshape field components to 2D structured mesh without transpose
u_x = U[:, 0].reshape(n_cells, n_cells)
u_y = U[:, 1].reshape(n_cells, n_cells)
u_mag = np.linalg.norm(U[:, :2], axis=1).reshape(n_cells, n_cells)
pressure_2d = p.reshape(n_cells, n_cells)

# Visualisation of velocity magnitude and streamlines
fig, ax = plt.subplots(figsize=hdg.figsize(column="double"))

contour = ax.contourf(X, Y, u_mag, levels=30, cmap="viridis")
cbar = fig.colorbar(contour, ax=ax)
hdg.plots.style_colorbar(cbar, label=r"$|\mathbf{U}|$ [m/s]")

ax.set_xlabel("x [m]")
ax.set_ylabel("y [m]")
ax.set_title("Velocity Field")
ax.set_aspect("equal")

hdg.plots.square(ax=ax)
hdg.plots.show(fig=fig)

Check the result¶

Inspect the run result and observation records for U and p, then load the final fields. The structured visualisation assumes the tutorial cavity grid; saved shapes and mesh dimensions should be checked against the actual case.