I am running a CPP version of the poisson equation and want to modify some of the HYPRE parameters. In petsc4py I could do something like this:
python myprogram.py -pc_hypre_boomeramg_strong_threshold 0.75
but when I try to do something like this in the CPP version of FEniCS:
./myprogram --petsc.pc_hypre_boomeramg_strong_threshold 0.75
It doesn't work. Here's my code:
#include <dolfin.h>
#include "CG1.h"
#include "papi.h"
#include <math.h>
using namespace dolfin;
class Source : public Expression
{
void eval(Array<double>& values, const Array<double>& x) const
{
values[0] = 12*M_PI*M_PI*sin(2*M_PI*x[0])*sin(2*M_PI*x[1])*sin(2*M_PI*x[2]);
}
};
class DirichletBoundary : public SubDomain
{
bool inside(const Array<double>& x, bool on_boundary) const
{
return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS or x[2] < DOLFIN_EPS or x[2] > 1.0 - DOLFIN_EPS;
}
};
int main(int argc, char* argv[])
{
// Initialize stuff
int seed=80;
//seed = atoi(argv[1]);
// Create mesh and function space
auto mesh = std::make_shared<UnitCubeMesh>(seed,seed,seed);
auto V = std::make_shared<CG1::FunctionSpace>(mesh);
// Define boundary condition
auto u0 = std::make_shared<Constant>(0.0);
auto boundary = std::make_shared<DirichletBoundary>();
auto bc = std::make_shared<DirichletBC>(V, u0, boundary);
// Define variational forms
CG1::BilinearForm a(V, V);
CG1::LinearForm L(V);
auto f = std::make_shared<Source>();
L.f = f;
// Define linear algebra
auto A = std::make_shared<Matrix>();
Vector b;
KrylovSolver solver("cg", "hypre_amg");
//KrylovSolver solver(argv[2],argv[3]);
parameters.parse(argc,argv);
solver.parameters["relative_tolerance"] = 1.0e-7;
PETScOptions::set("ksp_converged_reason");
PETScOptions::set("ksp_view");
//PETScOptions::set("pc_hypre_boomeramg_strong_threshold", 0.75);
Function u(V);
// Compute solution
assemble_system(*A, b, a, L, {bc});
solver.set_operators(A,A);
solver.solve(*u.vector(),b);
// Save solution in VTK format
File file("poisson.pvd");
file << u;
// Plot solution
//plot(u);
//interactive();
return 0;
}
Any ideas what to do? I am using DOLFIN version 2016.2.0, here's my UFL file if needed:
element = FiniteElement("Lagrange", tetrahedron, 1)
u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)
a = inner(grad(u), grad(v))*dx
L = f*v*dx