Local schemes

The local schemes are part of the Finite Volume module, enabled by

#include <samurai/schemes/fv.hpp>
#include <samurai/petsc.hpp> // optional, necessary for implicit shemes

They are characterized by a function that applies to a field, whose computation only involves information located in the current mesh cell. The C++ interface is built similarly to the flux-based Finite Volume schemes.

Implementing local scheme

Consider the implementation of a local operator \(\mathcal{A}\) that applies to a field \(u\). Let \(v\) denote the resulting field:

\[v = \mathcal{A}(u).\]

First of all, the structural information about the scheme must be declared. It contains:

• the input_field_type: the C++ type of the field \(u\).

• the output_field_size: the field size of the resulting field \(v\).

• the scheme_type, to be selected amongst the values of

enum class SchemeType
{
    NonLinear,
    LinearHeterogeneous,
    LinearHomogeneous
};

This configuration must be declared in a LocalCellSchemeConfig static structure. Here is an example:

auto u = samurai::make_field<...>("u", mesh);

using cfg  = samurai::LocalCellSchemeConfig<
                    SchemeType::NonLinear, // scheme_type
                    decltype(u)::size,     // output_field_size (here identical to the input field size)
                    decltype(u)>;          // input_field_type

Secondly, we create the operator from the configuration cfg:

auto A = samurai::make_cell_based_scheme<cfg>();

A.set_scheme_function(...);
A.set_jacobian_function(...); // only A is non-linear

The signature of the scheme function actually depends on the SchemeType declared in cfg (see sections below).

Once the operator is created and defined, it can be used in an explicit context

auto v = A(u);

or in an implicit context

auto b = samurai::make_field<...>("b", mesh);
samurai::petsc::solve(A, u, b); // solves the equation A(u) = b

Note that the solve function involves a linear or a non-linear solver according to the SchemeType declared in cfg.

Non-linear operators

The analytical formula of the operator is implemented as a lambda function.

A.set_scheme_function([&](auto& cell, const auto& field)
{
    // Local field value
    auto v = field[cell];

    // Use 'v' and captured parameters in your computation
    samurai::SchemeValue<cfg> result = ...;

    return result;
});

The parameters of the function are

• cell: the current local cell;

• field: the input field, to which the operator applies. Its actual type is declared in cfg.

The return type SchemeValue<cfg> is a array-like structure of size output_field_size (declared in cfg). It is based on the xtensor library, so all xtensor functions and accessors can be used. The \(i\)-th component can be accessed with result(i).

Note

If output_field_size is set to 1, SchemeValue<cfg> reduces to a scalar type (typically double).

If the operator is to be implicited, its jacobian function must also be defined. If only explicit applications of the operator shall be used, then this step is optional.

A.set_jacobian_function([&](auto& cell, const auto& field)
{
    // Local field value
    auto v = field[cell];

    samurai::JacobianMatrix<cfg> jac = ...
    return jac;
});

Warning

The type JacobianMatrix<cfg> is a matrix of size output_field_size x input_field_type. However, if output_field_size = input_field_size = 1, it reduces to a scalar type (typically double).