Logo
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ on Travis-CI Feel++ on Twitter Feel++ on YouTube Feel++ community
 All Classes Files Functions Variables Typedefs Pages
Using a backend

Introduction

After the discretization process, one may have to solve a (non) linear system. Feel++ interfaces with PETSc/SLEPc and Eigen3. Consider this system

\[A x = b \]

We call Backend an object that manages the solution strategy to solve it. Some explanation are available Solver and Preconditioner.

Feel++ provides a default backend that is mostly hidden to the final user. In many examples, you do not have to take care of the backend. You change the backend behavior via the command line or config files. For example

./feelpp_doc_mybackend --backend.pc-type=id

will use the identity matrix as a right preconditionner for the default backend. The size of the preconditionner will be defined from the size of the A matrix.

If you try to solve a different system \(A_1 y= c\) (in size) with the same backend or the default without rebuilding it, it will fail.

backend(_rebuild=true)->solve(_matrix=A1,_rhs=c,_sol=y);
Remarks
Each of that options can be retrieved via the –help-lib argument in the command line.

Non default backend

You may need to manage more than one backend in an application: you have different systems to solve and you want to keep some already computed objects such as preconditioners.

The default backend is in fact an unnamed backend: in order to distinguish between backend you have to name them. for example

po::options_description app_options( "MyBackend options" );
app_options.add(backend_options("myBackend"));
Environment env(_argc=argc, _argv=argv,
_desc = app_options,
_about=about(_name="mybackend",
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org"));

After that, you create the backend object:

// create a backend
boost::shared_ptr<Backend<double>> myBackend(backend(_name="myBackend"));
Remarks
Be careful, the backend' name has to match the name you gave at the options step.

Then, you load meshes, creates spaces etc. At solve time, or you solve with the default backend:

// solve a(u,v)=l(v)
if ( Environment::isMasterRank() )
std::cout << "With default backend\n";
a.solve(_rhs=l,_solution=u1);

One of the important backend option is to be able to monitor the residuals and iteration count

./feelpp_doc_mybackend --pc-type=id --ksp-monitor=true --myBackend.ksp-monitor=true

Finally you can create a named backend:

// solve a(u,v)=l(v)
if ( Environment::isMasterRank() )
std::cout << "With named backend\n";
a.solveb(_rhs=l,_solution=u2, _backend=myBackend);

The whole code example :

#include <feel/feel.hpp>
using namespace Feel;
int main(int argc, char**argv )
{
po::options_description app_options( "MyBackend options" );
app_options.add(backend_options("myBackend"));
Environment env(_argc=argc, _argv=argv,
_desc = app_options,
_about=about(_name="mybackend",
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org"));
// create a backend
boost::shared_ptr<Backend<double>> myBackend(backend(_name="myBackend"));
// create the mesh
auto mesh = loadMesh(_mesh=new Mesh<Simplex< 2 > > );
// function space
auto Vh = Pch<2>( mesh );
// element in Vh
auto u = Vh->element();
auto u1 = Vh->element(); //computed with default backend
auto u2 = Vh->element(); //computed with named backend
// left hand side
auto a = form2( _trial=Vh, _test=Vh );
a = integrate(_range=elements(mesh),
_expr=expr(soption("functions.alpha"))*trace(gradt(u)*trans(grad(u))) );
// right hand side
auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh),
_expr=expr(soption("functions.f"))*id(u));
// BC
a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u,
_expr=expr(soption("functions.g")));
// solve a(u,v)=l(v)
if ( Environment::isMasterRank() )
std::cout << "With default backend\n";
a.solve(_rhs=l,_solution=u1);
// solve a(u,v)=l(v)
if ( Environment::isMasterRank() )
std::cout << "With named backend\n";
a.solveb(_rhs=l,_solution=u2, _backend=myBackend);
// save results
auto e = exporter( _mesh=mesh );
e->step(0) -> add( "u", u1 );
e->step(1) -> add( "u", u2 );
e->save();
}