Hello dolfinx users, I have a problem that requires me to compute eigenvalues, so dolfinx wrappers and petsc4py are not enough, I need to leverage slepc4py.
Physically, I want to do resolvant analysis at frequency 0 on a simple shear flow. This requires linearising Navier-Stokes, but also using an extensor matrix and manually setting boundary conditions outside of FENICsx wrappers.
Mathematically, I’m solving a generalised eigenvalue problem of the like of : B^HJ^{-1H}NJBf=\sigma Mf
B is an extensor so that Bf=\left[\begin{array}{c}f\\0\end{array}\right]
J is the Jacobian of Navier-Stokes with a trivial baseflow.
N is a norm matrix that satisfies : q^HNq=[u^H\;p^H]N\left[\begin{array}{c}u\\p\end{array}\right]=\int |u|^2 dx
Likewise, M is a norm matrix. f^HMf=\int |f|^2 dx
I’ve had mixed success - my code runs in serial and produces a nice graph :
However, the same code breaks apart in parallel :
Here’s my MWE. Apologies in advance for length, I did my best to boil it down.
import ufl
import numpy as np
import dolfinx as dfx
from dolfinx.io import XDMFFile
from petsc4py import PETSc as pet
from slepc4py import SLEPc as slp
from mpi4py.MPI import COMM_WORLD as comm
mesh=dfx.UnitSquareMesh(comm,100,100)
FE_vector_1=ufl.VectorElement("Lagrange",mesh.ufl_cell(),1)
FE_vector_2=ufl.VectorElement("Lagrange",mesh.ufl_cell(),2)
FE_scalar =ufl.FiniteElement("Lagrange",mesh.ufl_cell(),1)
# Taylor Hood elements
V =dfx.FunctionSpace(mesh,FE_vector_1)
TH=dfx.FunctionSpace(mesh,FE_vector_2*FE_scalar)
# Extract velocities only
U, _ = TH.sub(0).collapse(collapsed_dofs=True)
q,l=ufl.TrialFunction(TH),ufl.TestFunction(TH)
ub=dfx.Function(V)
w,z=ufl.TrialFunction(U), ufl.TestFunction(U)
u,p=ufl.split(q); v,s=ufl.split(l)
# Baseflow
if comm.size==1:
x=mesh.geometry.x[:,:2]
ub.vector[::2]=-1+2*x[:,1]
viewer = pet.Viewer().createMPIIO("sanity_check_baseflow.dat", 'w', comm)
ub.vector.view(viewer)
else:
viewer = pet.Viewer().createMPIIO("sanity_check_baseflow.dat", 'r', comm)
ub.vector.load(viewer)
# Mass csv
J_form = ufl.inner(ufl.div(u), s)
# Momentum (different test functions and IBP)
J_form += ufl.inner(ufl.grad(ub)*u, v) # Convection
J_form += ufl.inner(ufl.grad(u)*ub, v)
J_form += ufl.inner(ufl.grad(u),ufl.grad(v))/100 # Diffusion (Re=100)
J_form -= ufl.inner(p, ufl.div( v)) # Pressure
J_form *= ufl.dx
# Velocity norm
N_form = ufl.inner(u,v)*ufl.dx
# Forcing norm
M_form = ufl.inner(w,z)*ufl.dx
# Extensor
B_form = ufl.inner(w,v)*ufl.dx
# Assembling matrices
J = dfx.fem.assemble_matrix(J_form); J.assemble() # Jacobian
N = dfx.fem.assemble_matrix(N_form); N.assemble() # Mass matrix for q
M = dfx.fem.assemble_matrix(M_form); M.assemble() # Mass matrix for f
B = dfx.fem.assemble_matrix(B_form); B.assemble() # Extensor
# Boundary conditions
def border(x): return np.isclose(x[0],0)+np.isclose(x[0],1)+np.isclose(x[1],0)+np.isclose(x[1],1)
dofs,_=dfx.fem.locate_dofs_geometrical((TH.sub(0),U),border)
J.zeroRows(dofs); B.zeroRows(dofs,0)
# Sizes
n=B.getSize()[1]
n_local=B.getLocalSize()[1]
def setupKSP(KSP):
KSP.setType('preonly')
# Preconditioner
PC = KSP.getPC(); PC.setType('lu')
PC.setFactorSolverType('mumps')
KSP.setFromOptions()
# Compute solvers instead of inverses
KSPs = []
for Mat in [J,J.copy().hermitianTranspose()]:
# Krylov subspace
KSP = pet.KSP().create(comm)
KSP.setOperators(Mat)
setupKSP(KSP)
KSPs.append(KSP)
w,z = J.getVecs()
# Resolvent operator
class LHS_class:
def mult(cls,A,x,y):
B.mult(x,w)
KSPs[0].solve(w,z)
N.mult(z,w)
KSPs[1].solve(w,z)
B.multTranspose(z,y)
LHS=pet.Mat().create(comm)
LHS.setSizes([[n_local,n]]*2)
LHS.setType(pet.Mat.Type.PYTHON)
LHS.setPythonContext(LHS_class())
LHS.setUp()
# Eigensolver
EPS = slp.EPS(); EPS.create(comm)
EPS.setOperators(LHS,M) # Solve B^T*L^-1H*L^-1*B*f=sigma^2*M*f (cheaper than a proper SVD)
EPS.setProblemType(slp.EPS.ProblemType.GHEP) # Specify that A is hermitian (by construction), & M is semi-positive
EPS.setDimensions(1,10) # Find k eigenvalues only with max number of Lanczos vectors
EPS.setTolerances(1e-6,100) # Set absolute tolerance and number of iterations
# Spectral transform
ST = EPS.getST(); ST.setType('shift'); ST.setShift(0)
# Krylov subspace
KSP = ST.getKSP()
setupKSP(KSP)
EPS.setFromOptions()
EPS.solve()
f=dfx.Function(U)
EPS.getEigenvector(0,f.vector)
q=dfx.Function(TH)
B.mult(f.vector,w)
KSPs[0].solve(w,q.vector)
u,_=q.split()
with XDMFFile(comm, "sanity_check_response.xdmf","w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(u)
Running the above code in the dolfinx 3.1.0 docker container in serial then in parallel (complex mode) should reproduce above behaviour.
I strongly suspect I’ve forgotten to ghost a vector or messed up my information-sharing between processors.
I couldn’t reproduce the issue with a single KSP so I’m stuck with considering how petsc4py and slepc4py interact…
