Instability when solving the linear system

Linear Algebra
1

Hello everyone,

I’m trying to implement the mixed method from this paper for the Linear Elasticity problem. When it comes to solving the resulting linear system, I’m using

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

which, from my understanding, employs the standard LU solver.

By doing this, I get a NaN solution. However, if I change the solver package to superlu

petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "superlu", "ksp_monitor": None}

the method now works as expected.

Can someone help me figure out what is happening? I’m using the 0.7.0 version of dolfinx installed through conda. Here is the full code for context

import dolfinx.mesh as mesh
from dolfinx.fem import (FunctionSpace, assemble_scalar, form)
from dolfinx.fem.petsc import LinearProblem

import ufl
from ufl import (inner, outer, sym, div, grad, tr, dx, Identity, TrialFunction, 
                 TestFunction, FiniteElement, VectorElement, MixedElement,
                 split, as_vector)

from mpi4py import MPI
import numpy as np

# Problem definition ==========================================================

# Domain (rectangle)
cornerA = [0, 0]
cornerB = [1, 1]

# Mesh
msh = mesh.create_rectangle(MPI.COMM_WORLD, [np.array(cornerA), np.array(cornerB)],
                            [30, 30], cell_type=mesh.CellType.quadrilateral)

# Lamé parameters
mu = 0.3
lamb = 1

# Elasticity tensor
def tensorC(u):
    return 2*mu*sym(grad(u)) + lamb*div(u)*Identity(2)

# Compliance tensor
def tensorA(s):
    return (1/(2*mu))*(s - (lamb/(2*lamb+2*mu))*tr(s)*Identity(2))

# Exact solutions
def u_aux(mod):
    return [ufl.sin(ufl.pi*x[0])*ufl.sin(ufl.pi*x[1]),
            ufl.sin(ufl.pi*x[0])*ufl.sin(ufl.pi*x[1])]

x  = ufl.SpatialCoordinate(msh)
u_ufl = ufl.as_vector(u_aux(x))
sig_ufl = tensorC(u_ufl)
r_ufl = (grad(u_ufl)[0,1] - grad(u_ufl)[1,0])/2

# Load function
load = div(sig_ufl)

# Mixed method with wak symmetry ==============================================

# Approximation spaces

degree_s = 2
degree_u = 1
degree_r = 1

el_sig = FiniteElement('RTCF', msh.ufl_cell(), degree_s)
el_u = VectorElement('DQ', msh.ufl_cell(), degree_u)
el_rot = FiniteElement('DPC',msh.ufl_cell(), degree_r)
W = FunctionSpace(msh,MixedElement([el_sig, el_sig, el_u, el_rot]))

# Test and trial functions
e1 = as_vector([1, 0])
e2 = as_vector([0, 1])
TrialF = TrialFunction(W)
TestF  = TestFunction(W)
(sig1, sig2, u, r) = split(TrialF)
(tau1, tau2, v, q) = split(TestF)
sig = outer(e1, sig1) + outer(e2,sig2)
tau = outer(e1, tau1) + outer(e2, tau2)

# Variational problem
a = inner(tensorA(sig), tau)*dx
a += inner(u,div(tau))*dx 
a += r*(tau[0,1]-tau[1,0])*dx
a += inner(v,div(sig))*dx
a += q*(sig[0,1]-sig[1,0])*dx
L = inner(load,v)*dx

# Linear system solver

petsc_options={"ksp_type": "preonly", "pc_type": "lu", "ksp_monitor": None}
#petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "superlu", "ksp_monitor": None}

problem = LinearProblem(a, L, bcs=[], petsc_options=petsc_options)
wh = problem.solve()

(sig1h, sig2h, uh, rh) = wh.split()
sigh = outer(e1, sig1h) + outer(e2, sig2h)

# Approximation errors ========================================================

esig = sig_ufl - sigh
eu = u_ufl - uh
er = r_ufl - rh

eL2sig = np.sqrt(assemble_scalar(form(inner(esig,esig)*dx)))
eL2u = np.sqrt(assemble_scalar(form(inner(eu,eu)*dx)))
eL2r = np.sqrt(assemble_scalar(form(inner(er,er)*dx)))

print(eL2sig, eL2u, eL2r)

# =============================================================================
2

I don’t have access to the paper you’re citing so I can’t be sure of the weak formulation you’re attempting to discretise. But it looks like your underlying discretised linear operator has a non trivial nullspace. Some solver implementations (including direct solvers) are not well equipped to be used as a “black box” in this case. You’d have to carefully inform the solver of the nullspace. In the case of superlu see, e.g., MATSOLVERSUPERLU — PETSc v3.20.3-550-gc9d73d4ad9 documentation.

You could also try the linear solver MUMPS which I’ve used to great success as a black box for solving problems with a non trivial nullspace MATSOLVERMUMPS — PETSc 3.20.3 documentation.

2 Likes
3

Hey Nate, thanks for the quick reply!

The paper I mentioned is also available on the personal website of one of the authors (https://www-users.cse.umn.edu/~arnold//papers/quadelas.pdf).

From my understanding, the resulting discrete problem should not have a non-trivial nullspace (if the proper approximation spaces are being used), but maybe I’m missing something.

I double-checked the variational form, and it appears to be right. I’m not so sure about the implementation of the approximation spaces though.

In the meantime, I will look into the documentation of superlu and mumps as you advised.