How to get the residual from KrylovSolver for linear elasticity

Linear Algebra
1

I’m solving a linear elasticity problem and I want to get the residual from the solver. I created a function that calculates it from the assembled matrices and the solution, but I think that it’s something that I should be able to get directly from the KrylovSolver object without doing a manual calculation. I didn’t find anything in the documentation so I thought I’d ask here.

from fenics import *
from fenics_adjoint import *
import numpy as np 
from ufl import nabla_div

# Define parameters
nelx = 30
nely = 15
nelz = 15
lx = float(nelx)
ly = float(nely)
lz = float(nelz)
E = 1.04
v = 0.3
mu = Constant(E / (2*(1+v)))
lmbda = Constant(E*v / ((1+v)*(1-2*v))) 

# Create a mesh
mesh = BoxMesh.create(
    [Point(0.0, 0.0, 0.0), Point(lx, ly, lz)], # define opposing corners
    [nelx, nely, nelz], # number of elements in each direction
    CellType.Type.hexahedron)
mesh = Mesh(mesh)

# FE function space and functions
V = VectorFunctionSpace(mesh, "CG", 1)
u_sol = Function(V) 
u = TrialFunction(V)
v = TestFunction(V)

# boundary condition
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.)
left = Left()
bc = DirichletBC(V, Constant((0., 0., 0.)), left)

# Stress and strain functions
def sigma(u):
    n = len(u) # size of u
    return lmbda*nabla_div(u)*Identity(n) + 2*mu*epsilon(u)

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

# Residual
def get_residual():
    A_np = A.array()
    u_np = u_sol.vector().get_local()
    b_np = b.get_local()
    res = np.dot(A_np, u_np) - b_np
    return res

# Define forms
f = Constant((0., 0., -1.))
L = dot(f, v) * dx
a = inner(sigma(u), nabla_grad(v)) * dx

# Solve
A, b = assemble_system(a, L, bc)
solver = KrylovSolver("cg", "ilu")
solver.solve(A, u_sol.vector(), b)
res = get_residual()

So by residual, I mean res = A*u - b where A is the stiffness matrix, b is the force vector and u is the solution. Thanks in advance.

2

Maybe this is what you are looking for ?
https://fenicsproject.org/qa/6482/there-return-array-residual-norms-from-petsckrylovsolver/

3

Thanks for the reply. :blush:

I am trying to do something very similar to that. It seems like the getConvergenceHistory() method that they reference returns the residual norm, but I want the want the residual vector. I tried looking for documentation on getConvergenceHistory() to see if it can be configured to return the vector instead of just the norm, but I didn’t find anything. Any idea where that documentation might be?

If not, that’s okay, this has been helpful.