Direct solve fails for a mixed system

[ctjacobs]

I'm trying to use a direct solver on a mixed system:

from firedrake import *

mesh = UnitSquareMesh(2, 2)

# Define the mixed function space
U = VectorFunctionSpace(mesh, "CG", 2)
H = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([U, H])

# The solution field defined on the mixed function space
solution = Function(W)
u, h = split(solution)
w, v = TestFunctions(W)

# Define the compulsory shallow water fields
solution_old = Function(W)
solution_old.interpolate(Expression(("0.1*sin(pi*x[0])", "1e-16", "cos(pi*x[0])")))
u_old, h_old = split(solution_old)

# The solution should first hold the initial condition.
solution.assign(solution_old)

# Mean free surface height
h_mean = Function(W.sub(1)).interpolate(Expression("50"))

# The total height of the free surface.
h_total = h_mean + h

# Time-step size
dt = 0.05

# The form, F
F = 0
M_momentum = (1.0/dt)*(inner(w, u) - inner(w, u_old))*dx
F += M_momentum

A_momentum = inner(dot(grad(u), u), w)*dx
F += A_momentum

C_momentum = -9.81*inner(w, grad(h))*dx
F += C_momentum

M_continuity = (1.0/dt)*(inner(v, h) - inner(v, h_old))*dx
F += M_continuity

Ct_continuity = inner(v, div(h_total*u))*dx
F += Ct_continuity

# Construct the solver objects
problem = NonlinearVariationalProblem(F, solution, bcs=[])
solver = NonlinearVariationalSolver(problem, solver_parameters={'ksp_type':'preonly', 'pc_type':'lu'})
solver.solve()

but I'm getting the following PETSc error:

Traceback (most recent call last):
  File "test_lu.py", line 52, in <module>
    solver.solve()
  File "<string>", line 2, in solve
  File "/usr/local/lib/python2.7/dist-packages/PyOP2-0.11.0_148_g8a1955c-py2.7-linux-x86_64.egg/pyop2/profiling.py", line 197, in wrapper
    return f(*args, **kwargs)
  File "/data/firedrake/firedrake/variational_solver.py", line 252, in solve
    self.snes.solve(None, v)
  File "PETSc/SNES.pyx", line 511, in petsc4py.PETSc.SNES.solve (src/petsc4py.PETSc.c:155735)
petsc4py.PETSc.Error: error code 56
[0] SNESSolve() line 3740 in /tmp/pip-IWeiSF-build/src/snes/interface/snes.c
[0] SNESSolve_NEWTONLS() line 232 in /tmp/pip-IWeiSF-build/src/snes/impls/ls/ls.c
[0] KSPSolve() line 419 in /tmp/pip-IWeiSF-build/src/ksp/ksp/interface/itfunc.c
[0] KSPSetUp() line 306 in /tmp/pip-IWeiSF-build/src/ksp/ksp/interface/itfunc.c
[0] PCSetUp() line 902 in /tmp/pip-IWeiSF-build/src/ksp/pc/interface/precon.c
[0] PCSetUp_LU() line 116 in /tmp/pip-IWeiSF-build/src/ksp/pc/impls/factor/lu/lu.c
[0] MatGetOrdering() line 260 in /tmp/pip-IWeiSF-build/src/mat/order/sorder.c
[0] MatGetOrdering_ND() line 19 in /tmp/pip-IWeiSF-build/src/mat/order/spnd.c
[0] No support for this operation for this object type
[0] Cannot get rows for matrix type nest

Do I need to pass in some additional options when building PETSc, or is this unsupported in Firedrake? If it's the latter, please tag this as a 'wishlist' item.

[kynan]

Direct solvers are not supported on mixed systems by PETSc afair, so it's not a Firedrake limitation - correct me if that's wrong, @wence-.

[wence-]

It's not possible to use a direct solve with a MatNest, I have a PyOP2 branch that allows assembly into a monolithic matrix which would allow a direct solve (baij-matrices). Unfortunately it doesn't work with vector function spaces due to a petsc bug that I have reported without any response.

[ctjacobs]

