Finite Volume schemes#

The Finite Volume module provides a framework for the implementation of Finite Volume schemes on adapted meshes. It allows you to make use of readily-available Finite Volume operators, such as diffusion, convection and divergence operators, and implement your own flux functions.

To enable it, use

#include <samurai/schemes/fv.hpp>

In the Finite Volume method (FVM), the computed values correspond to average values, in control volumes, of the physical variables involved. The mathematical system of equations is therefore integrated locally over each control volume \(V\). Let \(u\) denote a generic, conservative data field, and \(\mathcal{D}(u)\) a differential term to be integrated over the volume \(V\). In order to ensure flux conservation, the integral over \(V\) is rewritten as a sum of fluxes passing through the faces of \(V\), i.e.

\[\int_V \mathcal{D}(u) = \int_{\partial V} \mathcal{F}(u),\]

where \(\mathcal{F}\) is the flux function associated to the differential operator \(\mathcal{D}\). For each differential operator, the expression of \(\mathcal{F}\) is given by Green’s theorem or its derived formulas. For instance, if \(\mathcal{D}(u) = \Delta u\), then \(\mathcal{F}(u) = \nabla u\cdot \mathbf{n}\), where \(\mathbf{n}\) denotes the unit vector normal to \(\partial V\) pointing out of \(V\). In the discrete setting, \(\partial V\) decomposes into a set of faces, and for each face \(F\), the flux \(\mathcal{F}(u)\) is approximated by a discrete counterpart \(\mathcal{F}_h(u_h)\). This approximation is assumed constant on each face, thus simplifying the integral into

\[\int_{\partial V} \mathcal{F}(u) = \sum_{F\subset \partial V}\int_F \mathcal{F}(u) \approx \sum_{F\subset \partial V} |F| \mathcal{F}_h(u_h)_{|F}.\]

In a nutshell, from a samurai mesh and an implementation of \(\mathcal{F}_h(u_h)\), the Finite Volume module builds a discrete differential operator that can be called like a function in an explicit context, or passed into a solver in an implicit context. Specifically, upon execution, it iterates over all cell interfaces \(F\) and computes the flux \(\mathcal{F}_h(u_h)_{|F}\) going through.

../_images/flux.svg

Denoting \(V_L\) and \(V_R\) the two cells sharing the face \(F\), and ordered in the direction of the corresponding Cartesian vector (i.e., in the x-direction, \(V_L\) and \(V_R\) are the left and right cells, respectively), the scheme contributions are defined as

\[\begin{split}\mathcal{C}_L &:= \;\;\; \frac{|F|}{|V_L|} \mathcal{F}_h(u_h)_{|F}, \\ \mathcal{C}_R &:= - \frac{|F|}{|V_R|} \mathcal{F}_h(u_h)_{|F}.\end{split}\]

The flux \(\mathcal{F}_h(u_h)_{|F}\) is computed only once and used for both contributions.

Note

Flux conservation imposes that the flux in one direction be the opposite of that in the opposite direction. Defining two fluxes \(\mathcal{F}_h^+(u_h)_{|F}\) and \(\mathcal{F}_h^-(u_h)_{|F}\) such that \(\mathcal{F}_h^-(u_h)_{|F} \neq -\mathcal{F}_h^+(u_h)_{|F}\) is also possible. Refer to non conservative schemes.

Where level jumps occur, ghost cells are involved in the computation of \(\mathcal{F}_h(u_h)_{|F}\), so that it is always computed between two cells of same length.

../_images/flux_level_jump.svg

\(V_L\) and \(V_R\) still refer to real cells in the formulas of \(\mathcal{C}_L\) and \(\mathcal{C}_R\).

Note

You can remark that the contributions are divided by the measure of the corresponding cell. Instead of \(\int_V \mathcal{D}(u)\), we actually compute \(\frac{1}{|V|} \int_V \mathcal{D}(u)\). This choice is designed to avoid multiplying the discrete source terms by the cell measures.

Implementing a Finite Volume scheme#

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

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

Static configuration#

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\). For instance, the size of \(\nabla u\), if \(u\) is a scalar field, is the space dimension.

  • the stencil_size: the size of the stencil required to compute the flux. A typical low-order scheme requires a stencil of two cells.

  • the scheme_type, to be selected amongst the values of

enum class SchemeType
{
    NonLinear,
    LinearHeterogeneous,
    LinearHomogeneous
};

This configuration must be declared in a FluxConfig static structure. Here is an example for the vectorial Laplace operator:

// Creation of a field 'u' with 2 components
auto u = samurai::make_vector_field<2>("u", mesh);

// Configuration for the Laplace operator
using cfg = samurai::FluxConfig<samurai::SchemeType::LinearHomogeneous, // scheme_type
                                2,             // stencil_size (for the Laplacian of order 2)
                                decltype(u),   // output_field_type (here identical to the input field)
                                decltype(u)>;  // input_field_type

Stencil configuration#

From cfg, you can now declare a FluxDefinition object to hold how the flux is computed.

samurai::FluxDefinition<cfg> my_flux;

