Ideal workflow for saving parallel data?

dolfinx
1

In this example code, after solving parallelly I am plotting displacement field for a submesh.
What is the best way to store the parallel mesh and solution so that I can get the required plots for any other submesh later on??

import gmsh
import sys
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

import matplotlib.pyplot as plt
import dolfinx

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

gmsh.initialize()
gmsh.option.setNumber("General.Verbosity", 0)
gmsh.model.add("square")

gmsh.model.occ.addRectangle(0, 0, 0, 1, 1, 1)
gmsh.model.occ.addRectangle(1, 0, 0, 1, 1, 2)
gmsh.model.occ.addRectangle(0, 1, 0, 1, 1, 3)
gmsh.model.occ.addRectangle(1, 1, 0, 1, 1, 4)

gmsh.model.occ.fragment(
    [(2, 1)], [(2, 2), (2, 3), (2, 4)], removeObject=True, removeTool=True
)

gmsh.model.occ.synchronize()
lines = gmsh.model.getEntities(1)
surf = gmsh.model.getEntities(2)

for i in lines:
    gmsh.model.mesh.setTransfiniteCurve(i[1], 20)

for i in surf:
    gmsh.model.mesh.setTransfiniteSurface(i[1])
    gmsh.model.mesh.setRecombine(i[0], i[1])


gmsh.model.addPhysicalGroup(2, [1], 101)
gmsh.model.setPhysicalName(2, 101, "d1")

gmsh.model.addPhysicalGroup(2, [2], 102)
gmsh.model.setPhysicalName(2, 102, "d2")

gmsh.model.addPhysicalGroup(2, [3], 103)
gmsh.model.setPhysicalName(2, 103, "d3")

gmsh.model.addPhysicalGroup(2, [4], 104)
gmsh.model.setPhysicalName(2, 104, "d4")


gmsh.model.mesh.generate(2)

# if '-nopopup' not in sys.argv:
#     gmsh.fltk.run()
# gmsh.write("mesh.msh")
domain, markers, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)

gmsh.finalize()


####################################################################################################################


V = fem.functionspace(domain, ("Lagrange", 1))

uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0] ** 2 + 2 * x[1] ** 2)


# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)

Q = fem.functionspace(domain, ("DG", 0))
kappa = fem.Function(Q)
mA_cells = np.concatenate((markers.find(101), markers.find(102)))
kappa.x.array[mA_cells] = np.full_like(mA_cells, 1, dtype=default_scalar_type)
mB_cells = np.concatenate((markers.find(103), markers.find(104)))
kappa.x.array[mB_cells] = np.full_like(mB_cells, 0.1, dtype=default_scalar_type)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

f = fem.Constant(domain, default_scalar_type(-6))
a = ufl.inner(kappa * ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx

problem = LinearProblem(
    a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
uh = problem.solve()


####################################################################################################################


combined_indices = []
for tag in [101, 104]:
    combined_indices.append(markers.find(tag))
combined_indices = np.sort(np.hstack(combined_indices))


submesh, entity_map, vertex_map, geom_map = mesh.create_submesh(
    domain, domain.topology.dim, combined_indices
)

Vs = fem.functionspace(submesh, ("Lagrange", 1))
uDs = fem.Function(Vs)
uDs.interpolate(
    uh, cells0=entity_map, cells1=np.arange(len(entity_map), dtype=np.int32)
)

# Create local VTK mesh data structures
topology, cell_types, geometry = plot.vtk_mesh(Vs)
num_cells_local = domain.topology.index_map(domain.topology.dim).size_local
num_dofs_local = Vs.dofmap.index_map.size_local * Vs.dofmap.index_map_bs
if topology.shape[0] > 0:
    num_dofs_per_cell = topology[0]
else:
    num_dofs_per_cell = 0
# Get only dof indices
topology_dofs = (np.arange(len(topology)) % (num_dofs_per_cell + 1)) != 0
# Map to global dof indices
global_dofs = Vs.dofmap.index_map.local_to_global(topology[topology_dofs].copy())
# Overwrite topology
topology[topology_dofs] = global_dofs
# Gather data
global_topology = domain.comm.gather(
    topology[: (num_dofs_per_cell + 1) * num_cells_local], root=0
)
global_geometry = domain.comm.gather(
    geometry[: Vs.dofmap.index_map.size_local, :], root=0
)
global_ct = domain.comm.gather(cell_types[:num_cells_local])

global_def = domain.comm.gather(uDs.x.array[:num_dofs_local])

if rank == 0:
    # Stack data

    root_geom = np.vstack(global_geometry)
    root_top = np.concatenate(global_topology)
    root_ct = np.concatenate(global_ct)

    root_def = np.concatenate(global_def)

    u_grid = pyvista.UnstructuredGrid(root_top, root_ct, root_geom)
    u_grid.point_data["u"] = root_def
    pyvista.start_xvfb()
    u_plotter = pyvista.Plotter(off_screen=True)
    u_plotter.add_mesh(u_grid, show_edges=True)

    u_plotter.show()
    u_plotter.screenshot(f"u_{domain.comm.rank}_{domain.comm.size}.png")

I came across this where I can store the serialized data which can then be used to generate plots later on but this is not recommended.