unable to use mumps as solver_parameter

Question

Hello,

I am working on a discontinuous Galerkin problem. My code seems like converge so far but when I increase n > 120 I get an error:

Unable to successfully call PETSc function 'KSPSolve'

Reason: PETSc error code is: 76.

I google to find an answer. And the most helping recommendation appears to use "mumps" as the linear solver.

But here is my issue: I was not able to change my code. I tried

solve(a == L, u, bc, solver_parameters={"linear_solver":"mumps"});

and

problem = VariationalProblem(a, L, u, bc)
solver = VariationalSolver(problem)
solver.parameters["linaer_solver"]="mumps"
solver.solve()

since my code is in cpp. But I think I am leaving out something. This time I am getting the error that I need to declare "problem" and "VariationalProblem". I tried to declare them in different ways but I cant figure out the correct way.

I am not familiar with cpp but this is an assignment I am solving for one of my courses at the university. I would be very happy if anyone can help me out to understand how to make the change so I can actually compile my code and get it running.

Thanks a lot :)

My cpp code reads as follows

    #include <dolfin.h>
    #include "Interior_DG.h"

    using namespace dolfin;

    #define DIM 2
    #define DIRICHLET_BOUNDARY_ID 1

    class Source : public Expression{
    private:
   double rho;
    public:
Source(double rho) : Expression(DIM) {
    this->rho = rho;
}

void eval(Array<double>& values, const Array<double>& x) const
{
    values[0] = 0.0;
    values[1] = -9.81*this->rho;
    if (DIM == 3) values[2] = 0.0;
}
    };

    class NeumannData : public Expression
    {
    private:
double gmax;
    public:
NeumannData(double gmax) : Expression(DIM) {
    this->gmax = gmax;
}

void eval(Array<double>& values, const Array<double>& x) const
{
    values[0] = 0.0;
    values[1] = x[1]*gmax;
    if (DIM == 3) values[2] = 0.0;
}
    };

    class DirichletBoundary : public SubDomain
    {
bool inside(const Array<double>& x, bool on_boundary) const
{
    return on_boundary && (x[0] <= 0.0 + DOLFIN_EPS);
}
    };

    class DirichletData : public Expression
    {
    public:
DirichletData() : Expression(DIM) { }

void eval(Array<double>& values, const Array<double>& x) const
{
    values[0] = 0.0;
    values[1] = 0.0;
    if (DIM == 3) values[2] = 0.0;
}
    };

    int main(int argc, char* argv[]){
Parameters params;
params.add("rho", 7850.0e-9);
params.add("gmax", 0.0);
    params.add("eta", -1.0);
    params.add("delta_e", 0.001);
    params.add("beta", 1.0/(DIM-1.0));
params.add("n", 100);
params.add("h", 0.1);
params.add("E", 210e3);
params.add("nu", 0.3);

params.parse(argc, argv);

info(params.str(true));

double h = params["h"];
std::size_t n = params["n"];

// Create mesh
    #if DIM == 2
RectangleMesh mesh(Point(0.0, 0.0), Point(1.0, h), n, (std::size_t) (n*h));
//RectangleMesh(p0, p1, number of cells in x direction, number of cells in y direction);
    #else
BoxMesh mesh(Point(0.0, 0.0, 0.0), Point(1.0, h, h), n, (std::size_t) (n*h), (std::size_t)(n*h));
    #endif

DirichletBoundary dirichlet_boundary;

// Create mesh_function
    MeshFunction<std::size_t> subdomain_ids(mesh, mesh.topology().dim()-1, 0);
subdomain_ids.set_all(0);
dirichlet_boundary.mark(subdomain_ids, DIRICHLET_BOUNDARY_ID);

// Define function space
Interior_DG::FunctionSpace V(mesh);

// Define variational forms
Interior_DG::BilinearForm a(V, V);
Interior_DG::LinearForm L(V);

// Assign source function and Neumann data
double rho = params["rho"];
double gmax = params["gmax"];
Source f(rho);
NeumannData g(gmax);
L.f = f;
L.g = g;

// Assign MeshFunction
a.ds = subdomain_ids;

// Assign DG constants
Constant eta(params["eta"]);
Constant delta_e(params["delta_e"]);
Constant beta(params["beta"]);
a.eta = eta;
a.delta_e = delta_e;
a.beta = beta;

// Assign Lame constants
Constant E(params["E"]);
Constant nu(params["nu"]);
Constant mu(E / (2*(1+nu)));
Constant lmbda(nu * E / ((1+nu)*(1-2*nu)));
a.lmbda = lmbda;
a.mu = mu;



// SubDomain for Dirichlet BC
DirichletData u0;
DirichletBC bc(V, u0, subdomain_ids, DIRICHLET_BOUNDARY_ID);

// Compute solution
Function u(V);


// Solve PDE
solve(a == L, u, bc);

    //  problem = VariationalProblem(a, L, u, bc)
    //  solver = VariationalSolver(problem)
    //  solver.parameters["linaer_solver"]="mumps"
    //  solver.solve()

// Plot solution
plot(u);
interactive();

return 0;
    }

Comments

Ugh, the indentation....
Other than that, it looks like a Hyperelasticity problem, have you looked at the demo? Whilst it doesn't use DG methods it still may be helpful with regard to the c++ API.
It may also help to post the UFL code that generates Interior_DG.h.

---

Thank you for your answer.

I looked at several demos. Unfortunately I could not find the answer.

Maybe it will help if I formulate my question differently. When I want to use

problem = LinearVariationalProblem(a, L, u, bc)
solver = LinearVariationalSolver(problem)
solver.parameters["linaer_solver"]="mumps"
solver.solve()

I am getting the error

'problem' was not declared in this scope
expected ';' before solver

How do I declare 'problem' ?

Thank you.

---

It seems you are trying to use the python API in c++. I'm not really familiar with the c++ API but perhaps something like

 LinearVariationalProblem problem(a,L,u,bc);
 LinearVariationalSolver solver(problem);

is what your looking for.