In each direction, the flux is computed between one cell and its neighbour following the associated Cartesian unit vector. Namely, in 2D, from left to right and from bottom to top. The stencil is defined as an array of directional vectors from the origin cell. The origin cell is defined as the left cell in the horizontal direction, or the bottom cell in the vertical direction. It is represented by the zero vector {0,0} (or {0,0,0} in 3D). The desired neighbours are selected by inserting into the stencil the direction vectors that capture them from the origin cell. Typically, for the x-direction, the stencil {{0,0}, {1,0}} captures the origin cell and its right neighbour. Here is an example of larger stencils:

../_images/flux_stencils.svg

In this figure, four cells are required to compute the flux at the red interface. The orientation of the flux is indicated by a blue arrow. In each cell, we specify the direction vector that must be inserted into the stencil to capture it. The corresponding stencils are configured by

samurai::FluxDefinition<cfg> my_flux;

// x-direction
my_flux[0].stencil = {{-1,0}, {0,0}, {1,0}, {2,0}};
my_flux[0].cons_flux_function = my_flux_function_x;
// y-direction
my_flux[1].stencil = {{0,-1}, {0,0}, {0,1}, {0,2}};
my_flux[1].cons_flux_function = my_flux_function_y;

Note that one stencil and associated flux function must be defined for each positive Cartesian direction. The set of values captured by the stencil will be passed as argument of the flux function in the form of an array, arranged according the order chosen in the configured stencil:

my_flux[0].stencil = {{-1,0}, {0,0}, {1,0}, {2,0}};
my_flux[0].cons_flux_function = [](auto& flux_value, const auto& data, const auto& field_values)
    {
        auto& value_L2 = field_values[0]; // {-1,0}
        auto& value_L1 = field_values[1]; // { 0,0}
        auto& value_R1 = field_values[2]; // { 1,0}
        auto& value_R2 = field_values[3]; // { 2,0}
    };

Helper functions allow you to easily build line stencils (such as the one above) for any direction:

// x-direction
my_flux[0].stencil = samurai::line_stencil<dim, 0>(-1, 0, 1, 2); // {{-1,0}, {0,0}, {1,0}, {2,0}}
// y-direction
my_flux[1].stencil = samurai::line_stencil<dim, 1>(-1, 0, 1, 2); // {{0,-1}, {0,0}, {0,1}, {0,2}}

In this code, the second template parameter is the direction index. In the dynamic arguments, 0 still represents the origin cell and the other numbers its neighbours following that direction. Alternatively, instead of enumerating the neighbours, you can set only the farthest neighbour in the opposite direction along with the stencil size:

// x-direction
my_flux[0].stencil = samurai::line_stencil_from<dim, 0, 4>(-1); // {{-1,0}, {0,0}, {1,0}, {2,0}}
// y-direction
my_flux[1].stencil = samurai::line_stencil_from<dim, 1, 4>(-1); // {{0,-1}, {0,0}, {0,1}, {0,2}}

This code is strictly equivalent to the preceding one. The second template parameter is still the direction index, and the third one is the stencil size. The dynamic argument, here -1, represents the first neighbour of the sequence.

These helper functions allow you to write \(n\)-dimensional code through a static loop over the dimensions:

