Hello All,
I am using using Fenicsx code for the past few days, and I appreciate that it is available for free use. I am attempting a buckling analysis, but I am encountering issues with visualization in Paraview. When I execute the code without MPI, the visualization appears fine; however, when I utilize multiple processors, the Paraview image becomes fragmented. This seems to resemble the pattern of mesh distribution across the MPI processors. This issue occurs with both VTK and VTX formats, but not with XDMF. However, the XDMF format has its limitations as it only supports first or second-order Lagrange elements, necessitating that the degree is equal to the order of the mesh geometry. It has been discussed on Interpolating Functions on higher order elements. So I want to use VTX or VTK.
To reproduce the problem with a simple example, I followed the buckling analysis on:
Linear Buckling Analysis of a 3D solid — Computational Mechanics Numerical Tours with FEniCSx as shown below and tried to export the eigen modes to view in Paraview.
Here is the maincode:
import numpy as np
import ufl
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc
from dolfinx import io, mesh, fem, default_scalar_type
from dolfinx.fem.petsc import LinearProblem, assemble_matrix
from eigenvalue_solver import solve_GEP_shiftinvert, EPS_get_spectrum
# Setup the work directory where the resuts are stored
work_dir = "./"
# Define geometry (All dimensions in meters)
b, w, length = 0.01, 0.03, 1.0
Nx, Ny, Nz = 5, 5, 50
domain = mesh.create_box(MPI.COMM_WORLD, [[0, 0, 0], [b, w, length]], [Nx, Ny, Nz], mesh.CellType.hexahedron)
gdim = domain.geometry.dim
# Define material properties, Assumed Isotropic
E = 210e9 # Young's modulus (Pa)
nu = 0.3 # Poisson's ratio
mu = fem.Constant(domain, default_scalar_type(E / 2 / (1 + nu)))
lmbda = fem.Constant(domain, default_scalar_type(E * nu / (1 + nu) / (1 - 2 * nu)))
def eps(v):
return ufl.sym(ufl.grad(v))
def sigma(v):
return lmbda * ufl.tr(eps(v)) * ufl.Identity(gdim) + 2 * mu * eps(v)
# Solving for prestressed state
N0 = 1
T = fem.Constant(domain, default_scalar_type((0, 0, -N0)))
#V = fem.functionspace(domain, ("P", 2, (gdim,)))
#V = fem.functionspace(domain, ("Lagrange", 2, (gdim,)))
V = fem.functionspace(domain, ("Lagrange", 1, (gdim,)))
v = ufl.TestFunction(V)
du = ufl.TrialFunction(V)
# Define boundary conditions
def bottom(x):
return np.isclose(x[2], 0.0)
def top(x):
return np.isclose(x[2], length)
fdim = domain.topology.dim - 1
top_facets = mesh.locate_entities_boundary(domain, fdim, top)
bottom_facets = mesh.locate_entities_boundary(domain, fdim, bottom)
u_L = np.array([0, 0, 0], dtype=default_scalar_type)
bcs = [
fem.dirichletbc(u_L, fem.locate_dofs_topological(V, fdim, bottom_facets), V),
fem.dirichletbc(
default_scalar_type(0),
fem.locate_dofs_topological(V.sub(0), fdim, top_facets),
V.sub(0),
),
fem.dirichletbc(
default_scalar_type(0),
fem.locate_dofs_topological(V.sub(1), fdim, top_facets),
V.sub(1),
),
]
# Mark right facets for Neumann boundary condition application
facet_tag = mesh.meshtags(domain, fdim, top_facets, np.full_like(top_facets, 1))
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
# Linear Elasticity Bilinear Form
a = ufl.inner(sigma(du), eps(v)) * ufl.dx
# External loading on the right end
l = ufl.dot(T, v) * ds(1)
# Solution of the trivial problem
u = fem.Function(V, name="Displacement")
problem = LinearProblem(
a, l, u=u, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
problem.solve()
# Stiffness matrix
K = assemble_matrix(fem.form(a), bcs=bcs, diagonal=1)
K.assemble()
write_option = 'VTX'
if write_option == 'VTX':
with io.VTXWriter(MPI.COMM_WORLD, f"{work_dir}/solution.bp", u) as vtx:
vtx.write(0)
elif write_option == 'XDMF':
with io.XDMFFile(MPI.COMM_WORLD, work_dir + "/solution.xdmf", 'w') as xdmf:
xdmf.write_mesh(domain)
xdmf.write_function(u, 0)
# EIGEN BUCKLING ANALYSIS
kgform = -ufl.inner(sigma(u), ufl.grad(du).T * ufl.grad(v)) * ufl.dx
KG = assemble_matrix(fem.form(kgform), bcs=bcs, diagonal=0)
KG.assemble()
# Requested number of eigenvalues
N_eig = 6
# Solve eigenvalue problem
eigensolver = solve_GEP_shiftinvert(
K,
KG,
problem_type=SLEPc.EPS.ProblemType.GHEP,
solver=SLEPc.EPS.Type.KRYLOVSCHUR,
nev=N_eig,
tol=1e-12,
shift=1e-2,
)
# Extract eigenpairs
(eigval, eigvec_r, eigvec_i) = EPS_get_spectrum(eigensolver, V)
# Write eigenmodes
for i in range(N_eig):
eigvec_r[i].name = "Eigenmode {}".format(i+1)
if write_option == 'VTX':
with io.VTXWriter(MPI.COMM_WORLD, f"{work_dir}/eigenmodes.bp", eigvec_r) as vtx:
vtx.write(0)
elif write_option == 'XDMF':
with io.XDMFFile(MPI.COMM_WORLD, f"{work_dir}/eigenmodes.xdmf", 'w') as xdmf:
xdmf.write_mesh(domain)
for i in range(N_eig):
xdmf.write_function(eigvec_r[i], i+1)
and here is the code to solve the eigenvalue problem (as it is) from Linear Buckling Analysis of a 3D solid — Computational Mechanics Numerical Tours with FEniCSx
from typing import List, Tuple
import dolfinx.fem as fem
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc
from dolfinx.fem import FunctionSpace
def print0(string: str):
"""Print on rank 0 only"""
if MPI.COMM_WORLD.rank == 0:
print(string)
def monitor_EPS_short(
EPS: SLEPc.EPS, it: int, nconv: int, eig: list, err: list, it_skip: int
):
"""
Concise monitor for EPS.solve().
Parameters
----------
eps
Eigenvalue Problem Solver class.
it
Current iteration number.
nconv
Number of converged eigenvalue.
eig
Eigenvalues
err
Computed errors.
it_skip
Iteration skip.
"""
if it == 1:
print0("******************************")
print0("*** SLEPc Iterations... ***")
print0("******************************")
print0("Iter. | Conv. | Max. error")
print0(f"{it:5d} | {nconv:5d} | {max(err):1.1e}")
elif not it % it_skip:
print0(f"{it:5d} | {nconv:5d} | {max(err):1.1e}")
def EPS_print_results(EPS: SLEPc.EPS):
"""Print summary of solution results."""
print0("\n******************************")
print0("*** SLEPc Solution Results ***")
print0("******************************")
its = EPS.getIterationNumber()
print0(f"Iteration number: {its}")
nconv = EPS.getConverged()
print0(f"Converged eigenpairs: {nconv}")
if nconv > 0:
# Create the results vectors
vr, vi = EPS.getOperators()[0].createVecs()
print0("\nConverged eigval. Error ")
print0("----------------- -------")
for i in range(nconv):
k = EPS.getEigenpair(i, vr, vi)
error = EPS.computeError(i)
if k.imag != 0.0:
print0(f" {k.real:2.2e} + {k.imag:2.2e}j {error:1.1e}")
else:
pad = " " * 11
print0(f" {k.real:2.2e} {pad} {error:1.1e}")
def EPS_get_spectrum(
EPS: SLEPc.EPS, V: FunctionSpace
) -> Tuple[List[complex], List[PETSc.Vec], List[PETSc.Vec]]:
"""Retrieve eigenvalues and eigenfunctions from SLEPc EPS object.
Parameters
----------
EPS
The SLEPc solver
V
The function space
Returns
-------
Tuple consisting of: List of complex converted eigenvalues,
lists of converted eigenvectors (real part) and (imaginary part)
"""
# Get results in lists
eigval = list()
eigvec_r = list()
eigvec_i = list()
for i in range(EPS.getConverged()):
vr = fem.Function(V)
vi = fem.Function(V)
eigval.append(EPS.getEigenpair(i, vr.vector, vi.vector))
eigvec_r.append(vr)
eigvec_i.append(vi) # Sort by increasing magnitude
idx = np.argsort(np.abs(np.array(eigval)), axis=0)
eigval = [eigval[i] for i in idx]
eigvec_r = [eigvec_r[i] for i in idx]
eigvec_i = [eigvec_i[i] for i in idx]
return (eigval, eigvec_r, eigvec_i)
def solve_GEP_shiftinvert(
A: PETSc.Mat,
B: PETSc.Mat,
problem_type: SLEPc.EPS.ProblemType = SLEPc.EPS.ProblemType.GNHEP,
solver: SLEPc.EPS.Type = SLEPc.EPS.Type.KRYLOVSCHUR,
nev: int = 10,
tol: float = 1e-7,
max_it: int = 10,
target: float = 0.0,
shift: float = 0.0,
) -> SLEPc.EPS:
"""
Solve generalized eigenvalue problem A*x=lambda*B*x using shift-and-invert
as spectral transform method.
Parameters
----------
A
The matrix A
B
The matrix B
problem_type
The problem type, for options see: https://bit.ly/3gM5pth
solver:
Solver type, for options see: https://bit.ly/35LDcMG
nev
Number of requested eigenvalues.
tol
Tolerance for slepc solver
max_it
Maximum number of iterations.
target
Target eigenvalue. Also used for sorting.
shift
Shift 'sigma' used in shift-and-invert.
Returns
-------
EPS
The SLEPc solver
"""
# Build an Eigenvalue Problem Solver object
EPS = SLEPc.EPS()
EPS.create(comm=MPI.COMM_WORLD)
EPS.setOperators(A, B)
EPS.setProblemType(problem_type)
# set the number of eigenvalues requested
EPS.setDimensions(nev=nev)
# Set solver
EPS.setType(solver)
# set eigenvalues of interest
EPS.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
EPS.setTarget(target) # sorting
# set tolerance and max iterations
EPS.setTolerances(tol=tol, max_it=max_it)
# Set up shift-and-invert
# Only work if 'whichEigenpairs' is 'TARGET_XX'
ST = EPS.getST()
ST.setType(SLEPc.ST.Type.SINVERT)
ST.setShift(shift)
EPS.setST(ST)
# set monitor
it_skip = 1
EPS.setMonitor(
lambda eps, it, nconv, eig, err: monitor_EPS_short(
eps, it, nconv, eig, err, it_skip
)
)
# parse command line options
EPS.setFromOptions()
# Display all options (including those of ST object)
# EPS.view()
EPS.solve()
EPS_print_results(EPS)
return EPS
Here is how the figure looks when running on a single processor for any format VTK, VTX or XDMF and also for multiple processors if using XDMF format. This is what is expected for the first eigenmode (scaling factor = 0.005)
But when I use multiple processors (mpirun -np 4) with VTK or VTX formats here is how it looks.
I understand that this is not a Fenicsx issue but a problem with the output formats. If I am doing something wrong or if others have faced similar issue, please help.