Using a backend

1. 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 at 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 in size \(A_1 y= c\) with the same backend or the default without rebuilding it, it will fail.

backend(_rebuild=true)->solve(_matrix=A1,_rhs=c,_sol=y);

Each of that options can be retrieved via the --help-lib argument in the command line.

1.1. 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
	auto myBackend = backend(_name="myBackend");

the backend’s 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)
	Feel::cout << "With default backend\n";
	a.solve(_rhs=l,_solution=u1); // Compute with default backend

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

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

Finally you can create a named backend:

	// solve a(u,v)=l(v)
	Feel::cout << "With named backend\n";
	a.solveb(_rhs=l,_solution=u2, _backend=myBackend); // Compute with myBackend

2. Implementation

The implementation reads as follows:

#include <feel/feelcore/environment.hpp>
#include <feel/feeldiscr/pch.hpp>
#include <feel/feelfilters/unitsquare.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelvf/vf.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
	auto myBackend = backend(_name="myBackend");

	// create the mesh
	auto mesh = unitSquare();

	// function space
	auto Vh = Pch<2>( mesh );

	// element in Vh
	auto u  = Vh->element();
	auto u1 = Vh->element();
	auto u2 = Vh->element();

	// 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)
	Feel::cout << "With default backend\n";
	a.solve(_rhs=l,_solution=u1); // Compute with default backend
	// solve a(u,v)=l(v)
	Feel::cout << "With named backend\n";
	a.solveb(_rhs=l,_solution=u2, _backend=myBackend); // Compute with myBackend

	// save results
	auto e = exporter( _mesh=mesh );
	e->step(0) -> add( "u", u1 );
	e->step(1) -> add( "u", u2 );
	e->save();
}

and the associated config file

ksp-monitor=true
ksp-rtol=1e-5
backend.verbose=true
pc-type=id

[functions]
alpha=1
f=x+y:x:y
g=sin(x):x

[myBackend]
backend.verbose=true
ksp-monitor=true
ksp-rtol=1e-5
pc-type=lu