samurai::static_for<0, dim>::apply( // for each Cartesian direction 'd'
    [&](auto _d)
    {
        static constexpr std::size_t d = _d();

        my_flux[d].stencil = samurai::line_stencil<dim, d>(-1, 0, 1, 2);
    }

If the stencil is not specified, a the line stencil corresponding to the stencil_size stencil is used. If stencil_size is even, then the selected neighbours are evenly distributed on both sides of the interface. If it is odd, there is one more neighbour in the positive direction.

Note

If the stencil depends on the value of a dynamic parameter (e.g., in the upwind or WENO schemes, the sign of the local velocity component), then a fixed stencil composed of all cells possibly used by the scheme must be configured. The selection of the cells actually used in the computation must be performed dynamically within the flux function. See the upwind operator for an example.

Flux definition#

The flux function must be defined for each positive Cartesian direction. Here in 2D:

samurai::FluxDefinition<cfg> my_flux;
my_flux[0].cons_flux_function = my_flux_function_x; // flux in the x-direction
my_flux[1].cons_flux_function = my_flux_function_y; // flux in the y-direction

In this generic code, the flux functions remain abstract: their signatures actually depend on the SchemeType declared in cfg, and are described in the next sections. If the flux functions only differ by the direction index, you can write an \(n\)-dimensional code by using a static loop instead of the above sequence of affectations:

static constexpr std::size_t dim = 2;

samurai::FluxDefinition<cfg> my_flux;
samurai::static_for<0, dim>::apply( // for (int d=0; d<dim; d++)
        [&](auto _d)
        {
            static constexpr std::size_t d = _d(); // get the static direction index

            my_flux[d].cons_flux_function = my_flux_function_d;
        });

If the flux function is strictly identical for all directions, you can simply pass it into the constructor:

samurai::FluxDefinition<cfg> my_flux(my_flux_function);

Operator creation and usage#

Once the FluxDefinition object is constructed, you can finally declare your operator by

auto D = samurai::make_flux_based_scheme(my_flux);

and use it in an explicit context

auto v = D(u);

or in an implicit context

auto rhs = samurai::make_scalar_field<...>("rhs", mesh); // right-hand side
samurai::petsc::solve(D, u, rhs); // solves the equation D(u) = rhs

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

Note

The provided cells and associated field values correspond to the computational stencil. This is not the couple of real cells around the considered face. Where a level jump occurs, at least one of the computational cells is a ghost cell. Consequently, values must be set in the ghost cells for the considered field, typically with the instruction

samurai::update_ghost_mr(u);

This is done automatically when D(u) is called, so you do not need to call it explicitly before the operator execution. However, if the flux function captures a parameter field, its ghosts must be updated manually before the operator execution.

The definition of actual flux functions according the selected SchemeType is described in the next sections.

Linear, homogeneous operators#

This section refers to operators configured with SchemeType::LinearHomogeneous. In order to handle explicit and implicit schemes with the same definition, the user-defined flux function does not directly compute the discrete flux \(\mathcal{F}_h(u_h)\). Instead, it must return a set of coefficients. Indeed, given a face \(F\) and the associated stencil \((V_i)_i\), the discrete flux writes as a linear combination of the field values in the stencil cells, i.e.

(1)#\[\mathcal{F}_h(u_h)_{|F} := \sum_i c_i*u_i, \qquad \text{ where } u_i := (u_h)_{|V_i}.\]

The role of the flux function is to return the coefficients \((c_i)_i\). Its implementation looks like

auto my_flux_function = [](samurai::FluxStencilCoeffs<cfg>& c, double h)
{
    // Assuming a 2-cell stencil:
    c[0] = ...;
    c[1] = ...;
};

The FluxStencilCoeffs<cfg> object is an array-like structure of fixed-size. Its size is the stencil_size declared in cfg. Each c[i] is a matrix of size output_n_comp x input_n_comp, where output_n_comp is set in cfg, and input_n_comp (resp. output_n_comp) corresponds to the size of the field type set in cfg as input_field_type (resp. output_field_type). The matrix type is in fact an xtensor object. You can then, among other things, access the \(k\)-th row of c[0] via xt::row(c[0], k), its \(l\)-th column via xt::col(c[0], l), or its coefficient at indices \((k, l)\) via c[0](k, l).

Note

When both output_field_type and input_field_types are scalar fields, the matrix type employed to store the coefficients reduces to a scalar type (typically double). In particular, no accessor or xtensor function is available. To write an \(n\)-dimensional program, a separate code for the special case where the matrix reduces to a scalar is usually necessary.

As the operator is declared homogeneous over the mesh, the coefficients do not depend on specific cell values. They can, however, depend on the mesh size \(h\), making \(h\) the only parameter of the flux function. The coefficients can then be computed only once per mesh level, and re-used for all interfaces. If other (constant!) parameters are needed, they can be captured by the lambda function.

This implementation as a set of coefficients instead of the direct computation of the flux via an explicit formula allows the framework to handle both explicit and implicit schemes. Namely, the same coefficients will be used differently according to the context:

  • in an explicit context, they will be used as coefficients in the linear combination of the stencil values (formula (1)).

  • in an implicit context, they will be inserted into the matrix that will be inverted.

As examples, we implement various standard operators in the next sections.

Scalar laplacian#

Since we have

\[\int_V \Delta u = \int_{\partial V} \nabla u\cdot \mathbf{n},\]

the flux function to implement is a discrete version of \(\nabla u\cdot \mathbf{n}\). Here, we choose the normal gradient of the first order, requiring a stencil of two cells. This is enough to write the static configuration:

auto u = samurai::make_scalar_field<double>("u", mesh);

using cfg = samurai::FluxConfig<samurai::SchemeType::LinearHomogeneous,
                                2,            // stencil_size
                                decltype(u),  // output_field_type
                                decltype(u)>; // input_field_type

Now, denoting by \(V_L\) (left) and \(V_R\) (right) the stencil cells and \(F\) their interface, the discrete flux from \(V_L\) to \(V_R\) writes

\[\mathcal{F}_h(u_h)_{|F} := \frac{u_R-u_L}{h},\]

where \(u_L\) and \(u_R\) are the finite volume approximations of \(u\) in the respective cells, and \(h\) is the cell length. Referring to formula (1), the coefficients in the linear combination of \((u_L, u_R)\) correspond to \((-1/h, 1/h)\). The flux function then writes:

samurai::FluxDefinition<cfg> gradient([](samurai::FluxStencilCoeffs<cfg>& c, double h)
    {
        static constexpr std::size_t L = 0; // left
        static constexpr std::size_t R = 1; // right

        c[L] = -1/h;
        c[R] =  1/h;
    });

First of all, remark that we have declared only one flux function for all directions. We could have written as many functions as directions: they would have been identical, except that we would have replaced the name of the constants L=0, R=1 with B=0, T=1 (bottom, top) and B=0, F=1 (back, front) to better reflect the actual direction currently managed. The indexes 0 and 1 actually refer to the configured stencil. In this case, no particular stencil has been defined, so the default ones are used: in the x-direction of a 3D space, it is {{0,0,0}, {1,0,0}}, i.e. the current cell at index 0 (which we call L) and its right neighbour at index 1 (which we call R).

Finally, the operator must be constructed from the flux definition by the instruction

auto laplacian = samurai::make_flux_based_scheme(gradient);

The Laplace operator can also be viewed as the divergence of the gradient, so the following definition using the divergence function is also possible:

auto laplacian = samurai::make_divergence(gradient);

Indeed, as the divergence theorem reads

\[\int_V \nabla\cdot f(u) = \int_{\partial V} f(u) \cdot \mathbf{n},\]

the function make_divergence is simply an alias of make_flux_based_scheme, only meant to improve code readability.

Vector laplacian#

Formula (1) straightforwardly generalizes to a vectorial field \(\mathbf{u}_h\) of size n_comp:

\[\mathcal{F}_h(\mathbf{u}_h)_{|F} := \sum_i \mathbf{c}_i*\mathbf{u}_i, \qquad \text{ where } \mathbf{u}_i := (\mathbf{u}_h)_{|V_i}\]

\((\mathbf{u}_i)_i\) are now vectors of size n_comp and \((\mathbf{c}_i)_i\) are now matrices of size n_comp\(\times\)n_comp. For n_comp = 2, the same scheme as the scalar laplacian writes,

\[\begin{split}\mathcal{F}_h(\mathbf{u}_h)_{|F} := \begin{bmatrix} -1/h & 0 \\ 0 & -1/h \end{bmatrix} \mathbf{u}_L + \begin{bmatrix} 1/h & 0 \\ 0 & 1/h \end{bmatrix} \mathbf{u}_R.\end{split}\]

The matrix data structure is based on the xtensor library. You can use the xtensor syntax and access functions to manipulate matrices. The implementation of the vector laplacian operator then writes

static constexpr std::size_t n_comp = 3;
auto u = samurai::make_vector_field<n_comp>("u", mesh); // vector field

using cfg = samurai::FluxConfig<samurai::SchemeType::LinearHomogeneous,
                                2,                // stencil_size
                                decltype(u),      // output_field_type
                                decltype(u)>;     // input_field_type

samurai::FluxDefinition<cfg> gradient([](samurai::FluxStencilCoeffs<cfg>& c, double h)
    {
        static constexpr std::size_t L = 0;
        static constexpr std::size_t R = 1;

        c[L].fill(0);
        c[R].fill(0);
        for (std::size_t i = 0; i < n_comp; ++i)
        {
            c[L](i, i) = -1/h;
            c[R](i, i) =  1/h;
        }
    });

auto laplacian = samurai::make_divergence(gradient);

Compared to the scalar laplacian code, in the flux function, the coefficients for each cell of the stencil are now matrices. Specifically, they are diagonal matrices with the same coefficients as in the scalar laplacian.

When u is a scalar field, the matrix type actually reduces to a scalar type, thus forbidding instructions such as c[L](i, i). In order to manage all cases with one code, you must write

if constexpr (decltype(u)::is_scalar)
{
    c[L] = -1/h;
    c[R] =  1/h;
}
else
{
    c[L].fill(0);
    c[R].fill(0);
    for (std::size_t i = 0; i < n_comp; ++i)
    {
        c[L](i, i) = -1/h;
        c[R](i, i) =  1/h;
    }
}

Gradient#

To implement the gradient of a scalar field as Finite Volume scheme, we write

\[\int_V \nabla u = \int_{\partial V} u\, \mathbf{n}.\]

The flux function to implement is then a discrete version of \(u\, \mathbf{n}\). Here, we choose

\[\mathcal{F}_h(u_h)_{|F} := \frac{u_L+u_R}{2} \, \mathbf{n}\]

in all directions, i.e., in 2D

\[\begin{split}\mathcal{F}_h(u_h)_{|F} := \begin{cases} \begin{bmatrix} \frac{u_L+u_R}{2} \\ 0 \end{bmatrix} & \text{if } \mathbf{n} = \begin{bmatrix} 1 \\ 0 \end{bmatrix} \\ \begin{bmatrix} 0 \\ \frac{u_B+u_T}{2} \end{bmatrix} & \text{if } \mathbf{n} = \begin{bmatrix} 0 \\ 1 \end{bmatrix} \end{cases},\end{split}\]

where \(B\) and \(T\) refer to bottom and top cells.

In the configuration, the number of components of the output field is set to the space dimension:

static constexpr std::size_t dim = decltype(mesh)::dim;

auto u = samurai::make_scalar_field<double>("u", mesh);

using input_field_t                       = decltype(u);
using output_field_t = VectorField<typename input_field_t::mesh_t, typename input_field_t::value_type, dim, detail::is_soa_v<input_field_t>>;

using cfg = samurai::FluxConfig<samurai::SchemeType::LinearHomogeneous,
                                2,              // stencil_size
                                output_field_t, // output_n_comp
                                input_field_t>; // input_field_type

This time, the flux functions are different in each direction. In 2D, they write:

static constexpr std::size_t x = 0;
static constexpr std::size_t y = 1;

samurai::FluxDefinition<cfg> flux;

flux[x].cons_flux_function = [](samurai::FluxStencilCoeffs<cfg>& c, double h)
{
    static constexpr std::size_t L = 0;
    static constexpr std::size_t R = 1;

    xt::row(c[L], x) = 0.5;
    xt::row(c[L], y) = 0;

    xt::row(c[R], x) = 0.5;
    xt::row(c[R], y) = 0;
};

flux[y].cons_flux_function = [](samurai::FluxStencilCoeffs<cfg>& c, double h)
{
    static constexpr std::size_t B = 0;
    static constexpr std::size_t T = 1;

    xt::row(c[B], x) = 0;
    xt::row(c[B], y) = 0.5;

    xt::row(c[T], x) = 0;
    xt::row(c[T], y) = 0.5;
};

Here, the type FluxStencilCoeffs<cfg> contains, for each cell in the stencil, a matrix of size dim x 1, where dim is the space dimension. The flux in the x-direction computes a gradient with non-zero value only in the \(x\) coordinate; the flux in the y-direction computes a gradient with non-zero value only in the \(y\) coordinate.

This code can be compacted into the \(n\)-dimensional code

samurai::FluxDefinition<cfg> flux;
samurai::static_for<0, dim>::apply(
    [&](auto _d)
    {
        static constexpr std::size_t d = _d(); // direction index

        flux[d].cons_flux_function = [](samurai::FluxStencilCoeffs<cfg>& c, double h)
        {
            static constexpr std::size_t L = 0;
            static constexpr std::size_t R = 1;

            c[L].fill(0);
            xt::row(c[L], d) = 0.5;

            c[R].fill(0);
            xt::row(c[R], d) = 0.5;
        };
    });

Finally, we create the operator:

auto grad = make_flux_based_scheme(flux);

Divergence#

To implement the divergence of a vector field \(\mathbf{u}\) as Finite Volume scheme, we write

\[\int_V \nabla \cdot \mathbf{u} = \int_{\partial V} \mathbf{u} \cdot \mathbf{n}.\]

The flux function to implement is then a discrete version of \(\mathbf{u} \cdot \mathbf{n}\). Similarly to the gradient operator above, we choose the average value to approximate \(\mathbf{u}\). In 2D, it yields

\[\begin{split}\mathcal{F}_h(u_h)_{|F} := \begin{cases} \frac{(u_L)_x+(u_R)_x}{2} & \text{if } \mathbf{n} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \\ \frac{(u_B)_y+(u_T)_y}{2} & \text{if } \mathbf{n} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}. \end{cases}\end{split}\]

