Compute vorticity from velocity snapshots obtained with ADIOS2 IO

dolfinx
, ,
1

Hi,

I am new to FEniCSx and FEM in general.

I am trying to compute the vorticity from velocity field data obtained by reading .bp files generated by VTK writing. The simulation in question is the 2D Navier-Stokes flow past cylinder tutorial, with a Reynolds number of 100.

The following code is how I extracted the snapshots for further analysis, where data is obtained using the read_adios2_data function.

def read_adios2_data(file_path, variable_name):
    with adios2.open(file_path, "r") as fh:
        steps = []
        for step in fh:
            var = step.read(variable_name)
            steps.append(var)
    return np.array(steps)

def extract_snapshots(data, start_time, end_time, num_snapshots, dt):
    start_idx = int(start_time / dt)
    end_idx = int(end_time / dt)
    total_steps = end_idx - start_idx + 1

    # Calculate the stride to get num_snapshots
    stride = max(1, total_steps // num_snapshots)

    # Extract snapshots with the calculated stride
    snapshots = data[start_idx : end_idx + 1 : stride]

    # If we have more than num_snapshots, truncate
    snapshots = snapshots[:num_snapshots]

    snapshots_2d = snapshots[:, :, :2]
    n_timesteps = snapshots_2d.shape[0]
    u_component = snapshots_2d[:, :, 0].reshape(n_timesteps, -1).T
    v_component = snapshots_2d[:, :, 1].reshape(n_timesteps, -1).T
    return np.vstack((u_component, v_component))

Subsequently, with the velocity components stacked (due to my use case), I will need to compute the vorticity from these snapshots, hence I assigned them back to a CG2 function space and subsequently compute the vorticity. I am not sure if this is the right way to compute vorticity from external velocity data (I may have thousands of velocity samples generated from reduced-order-modeling). I have interleaved the x,y component of velocity into u.

# Load the mesh
with io.XDMFFile(MPI.COMM_WORLD, "/path/to/mesh/mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="mesh")

domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim - 1)

# Create function spaces
v_cg2 = element("Lagrange", domain.topology.cell_name(), 2, shape=(domain.geometry.dim,))
q_dg1 = element("DG", domain.topology.cell_name(), 1)

V = functionspace(domain, v_cg2)
Q = functionspace(domain, s_cg1)

# Create functions for velocity and vorticity
u = fem.Function(V)
omega = fem.Function(Q)

# Get the number of spatial points
dofs = V.tabulate_dof_coordinates()
num_spatial_points = len(dofs)  

# Reshape to match the function space structure
u_values = np.zeros(len(dofs) * 2)

u_values[0::2] = u_snapshot[:num_spatial_points, -1]  # x-component
u_values[1::2] = u_snapshot[num_spatial_points:, -1]  # y-component

# Assign values to u
u.x.array[:] = u_values
ux, uy = u.split()
vortex = ux.dx(1)-uy.dx(0)

omega.interpolate(fem.Expression(vortex, Q.element.interpolation_points()))

# Extract vorticity values and corresponding coordinates
vorticity_values = omega.x.array
vorticity_coords = Q.tabulate_dof_coordinates()

However the vorticity I obtained is not what I expected:

aea96b04-96d1-4eb5-ba81-3a98faba0bef
aea96b04-96d1-4eb5-ba81-3a98faba0bef997×259 76.8 KB

Any help is appreciated!

2

Please make a reproducible example, as you have not shown your whole workflow (on a simple built-in mesh) so that anyone can reproduce your results.

3

Thank you for your reply! Attached is the simple channel flow problem from the (your) FEniCSx tutorial.

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista

from dolfinx.fem import (
    Constant,
    Function,
    functionspace,
    assemble_scalar,
    dirichletbc,
    form,
    locate_dofs_geometrical,
)
from dolfinx.fem.petsc import (
    assemble_matrix,
    assemble_vector,
    apply_lifting,
    create_vector,
    set_bc,
)
from dolfinx.io import VTXWriter
from dolfinx.mesh import create_unit_square
from dolfinx.plot import vtk_mesh
from basix.ufl import element
from ufl import (
    FacetNormal,
    Identity,
    TestFunction,
    TrialFunction,
    div,
    dot,
    ds,
    dx,
    inner,
    lhs,
    nabla_grad,
    rhs,
    sym,
)

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)


