Local schemes#

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

#include <samurai/schemes/fv.hpp>

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_type: the C++ type 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_scalar_field<double>("u", mesh);

using cfg  = samurai::LocalCellSchemeConfig<
                    SchemeType::NonLinear, // scheme_type
                    decltype(u),           // output_field_type (here identical to the input field)
                    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_scalar_field<double>("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([&](samurai::SchemeValue<cfg>& result, const auto& cell, const auto& field)
{
    // Local field value
    auto v = field[cell];

    // Use 'v' and captured parameters in your computation
    result = ...;
});

The parameters of the function are

  • result: the result of the operator application, to be filled in the function;

  • 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_fied_type::n_comp (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_type is a scalar field, 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([&](samurai::JacobianMatrix<cfg>& jac, const auto& cell, const auto& field)
{
    // Local field value
    auto v = field[cell];

    jac = ...; // Fill the jacobian matrix
});

Warning

The type JacobianMatrix<cfg> is a matrix of size output_n_comp x input_field_type. However, if both output_field_type and input_field_type are scalar fields, it reduces to a scalar type (typically double).