The linear combination (1) reads

\[\begin{split}\mathcal{F}_h(u_h)_{|F} := \begin{cases} \mathbf{c}_L*\mathbf{u}_L + \mathbf{c}_R*\mathbf{u}_R & \text{if } \mathbf{n} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \\ \mathbf{c}_B*\mathbf{u}_B + \mathbf{c}_T*\mathbf{u}_T & \text{if } \mathbf{n} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}. \end{cases}\end{split}\]

In this formula, \(\mathcal{F}_h(u_h)_{|F}\) is a scalar value and \(\mathbf{u}_L, \mathbf{u}_R\) (resp. \(\mathbf{u}_B, \mathbf{u}_T\)) are column-vectors. Consequently, \(\mathbf{c}_L\) and \(\mathbf{c}_R\) (resp. \(\mathbf{c}_B\) and \(\mathbf{c}_T\)) are row-vectors. Specifically, their C++ type corresponds to a matrix of size 1 x dim. The considered scheme then writes

\[\begin{split}\begin{cases} \mathbf{c}_L := \begin{bmatrix} 1/2 & 0 \end{bmatrix}, & \mathbf{c}_R := \begin{bmatrix} 1/2 & 0 \end{bmatrix}, & \text{if } \mathbf{n} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \\ \mathbf{c}_B := \begin{bmatrix} 0 & 1/2 \end{bmatrix}, & \mathbf{c}_T := \begin{bmatrix} 0 & 1/2 \end{bmatrix}, & \text{if } \mathbf{n} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}. \end{cases}\end{split}\]