from dolfinx import io

# After the mesh is created
with io.XDMFFile(
    MPI.COMM_WORLD,
    "/path/to/channel_flow_mesh.xdmf",
    "w",
) as xdmf:
    xdmf.write_mesh(mesh)


t = 0
T = 10
num_steps = 500
dt = T / num_steps

v_cg2 = element("Lagrange", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim,))
s_cg1 = element("Lagrange", mesh.topology.cell_name(), 1)
V = functionspace(mesh, v_cg2)
Q = functionspace(mesh, s_cg1)

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)


def walls(x):
    return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1))


wall_dofs = locate_dofs_geometrical(V, walls)
u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bc_noslip = dirichletbc(u_noslip, wall_dofs, V)


def inflow(x):
    return np.isclose(x[0], 0)


inflow_dofs = locate_dofs_geometrical(Q, inflow)
bc_inflow = dirichletbc(PETSc.ScalarType(8), inflow_dofs, Q)


def outflow(x):
    return np.isclose(x[0], 1)


outflow_dofs = locate_dofs_geometrical(Q, outflow)
bc_outflow = dirichletbc(PETSc.ScalarType(0), outflow_dofs, Q)
bcu = [bc_noslip]
bcp = [bc_inflow, bc_outflow]

u_n = Function(V)
u_n.name = "u_n"
U = 0.5 * (u_n + u)
n = FacetNormal(mesh)
f = Constant(mesh, PETSc.ScalarType((0, 0)))
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(1))
rho = Constant(mesh, PETSc.ScalarType(1))


# Define strain-rate tensor
def epsilon(u):
    return sym(nabla_grad(u))


# Define stress tensor
def sigma(u, p):
    return 2 * mu * epsilon(u) - p * Identity(len(u))


# Define the variational problem for the first step
p_n = Function(Q)
p_n.name = "p_n"
F1 = rho * dot((u - u_n) / k, v) * dx
F1 += rho * dot(dot(u_n, nabla_grad(u_n)), v) * dx
F1 += inner(sigma(U, p_n), epsilon(v)) * dx
F1 += dot(p_n * n, v) * ds - dot(mu * nabla_grad(U) * n, v) * ds
F1 -= dot(f, v) * dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))

A1 = assemble_matrix(a1, bcs=bcu)
A1.assemble()
b1 = create_vector(L1)

# Define variational problem for step 2
u_ = Function(V)
a2 = form(dot(nabla_grad(p), nabla_grad(q)) * dx)
L2 = form(dot(nabla_grad(p_n), nabla_grad(q)) * dx - (rho / k) * div(u_) * q * dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)

# Define variational problem for step 3
p_ = Function(Q)
a3 = form(rho * dot(u, v) * dx)
L3 = form(rho * dot(u_, v) * dx - k * dot(nabla_grad(p_ - p_n), v) * dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)


# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")

# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.BCGS)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")

# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)

from pathlib import Path

folder = Path("results")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, folder / "poiseuille_u.bp", u_n, engine="BP4")
vtx_p = VTXWriter(mesh.comm, folder / "poiseuille_p.bp", p_n, engine="BP4")
vtx_u.write(t)
vtx_p.write(t)


def u_exact(x):
    values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
    values[0] = 4 * x[1] * (1.0 - x[1])
    return values


u_ex = Function(V)
u_ex.interpolate(u_exact)

L2_error = form(dot(u_ - u_ex, u_ - u_ex) * dx)

