Unable to set given (local) rows to identity matrix

1

Hello,In the following minimal running example, I am attempting to run the weak form over a rectangular domain, where the left and right edges have prescribed values. u is a vector field while P is a scalar field, I have cut off the piece associated with P in order to make the code minimal, so now it is essentially a 2D elastic stretch problem. I ran into an issue associated with applying BCs on u, the error is listed after the example code, could anyone help me with this issue? Thank you!

from dolfin import *
set_log_level(LogLevel.ERROR)

# Class representing the intial conditions
class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = 0.0 
        values[1] = 0.0
        values[2] = 0.0
    def value_shape(self):
        return (3,)

# Class for interfacing with the Newton solver
class Equation(NonlinearProblem):
    def __init__(self, L, a, bc_u_left, bc_u_right):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
        self.bc_u_left = bc_u_left
        self.bc_u_right = bc_u_right

    def F(self, b, x):
        assemble(self.L, tensor=b)
        self.bc_u_left.apply(b,x)
        self.bc_u_right.apply(b,x)

    def J(self, A, x):
        assemble(self.a, tensor=A)
        self.bc_u_left.apply(A)
        self.bc_u_right.apply(A)

# Boundary conditions
def left(x, on_boundary):
    return near(x[0],0.)
def right(x, on_boundary):
    return near(x[0],1.)

mesh = RectangleMesh(Point(0.,0.),Point(1.0, 0.1),100,10)

P1 = VectorElement("CG", mesh.ufl_cell(), 2)
P2 = FiniteElement("CG", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, MixedElement( [P1, P2] ))

dSol  = TrialFunction(ME)
testfunc = TestFunction(ME)
v,w = split(testfunc)
Sol   = Function(ME)
Sol0  = Function(ME)
du, dP = split(dSol)
u, P = Sol.split(True)
u0, P0 = split(Sol0)
Sol_init = InitialConditions(degree=2)
Sol.interpolate(Sol_init)
Sol0.interpolate(Sol_init)

def eps(u):return sym(nabla_grad(u))
sigma = eps(u) + div(u)*Identity(u.geometric_dimension())
L = inner( sigma, eps(v) )*dx
a = derivative(L, Sol, dSol)

g_u  = Constant((0.0,0.0))
bc_u_left  = DirichletBC(ME.sub(0), g_u, left)
Ubc   = 0.0
u_bc  = Constant((Ubc,0.0))
bc_u_right  = DirichletBC(ME.sub(0), u_bc, right)

problem = Equation(L, a, bc_u_left, bc_u_right)

Ubc += 5.0*10.0**-4.81
u_bc.assign(Constant((Ubc,0.0)))
Sol0.vector()[:] = Sol.vector()
solver = NewtonSolver()
solver.solve(problem, Sol.vector())

Traceback (most recent call last):
  File "x.py", line 88, in <module>
    solver.solve(problem, Sol.vector())
  File "x.py", line 34, in J
    self.bc_u_left.apply(A)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to set given (local) rows to identity matrix.
*** Reason:  some diagonal elements not preallocated (try assembler option keep_diagonal).
*** Where:   This error was encountered inside PETScMatrix.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  15b823f2c45d3036b6af931284d0f8e3c77b6845
*** -------------------------------------------------------------------------

[1]+  Done                    clear
2

You can get around this by doing the following:

    def J(self, A, x):
        #assemble(self.a, tensor=A)
        assemble(self.a, tensor=A, keep_diagonal=True)
        self.bc_u_left.apply(A)
        self.bc_u_right.apply(A)
1 Like
3

Thanks, David! Would you please explain why this keep_diagonal fixes this issue? Was my BCs definition insufficient? or the weak form? The keep_diagonal option did fix the issue on my full code, but I ended up with a convergence error, so I wanted to make sure it doesn’t come from the Bus

4

The reason this error happens is usually that there’s a zero block on the diagonal of your matrix. In this case, one’s initial guess would be that this happens because L doesn’t involve a P–w term. However, adding P*w*dx to L doesn’t fix the problem. The reason this doesn’t work is the use of Sol.split(True) instead of split(Sol) to define u and P from Sol. The first creates new Function objects that UFL won’t recognize as components of Sol (so derivative(L, Sol, dSol) will just be zero), while the second gives you symbols referring to components of Sol, which is what you want. (See the discussion here.)

1 Like