The following code corresponds directly to the \(n\)-dimensional version:

static constexpr std::size_t dim = decltype(mesh)::dim;

auto u = samurai::make_vector_field<double, dim>("u", mesh);

using input_field_t  = decltype(u);
using output_field_t = ScalarField<typename input_field_t::mesh_t, typename input_field_t::value_type>;
using cfg = samurai::FluxConfig<samurai::SchemeType::LinearHomogeneous,
                                2,               // stencil_size
                                output_field_t,
                                input_field_t>;

samurai::FluxDefinition<cfg> flux;

samurai::static_for<0, dim>::apply(
    [&](auto _d)
    {
        static constexpr std::size_t d = _d();

        flux[d].cons_flux_function = [](samurai::FluxStencilCoeffs<cfg>& c, double h)
        {
            static constexpr std::size_t L = 0;
            static constexpr std::size_t R = 1;

            if constexpr (input_field_t::is_scalar)
            {
                c[L] = 0.5;
                c[R] = 0.5;
            }
            else
            {
                c[L].fill(0);
                xt::col(c[L], d) = 0.5;

                c[R].fill(0);
                xt::col(c[R], d) = 0.5;
            }
        };
    });

auto div = samurai::make_flux_based_scheme(flux);

Note that it corresponds exactly to the code of the gradient operator, where the xt:row functions have been replaced with xt:col.

Linear, heterogeneous operators#

This section refers to operators configured with SchemeType::LinearHeterogeneous. Heterogeneous operators are meant to implement linear fluxes that depend on parameters that are not constant in the domain of study. The flux function resembles that of the homogeneous linear operators, except that it receives the stencil cells as arguments instead of the sole mesh size.

auto param = samurai::make_scalar_field<double>("param", mesh);

auto my_flux_function = [&](samurai::FluxStencilCoeffs<cfg>& c, const samurai::StencilData<cfg>& data)
{
    // Assuming a 2-cell stencil:
    c[0] = ... param[data.cells[0]] ...
    c[1] = ... param[data.cells[1]] ...
};