for i in range(num_steps):
    # Update current time step
    t += dt

    # Step 1: Tentative veolcity step
    with b1.localForm() as loc_1:
        loc_1.set(0)
    assemble_vector(b1, L1)
    apply_lifting(b1, [a1], [bcu])
    b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b1, bcu)
    solver1.solve(b1, u_.vector)
    u_.x.scatter_forward()

    # Step 2: Pressure corrrection step
    with b2.localForm() as loc_2:
        loc_2.set(0)
    assemble_vector(b2, L2)
    apply_lifting(b2, [a2], [bcp])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b2, bcp)
    solver2.solve(b2, p_.vector)
    p_.x.scatter_forward()

    # Step 3: Velocity correction step
    with b3.localForm() as loc_3:
        loc_3.set(0)
    assemble_vector(b3, L3)
    b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    solver3.solve(b3, u_.vector)
    u_.x.scatter_forward()
    # Update variable with solution form this time step
    u_n.x.array[:] = u_.x.array[:]
    p_n.x.array[:] = p_.x.array[:]

    # Write solutions to file
    vtx_u.write(t)
    vtx_p.write(t)

    # Compute error at current time-step
    error_L2 = np.sqrt(mesh.comm.allreduce(assemble_scalar(L2_error), op=MPI.SUM))
    error_max = mesh.comm.allreduce(
        np.max(u_.vector.array - u_ex.vector.array), op=MPI.MAX
    )
    # Print error only every 20th step and at the last step
    if (i % 20 == 0) or (i == num_steps - 1):
        print(f"Time {t:.2f}, L2-error {error_L2:.2e}, Max error {error_max:.2e}")
# Close xmdf file
vtx_u.close()
vtx_p.close()
b1.destroy()
b2.destroy()
b3.destroy()
solver1.destroy()
solver2.destroy()
solver3.destroy()

And the following is my workflow for attempting to compute the vorticity of the final timestep:

import numpy as np
import matplotlib.pyplot as plt
from dolfinx import mesh, fem, io
from mpi4py import MPI
import ufl

import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import tqdm.autonotebook

from mpi4py import MPI
from petsc4py import PETSc

from basix.ufl import element

from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx.fem import (
    Constant,
    Function,
    functionspace,
    assemble_scalar,
    dirichletbc,
    form,
    locate_dofs_topological,
    set_bc,
)
from dolfinx.fem.petsc import (
    apply_lifting,
    assemble_matrix,
    assemble_vector,
    create_vector,
    create_matrix,
    set_bc,
)
from dolfinx.graph import adjacencylist
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import VTXWriter, distribute_entity_data, gmshio
from dolfinx.mesh import create_mesh, meshtags_from_entities
from ufl import (
    FacetNormal,
    Identity,
    Measure,
    TestFunction,
    TrialFunction,
    as_vector,
    div,
    dot,
    ds,
    dx,
    inner,
    lhs,
    grad,
    nabla_grad,
    rhs,
    sym,
    system,
)

import numpy as np
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
from dolfinx import mesh, fem, io

import adios2
import hydra
from omegaconf import DictConfig, OmegaConf
from pathlib import Path


def read_adios2_data(file_path, variable_name):
    with adios2.open(file_path, "r") as fh:
        steps = []
        for step in fh:
            var = step.read(variable_name)
            steps.append(var)
    return np.array(steps)


