Hello,
I am running into an issue regarding function IO. I have a mesh saved using adios4dolfinx with the BP4 writer. I then have a program which solves a stochastic partial differential equation (SPDE) for a number of samples, and saves a fem.Function for each sample, also using adios4dolfinx and the BP4 writer. I would like to be able to add new sample results to an existing *.bp output in case I need more. Let me elaborate with a minimal working example.
First, we generate a mesh, for the MWE this is a 2D square mesh with the same number of cells in each direction. The script for generating it is shown below, let us call the file mesh.py. To generate the mesh run python3 mesh.py 30 ./mesh.bp.
mesh.py
from mpi4py import MPI
import sys
from pathlib import Path
import adios4dolfinx as adx
import dolfinx.mesh as mesh
if __name__ == "__main__":
assert (
len(sys.argv) == 3
), f"Usage: python3 {sys.argv[0]} <Number of cells for each direction> <Path to save mesh as *.bp>"
n_cells: int = int(sys.argv[1])
save_path: Path = Path(sys.argv[2])
domain: mesh.Mesh = mesh.create_rectangle(
comm=MPI.COMM_WORLD, points=[[-1.0, -1.0], [1.0, 1.0]], n=[n_cells, n_cells]
)
adx.write_mesh(filename=save_path, mesh=domain)
Next, we generate artificial fem.Function output. Here I picked the function f: \mathbb{R}^2 \to \mathbb{R}: (x, y) \mapsto x + y. We will do this twice, emulating two runs which generate one sample each by running mpiexec -np 4 python3 function.py ./mesh.bp FirstRun ./out.bp and mpiexec -np 4 python3 function.py ./mesh.bp SecondRun ./out.bp. The code for the function.py script is shown below. Note that on the first run, we read the mesh from mesh.bp and save it, together with the partitioning, in out.bp. On the second run we read the mesh from out.bp using the saved partitioning. The functions from both runs are now saved in out.bp under the names “FirstRun” and “SecondRun”.
function.py
from mpi4py import MPI
import sys
from pathlib import Path
import adios4dolfinx as adx
import dolfinx.mesh as mesh
import dolfinx.fem as fem
import numpy as np
import numpy.typing as npt
if __name__ == "__main__":
assert (
len(sys.argv) == 4
), f"Usage: python3 {sys.argv[0]} <Path to mesh *.bp directory> <Name of function> <Path to output *.bp directory>"
mesh_path: Path = Path(sys.argv[1])
fnc_name: str = str(sys.argv[2])
out_path: Path = Path(sys.argv[3])
if out_path.exists():
domain: mesh.Mesh = adx.read_mesh(
filename=out_path, comm=MPI.COMM_WORLD, read_from_partition=True
)
else:
domain: mesh.Mesh = adx.read_mesh(filename=mesh_path, comm=MPI.COMM_WORLD)
adx.write_mesh(filename=out_path, mesh=domain, store_partition_info=True)
element: fem.ElementMetaData = fem.ElementMetaData(family="Lagrange", degree=1)
V: fem.FunctionSpace = fem.functionspace(mesh=domain, element=element)
fnc = fem.Function(V)
assert isinstance(fnc, fem.Function)
def fnc_eval(x: npt.NDArray) -> npt.NDArray:
return np.sum(x, axis=0)
fnc.interpolate(u0=fnc_eval)
fnc.name = fnc_name
adx.write_function(filename=out_path, u=fnc)
Lastly, we visualize these two functions using pyvista by executing python3 view.py ./out.bp FirstRun SecondRun. The view.py script is shown below. What we then see is the function saved from the first run on the left and of the second run on the right. The first looks as you would expect, and the second looks wrong, though it still has a resemblance to what it should look like.
view.py
from mpi4py import MPI
import sys
from pathlib import Path
import adios4dolfinx as adx
import dolfinx.mesh as mesh
import dolfinx.fem as fem
from dolfinx.plot import vtk_mesh
import numpy as np
import pyvista as pv
if __name__ == "__main__":
assert (
len(sys.argv) == 4
), f"Usage: python3 {sys.argv[0]} <Path to mesh & data *.bp directory> <Name of first function> <Name of second function>"
mesh_path: Path = Path(sys.argv[1])
fnc_name_1: str = str(sys.argv[2])
fnc_name_2: str = str(sys.argv[3])
domain: mesh.Mesh = adx.read_mesh(filename=mesh_path, comm=MPI.COMM_WORLD)
element: fem.ElementMetaData = fem.ElementMetaData(family="Lagrange", degree=1)
V: fem.FunctionSpace = fem.functionspace(mesh=domain, element=element)
fnc_1, fnc_2 = fem.Function(V), fem.Function(V)
assert isinstance(fnc_1, fem.Function) and isinstance(fnc_2, fem.Function)
adx.read_function(filename=mesh_path, u=fnc_1, name=fnc_name_1)
adx.read_function(filename=mesh_path, u=fnc_2, name=fnc_name_2)
pv.set_plot_theme("paraview")
pl = pv.Plotter(shape=(1, 2))
grid = pv.UnstructuredGrid(*vtk_mesh(V))
grid_1 = grid.copy()
grid_1.point_data[fnc_name_1] = np.real(fnc_1.x.array)
grid_2 = grid.copy()
grid_2.point_data[fnc_name_2] = np.real(fnc_2.x.array)
pl.subplot(0, 0)
pl.show_axes()
pl.camera_position = "xy"
pl.zoom_camera(0.25)
pl.add_mesh(grid_1, show_edges=True, show_scalar_bar=True)
pl.subplot(0, 1)
pl.show_axes()
pl.camera_position = "xy"
pl.zoom_camera(0.25)
pl.add_mesh(grid_2, show_edges=True, show_scalar_bar=True)
pl.show()
So the question is, what am I doing wrong here? Thanks in advance for the help.
I am using dolfinx version 0.10.0.post1 and adios4dolfinx version 0.10.0.post0.
For convenience, here is a run.sh bash script for running the MWE assuming the *.py files mentioned are present.
run.sh
#!/bin/bash
python3 mesh.py 30 ./mesh.bp
mpiexec -np 4 python3 function.py ./mesh.bp FirstRun ./out.bp
mpiexec -np 4 python3 function.py ./mesh.bp SecondRun ./out.bp
python3 view.py ./out.bp FirstRun SecondRun