Here, the stencil cells of the computational stencil are provided, to be used to retrieved the values of the parameters in those cells. The parameter is captured by reference by the lambda function.

Warning

The provided cells correspond to the computational stencil. This is not the couple of real cells around the considered face. Where a level jump occurs, at least one of the computational cells is a ghost cell. Consequently, make sure that values are set for your parameters in the ghost cells. For a material parameter, you can do, e.g.

samurai::for_each_cell(mesh[decltype(mesh)::mesh_id_t::reference],
    [&](auto& cell)
    {
        if (cell.center(0) < 0 && cell.center(1) > 0.5)
        {
            param[cell] = ...;
        }
        else
        {
            param[cell] = ...;
        }
    });

For a computed field used as a parameter (like, e.g., a velocity field computed from the Navier-Stokes equation), you can simply update the ghost values by

samurai::update_ghost_mr(param);

Here is an example:

Upwind convection#

Let \(u\) be a scalar field and \(\mathbf{a}\) a velocity field. We implement the operator \(\mathbf{a} \cdot \nabla u\). Assuming that \(\mathbf{a}\) is divergence-free (i.e. \(\nabla \cdot \mathbf{a} = 0\)), we write

\[\int_V \mathbf{a} \cdot \nabla u = \int_{\partial V} u(\mathbf{a} \cdot \mathbf{n}).\]

The flux function to implement is then a discrete version of \(u(\mathbf{a} \cdot \mathbf{n})\). Here, we choose the upwind scheme. In each direction \(d\),

\[\begin{split}\mathcal{F}_h(u_h)_{|F} := \begin{cases} u_L & \text{if } (\mathbf{a})_d \geq 0, \\ u_R & \text{otherwise,} \end{cases}\end{split}\]

where \((\mathbf{a})_d\) is the \(d\)-th component of \(\mathbf{a}\). Here, \(\mathbf{a}\) is constant heterogeneous. We store it in a field of size the space dimension, and set its value for each cell and ghost.

static constexpr std::size_t dim = 2;

auto a = samurai::make_vector_field<dim>("velocity", mesh);

samurai::for_each_cell(mesh[decltype(mesh)::mesh_id_t::reference],
        [](auto& cell)
        {
            if (cell.center(0) < 0)
            {
                a[cell] = {1, -1};
            }
            else
            {
                a[cell] = {-1, 1};
            }
        });

It is important to also set values into ghosts, since ghost cell may be used in the computational stencil. The construction of the operator now reads

auto u = samurai::make_scalar_field<double>("u", mesh); // scalar field

using cfg = samurai::FluxConfig<samurai::SchemeType::LinearHeterogeneous,
                                2,            // stencil_size
                                decltype(u),  // output_field_type
                                decltype(u)>; // input_field_type

samurai::FluxDefinition<cfg> upwind;

samurai::static_for<0, dim>::apply(
    [&](auto _d)
    {
        static constexpr std::size_t d = _d();

        upwind[d].cons_flux_function = [&](samurai::FluxStencilCoeffs<cfg>& c, const samurai::StencilData<cfg>& data)
        {
            static constexpr std::size_t L = 0;
            static constexpr std::size_t R = 1;

            auto cell = data.cells[L]; // arbitrary choice

            if (a[cell](d) >= 0)
            {
                c[L] = a[cell](d);
                c[R] = 0;
            }
            else
            {
                c[L] = 0;
                c[R] = a[cell](d);
            }
        };
    });

auto conv = samurai::make_flux_based_scheme(upwind);

Non-linear operators#

This section refers to operators configured with SchemeType::NonLinear. Here, the analytical formula computing the flux must be implemented. The flux function is

auto my_flux_function = [](samurai::FluxValue<cfg>& flux, const samurai::StencilData<cfg>& data, const samurai::StencilValues<cfg>& field)
{
    // Compute your flux using the stencil values of the field
    flux_value = ...;
};

Here, FluxValue<cfg> is an array-like structure of size cfg::output_n_comp. If cfg::output_n_comp = 1, it collapses to a simple scalar.

Flux divergence#

We have

\[\int_V \nabla \cdot f(u) = \int_{\partial V} f(u)\cdot \mathbf{n},\]

so we have to implement

\[\mathcal{F}_h(u_h)_{|F} := f_h(u_h),\]

where \(f_h\) is a discrete version of \(f(\cdot)\cdot \mathbf{n}\). The simple centered scheme writes

\[f_h(u_h) := \frac{f(u_L) + f(u_R)}{2}.\]

The associated code yields

samurai::FluxDefinition<cfg> f_h;

samurai::static_for<0, dim>::apply(
    [&](auto _d)
    {
        static constexpr std::size_t d = _d();

        auto f = [](auto v)
        {
            samurai::FluxValue<cfg> f_v = ...;
            return f_v;
        };

        f_h[d].cons_flux_function = [f](samurai::FluxValue<cfg>& flux, const samurai::StencilData<cfg>& data, const samurai::StencilValues<cfg>& u)
        {
            static constexpr std::size_t L = 0;
            static constexpr std::size_t R = 1;

            flux = (f(u[L]) + f(u[R])) / 2;
        };
    });

Warning

