Parameters in FEniCS 2018/2019 equivalent to LUSolver('mumps') from 2017

1

I have been using DOLFIN 2017.2.0 for a long while and it seems I might have to switch to one of the later releases, however, my code breaks when I run it under the later versions. It is kind of hard to tell how much of it breaks, since I will have to fix each part to go on, but here i my first problem:

I used to create a solver

solver = LUSolver('mumps')

which now breaks.
What is/are the equivalent statement(s) in DOLFIN 2018.1.0 and/or 2019.1.0?

2

If you look into the dolfin C++ code in dolfin-2018.1.0/dolfin/la/PETScLUSolver.h on Line 52, it reads as:

PETScLUSolver(std::string method="default")

And in the Pybind11 interface file (dolfin-2018.1.0/python/src/la.cpp) on Line 852, it reads as:

.def(py::init<std::string>(), py::arg("method")="default")

as part of the available python interface for PETScLUSolver for Python.

I’ve not dived into dolfin.LUSolver. But if you are using PETSc as the backend, maybe give PETScLUSolver directly a try. The only thing to note is it seems like there’s no interface for set_operator for PETScLUSolver, so you might need to do

solver = dolfin.PETScLUSolver(A, 'mumps')

as indicated by Lines 855~856 of Pybind11 interface code.

Hope this help.
Victor

2 Likes
3

If you want to have more specific control over petsc, take a look at my answer here:

There it is shown how you can use petsc directly, and how you can define which solver is chosen using “command line options” for petsc.

4

Hi, Victor,

Would you mind to share more about PETScLUSolver? I used to use PETScLUSolver.set_operator(A), A is la.GenericMatrix. However, I found PETScLUSolver(A) requires an input of PETScMatrix. Do you have any idea how la.GenericMatrix can be convert to PETScMatrix?

5

Now that some time has passed and I have used Victor’s advice I believe I have some qualification to help shed some more light on this too. At least I can show you how I have approached this. Everything here strictly operates on the native DOLFIN matrix type, typically by converting it to a PETScMatrix. Note that this is just copy-pasta from one of my project files.

from dolfin import __version__ as DOLFINversion

from dolfin import as_backend_type
from dolfin import Matrix, PETScMatrix
from dolfin import Vector

if DOLFINversion == '2017.2.0':
    from dolfin import mpi_comm_self
else:
    from dolfin import MPI
from scipy.sparse import csr_matrix
from petsc4py import PETSc

def transpose(M):
    """Transpose for matrix of DOLFIN type
    """
    M_ = as_backend_type(M).mat()
    Q = M_.copy()
    Q.transpose()
    return Matrix(PETScMatrix(Q))

def zeroVector(n):
    """Zero-vector of DOLFIN type
    """
    if DOLFINversion == '2017.2.0':
        return Vector(mpi_comm_self(),n)
    else:
        return Vector(MPI.comm_self,n)
    # use MPI.comm_self() in 2018.1.0 and later

def eye(n):
    """Identity Matrix of DOLFIN type
    """
    li = list(range(n+1))
    I = PETSc.Mat().createAIJWithArrays((n,n),(li,li[:-1],[1 for _ in li[1:]]))
    return Matrix(PETScMatrix(I))

def ToScipy(M):
    """Creates a sparse Scipy matrix from a DOLFIN type matrix
    """
    M_ = as_backend_type(M).mat()
    return csr_matrix(M_.getValuesCSR()[::-1],shape=M_.size)

def FromScipy(M):
    """Creates a DOLFIN type matrix from a sparse Scipy matrix
    """
    M = M.tocsr()
    IM = M.indptr
    JM = M.indices
    DM = M.data
    shape = M.shape
    M_ = PETSc.Mat().createAIJWithArrays(shape,(IM,JM,DM))
    return Matrix(PETScMatrix(M_))

def Solver(A):
    s = PETScLUSolver(as_backend_type(A),'mumps')
    return s
4 Likes