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>
#include <samurai/petsc.hpp> // optional, necessary for implicit schemes
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.
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
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.
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
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.
\(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:
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_size
: the field size 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' of size 2
auto u = samurai::make_field<2>("u", mesh);
// Configuration for the Laplace operator
using cfg = samurai::FluxConfig<SchemeType::LinearHomogeneous, // scheme_type
decltype(u)::size, // output_field_size (here identical to the input field size)
2, // stencil_size (for the Laplacian of order 2)
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:
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 cells captured by the stencil will be passed as argument of the flux function in the form of a cell 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& cells, ...)
{
auto& L2 = cells[0]; // {-1,0}
auto& L1 = cells[1]; // { 0,0}
auto& R1 = cells[2]; // { 1,0}
auto& R2 = cells[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 integral_constant_d)
{
static constexpr std::size_t d = decltype(integral_constant_d)::value;
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 integral_constant_d)
{
static constexpr std::size_t d = decltype(integral_constant_d)::value; // 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_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
.
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.
The role of the flux function is to return the coefficients \((c_i)_i\). Its implementation looks like
auto my_flux_function = [](double h)
{
FluxStencilCoeffs<cfg> c;
// Assuming a 2-cell stencil:
c[0] = ...
c[1] = ...
return c;
};
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_field_size x input_field_size
,
where output_field_size
is set in cfg
,
and input_field_size
corresponds to the size of the field type set in cfg
as input_field_type
.
The matrix type is in fact an xtensor
object.
You can then, amongs 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_size
and input_field_size
equal 1,
the matrix type employed to store the coefficents 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 homogenous 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
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_field<1>("u", mesh); // scalar field
using cfg = samurai::FluxConfig<SchemeType::LinearHomogeneous,
1, // output_field_size
2, // stencil_size
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
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([](double h)
{
static constexpr std::size_t L = 0; // left
static constexpr std::size_t R = 1; // right
samurai::FluxStencilCoeffs<cfg> c;
c[L] = -1/h;
c[R] = 1/h;
return c;
});
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
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 field_size
:
\((\mathbf{u}_i)_i\) are now vectors of size field_size
and \((\mathbf{c}_i)_i\) are now matrices of size field_size
\(\times\)field_size
.
For field_size = 2
, the same scheme as the scalar laplacian writes,
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 field_size = 3;
auto u = samurai::make_field<field_size>("u", mesh); // vector field
using cfg = samurai::FluxConfig<SchemeType::LinearHomogeneous,
field_size, // output_field_size
2, // stencil_size
decltype(u)>; // input_field_type
samurai::FluxDefinition<cfg> gradient([](double h)
{
static constexpr std::size_t L = 0;
static constexpr std::size_t R = 1;
samurai::FluxStencilCoeffs<cfg> c;
c[L].fill(0);
c[R].fill(0);
for (std::size_t i = 0; i < field_size; ++i)
{
c[L](i, i) = -1/h;
c[R](i, i) = 1/h;
}
return c;
});
auto laplacian = samurai::make_divergence(gradient);
Compared to the scalar laplacian code:
in the configuration, the
output_field_size
now equals thefield_size
,in the flux fonction, 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 field_size = 1
, 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 (field_size == 1)
{
c[L] = -1/h;
c[R] = 1/h;
}
else
{
c[L].fill(0);
c[R].fill(0);
for (std::size_t i = 0; i < field_size; ++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
The flux function to implement is then a discrete version of \(u\, \mathbf{n}\). Here, we choose
in all directions, i.e., in 2D
where \(B\) and \(T\) refer to bottom and top cells.
In the configuration, output_field_size
is set to the space dimension:
static constexpr std::size_t dim = decltype(mesh)::dim;
auto u = samurai::make_field<1>("u", mesh); // scalar field
using cfg = samurai::FluxConfig<SchemeType::LinearHomogeneous,
dim, // output_field_size
2, // stencil_size
decltype(u)>; // 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 = [](double h)
{
static constexpr std::size_t L = 0;
static constexpr std::size_t R = 1;
samurai::FluxStencilCoeffs<cfg> c;
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;
return c;
};
flux[y].cons_flux_function = [](double h)
{
static constexpr std::size_t B = 0;
static constexpr std::size_t T = 1;
samurai::FluxStencilCoeffs<cfg> c;
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;
return c;
};
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 integral_constant_d)
{
static constexpr std::size_t d = decltype(integral_constant_d)::value; // direction index
flux[d].cons_flux_function = [](double h)
{
static constexpr std::size_t L = 0;
static constexpr std::size_t R = 1;
samurai::FluxStencilCoeffs<cfg> c;
c[L].fill(0);
xt::row(c[L], d) = 0.5;
c[R].fill(0);
xt::row(c[R], d) = 0.5;
return c;
};
});
When dim = 1
, the matrix type in FluxStencilCoeffs<cfg>
actually reduces to a scalar type.
In order to manage all cases with one code, the content of the flux function evolves into
if constexpr (dim == 1)
{
c[L] = 0.5;
c[R] = 0.5;
}
else
{
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
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
The linear combination (1) reads
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
The following code corresponds directly to the \(n\)-dimensional version:
static constexpr std::size_t dim = decltype(mesh)::dim;
auto u = samurai::make_field<dim>("u", mesh); // vector field
using cfg = samurai::FluxConfig<SchemeType::LinearHomogeneous,
1, // output_field_size
2, // stencil_size
decltype(u)>; // input_field_type
samurai::FluxDefinition<cfg> flux;
samurai::static_for<0, dim>::apply(
[&](auto integral_constant_d)
{
static constexpr std::size_t d = decltype(integral_constant_d)::value;
flux[d].cons_flux_function = [](double)
{
static constexpr std::size_t L = 0;
static constexpr std::size_t R = 1;
samurai::FluxStencilCoeffs<cfg> c;
if constexpr (dim == 1)
{
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;
}
return c;
};
});
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_field<...>("param", mesh)
auto my_flux_function = [&](const auto& cells)
{
FluxStencilCoeffs<cfg> c;
// Assuming a 2-cell stencil:
c[0] = ... param[cell[0]] ...
c[1] = ... param[cell[1]] ...
return c;
};
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
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\),
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_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_field<1>("u", mesh); // scalar field
using cfg = samurai::FluxConfig<SchemeType::LinearHeterogeneous,
1, // output_field_size
2, // stencil_size
decltype(u)>; // input_field_type
samurai::FluxDefinition<cfg> upwind;
samurai::static_for<0, dim>::apply(
[&](auto integral_constant_d)
{
static constexpr std::size_t d = decltype(integral_constant_d)::value;
upwind[d].cons_flux_function = [&](const auto& cells)
{
static constexpr std::size_t L = 0;
static constexpr std::size_t R = 1;
auto cell = cells[L]; // arbitrary choice
samurai::FluxStencilCoeffs<cfg> c;
if (a[cell](d) >= 0)
{
c[L] = a[cell](d);
c[R] = 0;
}
else
{
c[L] = 0;
c[R] = a[cell](d);
}
return c;
};
});
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 = [](const auto& cells, const auto& field)
{
samurai::FluxValue<cfg> flux_value;
// Compute your flux using the field and the stencil cells
return flux_value;
};
Here, FluxValue<cfg>
is an array-like structure of size cfg::output_field_size
.
If cfg::output_field_size = 1
, it collapses to a simple scalar.
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 in the ghost cells for the considered field, typically with the instruction
samurai::update_ghost_mr(u);
Flux divergence#
We have
so we have to implement
where \(f_h\) is a discrete version of \(f(\cdot)\cdot \mathbf{n}\). The simple centered scheme writes
The associated code yields
samurai::FluxDefinition<cfg> f_h;
samurai::static_for<0, dim>::apply(
[&](auto integral_constant_d)
{
static constexpr std::size_t d = decltype(integral_constant_d)::value;
auto f = [](auto v)
{
samurai::FluxValue<cfg> f_v = ...;
return f_v;
};
f_h[d].cons_flux_function = [f](auto& cells, const Field& u)
{
auto& L = cells[0];
auto& R = cells[1];
return (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
Developped in 2D, where \(\mathbf{u} := [u\;v]\), it reads
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](auto& cells, const Field& u)
{
auto& L = cells[0]; // left
auto& R = cells[1]; // right
return u[L](0) >= 0 ? f_x(u[L]) : f_x(u[R]);
};
// y-direction
upwind_f[1].cons_flux_function = [f_y](auto& cells, const Field& u)
{
auto& B = cells[0]; // bottom
auto& T = cells[1]; // top
return 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
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 integral_constant_d)
{
static constexpr std::size_t d = decltype(integral_constant_d)::value;
auto f_d = [](auto u) -> samurai::FluxValue<cfg>
{
return u(d) * u;
};
upwind_f[d].cons_flux_function = [f_d](auto& cells, const Field& u)
{
auto& L = cells[0];
auto& R = cells[1];
return 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. Exemples 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
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
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 = [](auto& cells, const Field& u)
{
samurai::FluxValuePair<cfg> flux;
flux[0] = ...; // left --> right (direction '+')
flux[1] = ...; // right --> left (direction '-')
return flux;
};
Alternatively, you can also write
my_flux[0].flux_function = [](auto& cells, const Field& u) -> samurai::FluxValuePair<cfg>
{
samurai::FluxValue<cfg> fluxLR = ...; // left --> right (direction '+')
samurai::FluxValue<cfg> fluxRL = ...; // right --> left (direction '-')
return {fluxLR, fluxRL};
};
For instance, conservativity can be enforced by
my_flux[0].flux_function = [](auto& cells, const Field& u) -> samurai::FluxValuePair<cfg>
{
samurai::FluxValue<cfg> flux = ...;
return {flux, -flux};
};
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_field<samurai::DiffCoeff<dim>, 1>("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_field<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_field<1>("u", mesh);
auto unp1 = samurai::make_field<1>("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_field<dim>("velocity", mesh);
auto pressure = samurai::make_field<1>("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_field<dim>("f", mesh);
auto z = samurai::make_field<1>("z", mesh, 0.);
// Linear solver
auto stokes_solver = samurai::petsc::make_solver(stokes);
stokes_solver.set_unknowns(velocity, pressure);
stokes_solver.solve(f, z);