If the flux function f_h calls a continuous function f, make sure f still exists when f_h is called. Typically, if the non-linear flux operator is created and returned by a function, and if f is also defined inside that function’s scope, it will be deallocated prior to its use. In situations like this, make sure that f is captured by value by f_h. This is done in the above example by the capture instruction [f]. If f_h also requires to capture other fields or parameters by reference, you can write [f,&].

To build the operator, you can use indifferently

auto my_flux_op = samurai::make_flux_based_scheme(f_h);

or

auto my_flux_op = samurai::make_divergence(f_h);

The latter function is simply an alias of make_flux_based_scheme, proposed to improve code readability in the specific case where the operator is a flux divergence.

Convection#

We implement the operator \(\nabla \cdot \mathbf{u}\otimes\mathbf{u}\) for a vectorial field \(\mathbf{u}\). Recall that if \(\nabla \cdot \mathbf{u} = 0\), it is equivalent to \(\mathbf{u} \cdot \nabla\mathbf{u}\). We have

\[\int_V \nabla \cdot \mathbf{u}\otimes\mathbf{u} = \int_{\partial V} (\mathbf{u}\otimes\mathbf{u})\mathbf{n}.\]

Developed in 2D, where \(\mathbf{u} := [u\;v]\), it reads

\[\begin{split}(\mathbf{u}\otimes\mathbf{u})\mathbf{n} = \begin{bmatrix} u^2 & uv \\ uv & v^2 \end{bmatrix} \mathbf{n} = \begin{cases} \begin{bmatrix} u^2 \\ uv \end{bmatrix} & \text{if } \mathbf{n} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \\ \begin{bmatrix} uv \\ v^2 \end{bmatrix} & \text{if } \mathbf{n} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}. \end{cases}\end{split}\]

We implement both cases as functions:

auto f_x = [](auto u)
{
    samurai::FluxValue<cfg> f_u;
    f_u(0) = u(0) * u(0);
    f_u(1) = u(0) * u(1);
    return f_u;
};

auto f_y = [](auto u)
{
    samurai::FluxValue<cfg> f_u;
    f_u(0) = u(1) * u(0);
    f_u(1) = u(1) * u(1);
    return f_u;
};

We choose the upwind scheme, and implement:

samurai::FluxDefinition<cfg> upwind_f;

// x-direction
upwind_f[0].cons_flux_function = [f_x](samurai::FluxValue<cfg>& flux, const samurai::StencilData<cfg>& data, const samurai::StencilValues<cfg>& u)
{
    static constexpr std::size_t L = 0; // left
    static constexpr std::size_t R = 1; // right

    flux = u[L](0) >= 0 ? f_x(u[L]) : f_x(u[R]);
};

// y-direction
upwind_f[1].cons_flux_function = [f_y](samurai::FluxValue<cfg>& flux, const samurai::StencilData<cfg>& data, const samurai::StencilValues<cfg>& u)
{
    static constexpr std::size_t B = 0; // bottom
    static constexpr std::size_t T = 1; // top

    flux = u[B](1) >= 0 ? f_y(u[B]) : f_y(u[T]);
};

return samurai::make_flux_based_scheme(upwind_f);

Now, given that for each direction \(d\), we have

\[(\mathbf{u}\otimes\mathbf{u})\mathbf{n} = u_d\, \mathbf{u},\]

where \(u_d\) is the \(d\)-th component of \(\mathbf{u}\), the code can be generalized in \(n\) dimensions:

samurai::FluxDefinition<cfg> upwind_f;

samurai::static_for<0, dim>::apply(
    [&](auto _d)
    {
        static constexpr std::size_t d = _d();

        auto f_d = [](auto u) -> samurai::FluxValue<cfg>
        {
            return u(d) * u;
        };

        upwind_f[d].cons_flux_function = [f_d](samurai::FluxValue<cfg>& flux, const samurai::StencilData<cfg>& data, const samurai::StencilValues<cfg>& u)
        {
            static constexpr std::size_t L = 0;
            static constexpr std::size_t R = 1;

            flux = u[L](d) >= 0 ? f_d(u[L]) : f_d(u[R]);
        };
    });

return samurai::make_flux_based_scheme(upwind_f);

Implementing a non-conservative scheme#

Flux-based, non-conservative schemes also exist. Examples can be found in two-phase flow simulation: while the scheme remains conservative within each phase, non-conservative fluxes can be computed at the interface between phases.

We recall \(V_L\) and \(V_R\), the two cells sharing the face \(F\), and ordered in the direction of the corresponding Cartesian vector (i.e., in the x-direction, \(V_L\) and \(V_R\) are the left and right cells, respectively). In a conservative scheme, the respective contributions of \(F\) on \(V_L\) and \(V_R\) are defined as

\[\begin{split}\mathcal{C}_L &:= \;\;\; \frac{|F|}{|V_L|} \mathcal{F}_h(u_h)_{|F}, \\ \mathcal{C}_R &:= - \frac{|F|}{|V_R|} \mathcal{F}_h(u_h)_{|F},\end{split}\]

and flux \(\mathcal{F}_h(u_h)_{|F}\) is computed only once and used for both contributions.

Now, the contributions of a non-conservative scheme read

\[\begin{split}\mathcal{C}_L &:= \frac{|F|}{|V_L|} \mathcal{F}_h^+(u_h)_{|F}, \\ \mathcal{C}_R &:= \frac{|F|}{|V_R|} \mathcal{F}_h^-(u_h)_{|F},\end{split}\]

