I added .x.scatter_forward everywhere, it’s different but not really better…
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 = dfx.Function(TH), dfx.Function(TH)
s = dfx.Function(U)
# Resolvent operator
class LHS_class:
def mult(cls,A,x,y):
B.mult(x,w.vector)
w.x.scatter_forward()
KSPs[0].solve(w.vector,z.vector)
z.x.scatter_forward()
N.mult(z.vector,w.vector)
w.x.scatter_forward()
KSPs[1].solve(w.vector,z.vector)
z.x.scatter_forward()
B.multTranspose(z.vector,s.vector)
s.x.scatter_forward()
s.vector.copy(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.vector)
w.x.scatter_forward()
KSPs[0].solve(w.vector,q.vector)
q.x.scatter_forward()
u,_=q.split()
with XDMFFile(comm, "sanity_check_response.xdmf","w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(u)
Gives the spurious but different :
