# Compute vorticity from velocity snapshots obtained with ADIOS2 IO

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):
steps = []
for step in fh:
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.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:

Any help is appreciated!

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.

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,
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):

# 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)
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])
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])
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)
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.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,
rhs,
sym,
system,
)

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

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

steps = []
for step in fh:
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"

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

with io.XDMFFile(
MPI.COMM_WORLD,
"path/to/channel_flow_mesh.xdmf",
"r",
) as xdmf:

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

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