where we do not necessarily have \(\mathcal{F}_h^-(u_h)_{|F} = -\mathcal{F}_h^+(u_h)_{|F}\).

Implementation-wise, while conservative schemes implement \(\mathcal{F}_h(u_h)_{|F}\) through cons_flux_function, non-conservative ones return both values of \(\mathcal{F}_h^+(u_h)_{|F}\) and \(\mathcal{F}_h^-(u_h)_{|F}\) through flux_function. So far, only non-linear schemes are possible. The signature is the same as flux_function, except that it returns two values instead of one:

samurai::FluxDefinition<cfg> my_flux;

my_flux[0].flux_function = [](samurai::FluxValuePair<cfg>& flux, const samurai::StencilData<cfg>& data, const samurai::StencilValues<cfg>& u)
                           {
                               flux[0] = ...; // left --> right (direction '+')
                               flux[1] = ...; // right --> left (direction '-')
                           };

For instance, conservativity can be enforced by

my_flux[0].flux_function = [](samurai::FluxValuePair<cfg>& flux, const samurai::StencilData<cfg>& data, const samurai::StencilValues<cfg>& u)
                           {
                               flux[0] = ...;
                               flux[1] = -flux[0];
                           };

Available implementations#

The following operators and schemes are available in the framework:

Diffusion#

The diffusion operator \(-\nabla\cdot(K\nabla u)\) is implemented with the classical scheme of order 2. It applies to both scalar and vector fields. Here are some examples according to the type of coefficient:

  • Homogeneous, scalar coefficient:

double K = 10;
auto diff = samurai::make_diffusion_order2<FieldType>(K);

  • Homogeneous, diagonal tensor coefficient:

static constexpr std::size_t dim = 2;

samurai::DiffCoeff<dim> K;
K(0) = 1; // x-direction
K(1) = 2; // y-direction

auto diff = samurai::make_diffusion_order2<FieldType>(K);

  • Heterogeneous, diagonal tensor coefficient:

static constexpr std::size_t dim = 2;

auto K = samurai::make_scalar_field<samurai::DiffCoeff<dim>>("K", mesh);

samurai::for_each_cell(mesh[decltype(mesh)::mesh_id_t::reference],
        [&](auto& cell)
        {
            double x = cell.center(0);
            K[cell] = x < 0 ? {1, 2} : {2, 1};
        });

auto diff = samurai::make_diffusion_order2<FieldType>(K);

  • Laplacian operator (equivalent to the diffusion operator with no minus sign and constant coefficient 1):

auto lap = samurai::make_laplacian_order2<FieldType>();

Convection#

The convection operators are accessible via the function make_convection_SCHEME<FieldType>(...), where SCHEME must be replaced with the name of desired the discrete scheme. Two discrete schemes are implemented: upwind and weno5 (Jiang & Shu). The mathematical operator implemented is \(\nabla \cdot (a \otimes u)\), which corresponds to \(a\cdot\nabla u\) if \(a\) is divergence-free.

  • Linear convection with constant velocity:

samurai::VelocityVector<dim> a = {1, -2};
auto conv = samurai::make_convection_SCHEME<FieldType>(a);

  • Linear convection with a velocity field:

auto a = samurai::make_vector_field<double, dim>("velocity", mesh); // the size must correspond to the space dimension
auto conv = samurai::make_convection_SCHEME<FieldType>(a);

  • Non-linear convection \(\nabla \cdot (u \otimes u)\):

auto conv = samurai::make_convection_SCHEME<FieldType>();

Gradient#

The gradient operator is only implemented for scalar fields with the centered scheme of order 2.

auto grad = samurai::make_gradient_order2<FieldType>();

Divergence#

The gradient operator is only implemented for vector fields of size the space dimension, with the centered scheme of order 2.

auto div = samurai::make_divergence_order2<FieldType>();

Identity#

The identity operator is only used for the implementation of implicit time stepping schemes. Example of the heat equation

auto u    = samurai::make_scalar_field<double>("u", mesh);
auto unp1 = samurai::make_scalar_field<double>("unp1", mesh);

auto diff = samurai::make_diffusion_order2<decltype(u)>();
auto id   = samurai::make_identity<decltype(u)>();

samurai::petsc::solve(id + dt*diff, unp1, u); // solves the linear equation [id + dt*diff](unp1) = u

Zero operator#

The zero operator is used to build block matrices where one of the blocks is zero. Example of the Stokes system:

// Unknowns
auto velocity = samurai::make_vector_field<double, dim>("velocity", mesh);
auto pressure = samurai::make_scalar_field<double>("pressure", mesh);

// Stokes operator
auto diff = samurai::make_diffusion_order2<decltype(velocity)>();
auto grad = samurai::make_gradient_order2<decltype(pressure)>();
auto div  = samurai::make_divergence_order2<decltype(velocity)>();
auto zero = samurai::make_zero_operator<decltype(pressure)>();

auto stokes = samurai::make_block_operator<2, 2>(diff, grad,
                                                 -div, zero);
// Right-hand side
auto f = samurai::make_vector_field<double, dim>("f", mesh);
auto z = samurai::make_scalar_field<double>("z", mesh, 0.);

// Linear solver
auto stokes_solver = samurai::petsc::make_solver(stokes);
stokes_solver.set_unknowns(velocity, pressure);
stokes_solver.solve(f, z);