Diffusion equation model definition
In the file model.hh
, we define the model equations and
set all default model properties. The setup consist of three steps:
- Create a model type tag (used to specialize properties)
- Define the local residual class implementing the discrete equation
- Specialize important properties of the model such that Dumux knows how to assemble the system matrix
We start model.hh
with the necessary header includes:
Click to show includes
#include <dumux/common/math.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/common/volumevariables.hh>
#include <dumux/discretization/method.hh>
1. Property Tag
The property tag is simply an empty struct with the name DiffusionModel
Click to hide/show the file documentation (or inspect the source code)
namespace Dumux::Properties::TTag {
//! The diffusion model tag that we can specialize properties for
struct DiffusionModel {};
} // end namespace Dumux::Properties::TTag
2. The local (element-wise) residual
The local residual assembles the contribution to the residual for
all degrees of freedom associated with an element. Here, we use the
Box method which is based on P_1 basis functions (piece-wise linears)
and the degrees of freedom are on the nodes. Each node is associate with
exactly one sub control volume (scv
) per element and several (2 in \mathbb{R}^2)
sub control volume faces (scvf
). In the local residual, we can implement the
contribution for one scv
(storage and source terms) or one scvf
(flux terms).
Let's have a look at the class implementation.
Click to hide/show the file documentation (or inspect the source code)
The class DiffusionModelLocalResidual
inherits from something called BaseLocalResidual
.
This base class differs depending on the chosen discretization scheme. For the box method
(which is a control-volume finite element scheme) used in this example, the property
BaseLocalResidual
is specialized to CVFELocalResidual<TypeTag>
in dumux/discretization/box.hh.
Since this local residual only works with control-volume finite element schemes due to
the flux implementation, we could have also chosen to inherit from public CVFELocalResidual<TypeTag>
.
namespace Dumux {
template<class TypeTag>
class DiffusionModelLocalResidual
: public GetPropType<TypeTag, Properties::BaseLocalResidual>
{
Storage term: Evaluate the rate of change of all conserved quantities
NumEqVector computeStorage(const Problem& problem,
const SubControlVolume& scv,
const VolumeVariables& volVars) const
{
NumEqVector storage;
storage[Indices::massBalanceEqIdx] = volVars.priVar(Indices::concentrationIdx);
return storage;
}
Flux term: Evaluate the fluxes over a face of a sub control volume Here we evaluate the (integrated) flux
F = -D \sum_{B \in \mathcal{B}_K} \left( c_{h,B} \nabla N_B \right) \cdot\boldsymbol{n} \vert \sigma \vert
NumEqVector computeFlux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
static_assert(DiscretizationMethods::isCVFE<typename GridGeometry::DiscretizationMethod>,
"This local residual is hard-coded to control-volume finite element schemes");
// Compute ∇c at the integration point of this sub control volume face.
const auto& fluxVarCache = elemFluxVarsCache[scvf];
Dune::FieldVector<Scalar, dimWorld> gradConcentration(0.0);
for (const auto& scv : scvs(fvGeometry))
{
const auto& volVars = elemVolVars[scv];
// v.axpy(a, w) means v <- v + a*w
gradConcentration.axpy(
volVars.priVar(Indices::concentrationIdx),
fluxVarCache.gradN(scv.indexInElement())
);
}
NumEqVector flux;
// Compute the flux with `vtmv` (vector transposed times matrix times vector) or -n^T D ∇c A.
// The diffusion coefficient comes from the `problem` (see Part II of the example).
flux[Indices::massBalanceEqIdx] = -1.0*vtmv(
scvf.unitOuterNormal(), problem.diffusionCoefficient(), gradConcentration
)*scvf.area();
return flux;
}
};
} // end namespace Dumux
3. The model properties
By specializing properties for our type tag DiffusionModel
,
every other class that knows about the type tag (this will be
for example the assembler or the problem), can extract the
type information the we specify here.
Note that these types can be overwritten for specific problem definitions if this is needed (we will show this on the next page).
Click to hide/show the file documentation (or inspect the source code)
namespace Dumux::Properties {
The type of the local residual is the class defined above
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::DiffusionModel>
{ using type = DiffusionModelLocalResidual<TypeTag>; };
The default scalar type is double we compute with double precision floating point numbers
template<class TypeTag>
struct Scalar<TypeTag, TTag::DiffusionModel>
{ using type = double; };
The model traits specify some information about our equation system Here we have just one equation. We still specify indices so in the places where we access primary variables, we can do so with a named variable.
template<class TypeTag>
struct ModelTraits<TypeTag, TTag::DiffusionModel>
{
struct type
{
struct Indices
{
static constexpr int concentrationIdx = 0;
static constexpr int massBalanceEqIdx = 0;
};
static constexpr int numEq() { return 1; }
};
};
The primary variable vector has entries of type Scalar
and is
as large as the number of equations (here 1) but we keep it general.
template<class TypeTag>
struct PrimaryVariables<TypeTag, TTag::DiffusionModel>
{
using type = Dune::FieldVector<
GetPropType<TypeTag, Properties::Scalar>,
GetPropType<TypeTag, Properties::ModelTraits>::numEq()
>;
};
The BasicVolumeVariables
are the simples volume variables
They only store one instance of PrimaryVariables
for the
degree of freedom (here: vertex dof) that they are attached to.
template<class TypeTag>
struct VolumeVariables<TypeTag, TTag::DiffusionModel>
{
struct Traits
{
using PrimaryVariables
= GetPropType<TypeTag, Properties::PrimaryVariables>;
};
using type = BasicVolumeVariables<Traits>;
};
} // end namespace Dumux::Properties