def extract_snapshots(data, start_time, end_time, num_snapshots, dt):
    start_idx = int(start_time / dt)
    end_idx = int(end_time / dt)
    total_steps = end_idx - start_idx + 1

    # Calculate the stride to get num_snapshots
    stride = max(1, total_steps // num_snapshots)

    # Extract snapshots with the calculated stride
    snapshots = data[start_idx : end_idx + 1 : stride]

    # If we have more than num_snapshots, truncate
    snapshots = snapshots[:num_snapshots]

    snapshots_2d = snapshots[:, :, :2]
    n_timesteps = snapshots_2d.shape[0]
    u_component = snapshots_2d[:, :, 0].reshape(n_timesteps, -1).T
    v_component = snapshots_2d[:, :, 1].reshape(n_timesteps, -1).T
    return np.vstack((u_component, v_component))


# change working directory to script location
os.chdir(os.path.dirname(os.path.realpath(__file__)))

bp_file_path = Path("./results/poiseuille_u.bp")
output_base_path = Path("./results/")
variable = "u_n"
t = 0
T = 10
num_steps = 500
dt = T / num_steps
num_snapshots = 200

if not bp_file_path.exists():
    print(f"No .bp file found at {bp_file_path}. Skipping...")
    raise ValueError("No data found")

# Construct output path
output_folder = output_base_path
output_folder.mkdir(parents=True, exist_ok=True)
output_file = output_folder / f"{variable}.npy"

# Read the data
data = read_adios2_data(str(bp_file_path), variable)

snapshot_matrix = extract_snapshots(data, t, T, num_snapshots, dt)


# Load the mesh
with io.XDMFFile(
    MPI.COMM_WORLD,
    "path/to/channel_flow_mesh.xdmf",
    "r",
) as xdmf:
    mesh = xdmf.read_mesh(name="mesh")


v_cg2 = element("Lagrange", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim,))
dg1 = element("DG", mesh.topology.cell_name(), 1)
V = functionspace(mesh, v_cg2)
Q = functionspace(mesh, dg1)


# Create functions for velocity and vorticity
u = fem.Function(V)
omega = fem.Function(Q)

# Get the number of spatial points
dofs = V.tabulate_dof_coordinates()
num_spatial_points = len(dofs)

# Reshape to match the function space structure
u_values = np.zeros(len(dofs) * 2)
u_values[0::2] = snapshot_matrix[:num_spatial_points, -1]  # x-component
u_values[1::2] = snapshot_matrix[num_spatial_points:, -1]  # y-component

# Assign values to u
u.x.array[:] = u_values
ux, uy = u.split()
vortex = ux.dx(1) - uy.dx(0)


omega.interpolate(fem.Expression(vortex, Q.element.interpolation_points()))

# Extract vorticity values and corresponding coordinates
vorticity_values = omega.x.array
vorticity_coords = Q.tabulate_dof_coordinates()

fig, ax = plt.subplots(figsize=(12, 6))
plt.tricontourf(
    vorticity_coords[:, 0],
    vorticity_coords[:, 1],
    vorticity_values,
    cmap="RdBu_r",
    levels=500,
)
plt.savefig(
    "./results/vorticity.png",
    dpi=300,
    bbox_inches="tight",
)

The whole workflow is:

  1. Create mesh and save to xdmf.
  2. Solve the channel flow problem and write to .bp using VTK
  3. Load mesh and extract snapshots from .bp data
  4. Create function space V of CG2 for u and assign extracted values to u
  5. Create function space Q of DG1 for vorticity and interpolate the vorticity expression to compute vorticity.

The paraview computed vorticity

image
image1318×1022 27.5 KB

Thank you for your help!

4

And my computed vorticity:

vorticity
vorticity2811×1533 77.8 KB

5

I would strongly advice you to use adios4dolfinx to create your snapshots, as you are trying to reinvent the wheel, and is not taking into account that the mesh data is redistributed when read in.
See: ADIOS4DOLFINx - A framework for checkpointing in DOLFINx — ADIOS2Wrappers for various demos of different checkpoints.

6

Thank you! I have successfully computed my vorticity field by referring to the read_function example.

Another question I have will be, is it possible to read multiple time steps at once? Or to load a .bp file containing 200 snapshots, we need to access the file 200 times?

Thanks for the help!

7

I have not thought of this as a usecase, so it is currently not possible. It would be quite a memory consuming procedure, as you would have to store 200 x a single function.