Is the patch for this in petsc and petsc4py MAPDES-master? I just tried giving this a go with up-to-date packages, and locally merged master into the baij-matrices PyOP2 branch. I got a load of conflicts which I tried to fix manually. I'm still getting an error in petsc4py though:

Traceback (most recent call last):
  File "test.py", line 51, in <module>
    solver = NonlinearVariationalSolver(problem, solver_parameters={'ksp_type':'preonly', 'pc_type':'lu'})
  File "/data/firedrake/firedrake/variational_solver.py", line 115, in __init__
    ctx.set_jacobian(self.snes)
  File "/data/firedrake/firedrake/solving_utils.py", line 132, in set_jacobian
    snes.setJacobian(self.form_jacobian, J=self._jacs[-1]._M.handle,
  File "/usr/local/lib/python2.7/dist-packages/PyOP2-0.11.0_179_ga0030bb_dirty-py2.7-linux-x86_64.egg/pyop2/petsc_base.py", line 753, in handle
    self._init()
  File "/usr/local/lib/python2.7/dist-packages/PyOP2-0.11.0_179_ga0030bb_dirty-py2.7-linux-x86_64.egg/pyop2/petsc_base.py", line 351, in _init
    self._init_nest()
  File "/usr/local/lib/python2.7/dist-packages/PyOP2-0.11.0_179_ga0030bb_dirty-py2.7-linux-x86_64.egg/pyop2/petsc_base.py", line 548, in _init_nest
    mat.createNest([[m.handle for m in row_] for row_ in self._blocks],
  File "/usr/local/lib/python2.7/dist-packages/PyOP2-0.11.0_179_ga0030bb_dirty-py2.7-linux-x86_64.egg/pyop2/petsc_base.py", line 753, in handle
    self._init()
  File "/usr/local/lib/python2.7/dist-packages/PyOP2-0.11.0_179_ga0030bb_dirty-py2.7-linux-x86_64.egg/pyop2/petsc_base.py", line 355, in _init
    self._init_block()
  File "/usr/local/lib/python2.7/dist-packages/PyOP2-0.11.0_179_ga0030bb_dirty-py2.7-linux-x86_64.egg/pyop2/petsc_base.py", line 566, in _init_block
    (self.sparsity._rowptr, self.sparsity._colidx, self._array))
  File "PETSc/Mat.pyx", line 343, in petsc4py.PETSc.Mat.createAIJWithArrays (src/petsc4py.PETSc.c:109338)
ValueError: A matrix with 50 rows requires a row pointer of length 51 (given: 26)

[wence-]

There is a dangling change in petsc (not merged) and two open changes in pyop2 and firedrake, either of those give details.

[wence-]

Should now work, pass nest=False to the NLVProblem constructor (or nest=False to a solve(...) call) or, if you want to set the default parameters["matnest"] = False.

[ctjacobs]

Works for me. Thanks.

[agzimmerman]

I'm trying to use a direct solver with a mixed problem and the solution to this issue no longer seems correct.

firedrake.LinearVariationalProblem.__init__ does not have a nest keyword argument. Neither does firedrake.LinearVariationalSolver.solve. The same is true for the nonlinear counterparts. Setting firedrake.parameters["matnest"] = False does not change the behavior.

Here's a MWE based on the mixed Poisson tutorial:

import firedrake as fe


print("Using firedrake-" + fe.__version__)

mesh = fe.UnitSquareMesh(3, 3)

BDM = fe.FunctionSpace(mesh, "BDM", 1)

DG = fe.FunctionSpace(mesh, "DG", 0)

W = BDM * DG


sigma, u = fe.TrialFunctions(W)

tau, v = fe.TestFunctions(W)


x, y = fe.SpatialCoordinate(mesh)

f = fe.Function(DG).interpolate(
    10*fe.exp(-(pow(x - 0.5, 2) + pow(y - 0.5, 2)) / 0.02))

dot, div, dx = fe.dot, fe.div, fe.dx
    
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx

L = - f*v*dx


bc0 = fe.DirichletBC(W.sub(0), fe.as_vector([0.0, -fe.sin(5*x)]), 1)

bc1 = fe.DirichletBC(W.sub(0), fe.as_vector([0.0, fe.sin(5*y)]), 2)


w = fe.Function(W)

fe.parameters["matnest"] = False

problem = fe.LinearVariationalProblem(a, L, w, bcs=[bc0, bc1])

solver = fe.LinearVariationalSolver(
    problem, solver_parameters={'ksp_type':'preonly', 'pc_type':'lu'})

solver.solve()

This prints

Using firedrake-0.13.0+2351.gb7b1294
Traceback (most recent call last):
File "mixed_poisson.py", line 46, in
solver.solve()
File "/home/zimmerman/firedrake/src/firedrake/firedrake/variational_solver.py", line 239, in solve
work.copy(u)
File "/usr/lib/python3.5/contextlib.py", line 77, in exit
self.gen.throw(type, value, traceback)
File "/home/zimmerman/firedrake/src/firedrake/firedrake/petsc.py", line 199, in inserted_options
yield
File "/home/zimmerman/firedrake/src/firedrake/firedrake/variational_solver.py", line 239, in solve
work.copy(u)
File "/usr/lib/python3.5/contextlib.py", line 77, in exit
self.gen.throw(type, value, traceback)
File "/home/zimmerman/firedrake/src/PyOP2/pyop2/petsc_base.py", line 410, in vecscatter
yield self._vec
File "/home/zimmerman/firedrake/src/firedrake/firedrake/variational_solver.py", line 238, in solve
self.snes.solve(None, work)
File "PETSc/SNES.pyx", line 538, in petsc4py.PETSc.SNES.solve
petsc4py.PETSc.Error: error code 56
[0] SNESSolve() line 4396 in /tmp/pip-req-build-4g28hq0f/src/snes/interface/snes.c
[0] SNESSolve_KSPONLY() line 38 in /tmp/pip-req-build-4g28hq0f/src/snes/impls/ksponly/ksponly.c
[0] KSPSolve() line 723 in /tmp/pip-req-build-4g28hq0f/src/ksp/ksp/interface/itfunc.c
[0] KSPSetUp() line 391 in /tmp/pip-req-build-4g28hq0f/src/ksp/ksp/interface/itfunc.c
[0] PCSetUp() line 932 in /tmp/pip-req-build-4g28hq0f/src/ksp/pc/interface/precon.c
[0] PCSetUp_LU() line 84 in /tmp/pip-req-build-4g28hq0f/src/ksp/pc/impls/factor/lu/lu.c
[0] MatGetOrdering() line 250 in /tmp/pip-req-build-4g28hq0f/src/mat/order/sorder.c
[0] MatGetOrdering_ND() line 17 in /tmp/pip-req-build-4g28hq0f/src/mat/order/spnd.c
[0] No support for this operation for this object type
[0] Cannot get rows for matrix type nest

So with the latest Firedrake, how does one disable matrix nesting, or generally how does one use a direct solver with mixed FE problems?

[colinjcotter]

Dear Alex, For mixed problems with LU you need to use the PETSc option "mat_type":"aij" to tell PETSc to make one big matrix instead of storing the blocks separately. all the best

…

--cjc

________________________________ From: Alexander G. Zimmerman <notifications@github.com> Sent: 13 November 2018 08:54:12 To: firedrakeproject/firedrake Cc: Subscribed Subject: Re: [firedrakeproject/firedrake] Direct solve fails for a mixed system (#431) I'm trying to use a direct solver with a mixed problem and the solution to this issue no longer seems correct. If I try to pass firedrake.LinearVariationalProblem.__init__ does not have a nest keyword argument. Neither does firedrake.LinearVariationalSolver.solve. The same is true for the nonlinear counterparts. Setting fe.parameters["matnest"] = False does not change the behavior. Here's a MWE based on the mixed Poisson tutorial<https://www.firedrakeproject.org/demos/poisson_mixed.py.html>: import firedrake as fe print("Using firedrake-" + fe.__version__) mesh = fe.UnitSquareMesh(3, 3) BDM = fe.FunctionSpace(mesh, "BDM", 1) DG = fe.FunctionSpace(mesh, "DG", 0) W = BDM * DG sigma, u = fe.TrialFunctions(W) tau, v = fe.TestFunctions(W) x, y = fe.SpatialCoordinate(mesh) f = fe.Function(DG).interpolate( 10*fe.exp(-(pow(x - 0.5, 2) + pow(y - 0.5, 2)) / 0.02)) dot, div, dx = fe.dot, fe.div, fe.dx a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx L = - f*v*dx bc0 = fe.DirichletBC(W.sub(0), fe.as_vector([0.0, -fe.sin(5*x)]), 1) bc1 = fe.DirichletBC(W.sub(0), fe.as_vector([0.0, fe.sin(5*y)]), 2) w = fe.Function(W) fe.parameters["matnest"] = False problem = fe.LinearVariationalProblem(a, L, w, bcs=[bc0, bc1]) solver = fe.LinearVariationalSolver( problem, solver_parameters={'ksp_type':'preonly', 'pc_type':'lu'}) solver.solve() This prints Using firedrake-0.13.0+2351.gb7b1294 Traceback (most recent call last): File "mixed_poisson.py", line 46, in solver.solve() File "/home/zimmerman/firedrake/src/firedrake/firedrake/variational_solver.py", line 239, in solve work.copy(u) File "/usr/lib/python3.5/contextlib.py", line 77, in exit self.gen.throw(type, value, traceback) File "/home/zimmerman/firedrake/src/firedrake/firedrake/petsc.py", line 199, in inserted_options yield File "/home/zimmerman/firedrake/src/firedrake/firedrake/variational_solver.py", line 239, in solve work.copy(u) File "/usr/lib/python3.5/contextlib.py", line 77, in exit self.gen.throw(type, value, traceback) File "/home/zimmerman/firedrake/src/PyOP2/pyop2/petsc_base.py", line 410, in vecscatter yield self._vec File "/home/zimmerman/firedrake/src/firedrake/firedrake/variational_solver.py", line 238, in solve self.snes.solve(None, work) File "PETSc/SNES.pyx", line 538, in petsc4py.PETSc.SNES.solve petsc4py.PETSc.Error: error code 56 [0] SNESSolve() line 4396 in /tmp/pip-req-build-4g28hq0f/src/snes/interface/snes.c [0] SNESSolve_KSPONLY() line 38 in /tmp/pip-req-build-4g28hq0f/src/snes/impls/ksponly/ksponly.c [0] KSPSolve() line 723 in /tmp/pip-req-build-4g28hq0f/src/ksp/ksp/interface/itfunc.c [0] KSPSetUp() line 391 in /tmp/pip-req-build-4g28hq0f/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 932 in /tmp/pip-req-build-4g28hq0f/src/ksp/pc/interface/precon.c [0] PCSetUp_LU() line 84 in /tmp/pip-req-build-4g28hq0f/src/ksp/pc/impls/factor/lu/lu.c [0] MatGetOrdering() line 250 in /tmp/pip-req-build-4g28hq0f/src/mat/order/sorder.c [0] MatGetOrdering_ND() line 17 in /tmp/pip-req-build-4g28hq0f/src/mat/order/spnd.c [0] No support for this operation for this object type [0] Cannot get rows for matrix type nest So with the latest Firedrake, how does one disable matrix nesting, or generally how does one use a direct solver with mixed FE problems? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#431 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AF0q8rChAHgt05qNo1BJRyct4xMW3kxZks5uuoi0gaJpZM4DTzsg>.

[agzimmerman]

@colinjcotter , thanks! There was only one more thing missing, which was setting "pc_factor_mat_solver_type": "mumps". So it works if I instantiate the solver with

solver = fe.LinearVariationalSolver(
    problem, 
    solver_parameters = {
        "ksp_type": "preonly", 
        "pc_type": "lu", 
        "mat_type": "aij",
        "pc_factor_mat_solver_type": "mumps"})

Being able to search for "aij", I found this in the Navier-Stokes tutorial.