From 5cb1de0a5e1b368f31a94422b71c0f449eb4add3 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" <dennis.glaeser@iws.uni-stuttgart.de> Date: Tue, 24 Mar 2020 13:23:22 +0100 Subject: [PATCH] [example][1ptracer] edit docu --- examples/1ptracer/.doc_config | 2 + examples/1ptracer/README.md | 752 +++++++++++-------------- examples/1ptracer/doc/intro.md | 8 +- examples/1ptracer/main.cc | 214 +++---- examples/1ptracer/problem_1p.hh | 89 +-- examples/1ptracer/problem_tracer.hh | 130 +---- examples/1ptracer/properties_1p.hh | 112 ++++ examples/1ptracer/properties_tracer.hh | 156 +++++ examples/1ptracer/spatialparams_1p.hh | 53 +- 9 files changed, 770 insertions(+), 746 deletions(-) create mode 100644 examples/1ptracer/properties_1p.hh create mode 100644 examples/1ptracer/properties_tracer.hh diff --git a/examples/1ptracer/.doc_config b/examples/1ptracer/.doc_config index 028f23fcd6..8f6f18e059 100644 --- a/examples/1ptracer/.doc_config +++ b/examples/1ptracer/.doc_config @@ -3,8 +3,10 @@ "doc/intro.md", "spatialparams_1p.hh", "problem_1p.hh", + "properties_1p.hh", "spatialparams_tracer.hh", "problem_tracer.hh", + "properties_tracer.hh", "main.cc", "doc/results.md" ] diff --git a/examples/1ptracer/README.md b/examples/1ptracer/README.md index 02bb507e6c..b7fc9700b7 100644 --- a/examples/1ptracer/README.md +++ b/examples/1ptracer/README.md @@ -23,22 +23,22 @@ with the darcy velocity $` \textbf v `$, the permeability $` \textbf K`$, the dy Darcy's law is inserted into the mass balance equation: ```math -\phi \frac{\partial \varrho}{\partial t} + \text{div} \textbf v = 0 +\phi \frac{\partial \varrho}{\partial t} + \text{div} \textbf v = 0, ``` where $`\phi`$ is the porosity. -The equation is discretized using a cell-centered finite volume scheme as spatial discretization for the pressure as primary variable. For details on the discretization scheme, have a look at the dumux [handbook](https://dumux.org/handbook). +The equation is discretized using cell-centered finite volumes with two-point flux approximation as spatial discretization scheme for the pressure as primary variable. For details on the discretization schemes available in DuMuX, have a look at the [handbook](https://dumux.org/handbook). ### Tracer Model -The transport of the contaminant component $`\kappa`$ is based on the previously evaluated velocity field $`\textbf v`$ with the help of the following mass balance equation: +The transport of the contaminant component $`\kappa`$ occurs with the velocity field $`\textbf v`$ +that is computed with the __1p model__ (see above): ```math \phi \frac{ \partial \varrho X^\kappa}{\partial t} - \text{div} \left\lbrace \varrho X^\kappa {\textbf v} + \varrho D^\kappa_\text{pm} \textbf{grad} X^\kappa \right\rbrace = 0, ``` where $`X^\kappa`$ is the mass fraction of the contaminant component $`\kappa`$ and $` D^\kappa_\text{pm} `$ is the effective diffusivity. - The effective diffusivity is a function of the diffusion coefficient of the component $`D^\kappa`$ and the porosity and tortuosity $`\tau`$ of the porous medium (see [dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh)): ```math @@ -54,30 +54,26 @@ For the tracer model, this is done in the files `problem_tracer.hh` and `spatial ## The file `spatialparams_1p.hh` -In this file, we generate a random permeability field in the constructor of the `OnePTestSpatialParams` class. -For this, we use the random number generation facilities provided by the C++ standard library. +This file contains the __spatial parameters class__ which defines the +distributions for the porous medium parameters permeability and porosity +over the computational grid +In this example, we use a randomly generated and element-wise distributed +permeability field. For this, we use the random number generation facilitie +provided by the C++ standard library. ```cpp #include <random> ``` - -We use the properties for porous medium flow models, declared in the file `properties.hh`. - -```cpp -#include <dumux/porousmediumflow/properties.hh> -``` - -We include the spatial parameters class for single-phase models discretized by finite volume schemes. -The spatial parameters defined for this example will inherit from those. - +We include the spatial parameters class for single-phase models discretized +by finite volume schemes, from which the spatial parameters defined for this +example will inherit. ```cpp #include <dumux/material/spatialparams/fv1p.hh> namespace Dumux { ``` - -In the `OnePTestSpatialParams` class, we define all functions needed to describe the porous medium, e.g. porosity and permeability for the 1p_problem. - +In the `OnePTestSpatialParams` class, we define all functions needed to describe +the porous medium, e.g. porosity and permeability, for the 1p_problem. ```cpp template<class GridGeometry, class Scalar> class OnePTestSpatialParams @@ -85,9 +81,7 @@ class OnePTestSpatialParams OnePTestSpatialParams<GridGeometry, Scalar>> { ``` - -We declare aliases for types that we are going to need in this class. - +The following convenience aliases will be used throughout this class: ```cpp using GridView = typename GridGeometry::GridView; using FVElementGeometry = typename GridGeometry::LocalView; @@ -107,21 +101,20 @@ Here, we are using scalar permeabilities, but tensors are also supported. ```cpp using PermeabilityType = Scalar; +``` +### Generation of the random permeability field +We generate the random permeability field upon construction of the spatial parameters class +```cpp OnePTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry) : ParentType(gridGeometry), K_(gridGeometry->gridView().size(0), 0.0) { ``` - -### Generation of the random permeability field -We get the permeability of the domain and the lens from the `params.input` file. - +The permeability of the domain and the lens are obtained from the `params.input` file. ```cpp permeability_ = getParam<Scalar>("SpatialParams.Permeability"); permeabilityLens_ = getParam<Scalar>("SpatialParams.PermeabilityLens"); ``` - -Furthermore, we get the position of the lens, which is defined by the position of the lower left and the upper right corner. - +Furthermore, the position of the lens, which is defined by the position of the lower left and the upper right corners, are obtained from the input file. ```cpp lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft"); lensUpperRight_ =getParam<GlobalPosition>("SpatialParams.LensUpperRight"); @@ -154,14 +147,14 @@ Here, we use element-wise distributed permeabilities that were randomly generate const SubControlVolume& scv, const ElementSolution& elemSol) const { - return K_[scv.dofIndex()]; + return K_[scv.elementIndex()]; } ``` We set the porosity $`[-]`$ for the whole domain to a value of $`20 \%`$. ```cpp - Scalar porosityAtPos(const GlobalPosition &globalPos) const + Scalar porosityAtPos(const GlobalPosition& globalPos) const { return 0.2; } ``` @@ -173,27 +166,24 @@ We reference to the permeability field. This is used in the main function to wri private: ``` - -We have a convenient definition of the position of the lens. - +The following function returns true if a given position is inside the lens. +We use an epsilon of 1.5e-7 here for floating point comparisons. ```cpp - bool isInLens_(const GlobalPosition &globalPos) const + bool isInLens_(const GlobalPosition& globalPos) const { - for (int i = 0; i < dimWorld; ++i) { - if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_) + for (int i = 0; i < dimWorld; ++i) + { + if (globalPos[i] < lensLowerLeft_[i] + 1.5e-7 + || globalPos[i] > lensUpperRight_[i] - 1.5e-7) return false; } + return true; } - GlobalPosition lensLowerLeft_; - GlobalPosition lensUpperRight_; - + GlobalPosition lensLowerLeft_, lensUpperRight_; Scalar permeability_, permeabilityLens_; - std::vector<Scalar> K_; - - static constexpr Scalar eps_ = 1.5e-7; }; } // end namespace Dumux @@ -205,177 +195,25 @@ We have a convenient definition of the position of the lens. ## The file `problem_1p.hh` -Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties. -### Include files -We use `YaspGrid`, an implementation of the dune grid interface for structured grids: - -```cpp -#include <dune/grid/yaspgrid.hh> -``` - -The cell centered, two-point-flux discretization scheme is included: - -```cpp -#include <dumux/discretization/cctpfa.hh> -``` - -The one-phase flow model is included: - -```cpp -#include <dumux/porousmediumflow/1p/model.hh> -``` - -This is the porous medium problem class that this class is derived from: +This file contains the __problem class__ which defines the initial and boundary +conditions for the single-phase flow simulation. +### Include files +This header contains the porous medium problem class that this class is derived from: ```cpp #include <dumux/porousmediumflow/problem.hh> ``` - -The fluid properties are specified in the following headers (we use liquid water as the fluid phase): - -```cpp -#include <dumux/material/components/simpleh2o.hh> -#include <dumux/material/fluidsystems/1pliquid.hh> -``` - -The local residual for incompressible flow is included: - -```cpp -#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh> -``` - -We include the header that specifies all spatially variable parameters: - +This header contains the class that specifies all spatially variable parameters +related to this problem. ```cpp #include "spatialparams_1p.hh" ``` - -### Define basic properties for our simulation -We enter the namespace Dumux in order to import the entire Dumux namespace for general use - -```cpp -namespace Dumux { -``` - -The problem class is forward declared: - -```cpp -template<class TypeTag> -class OnePTestProblem; -``` - -We enter the namespace Properties, which is a sub-namespace of the namespace Dumux: - -```cpp -namespace Properties { -``` - -A `TypeTag` for our simulation is created which inherits from the one-phase flow model and the cell centered, two-point-flux discretization scheme. - -```cpp -namespace TTag { -struct IncompressibleTest { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; }; -} -``` - -We use a structured 2D grid: - -```cpp -template<class TypeTag> -struct Grid<TypeTag, TTag::IncompressibleTest> { using type = Dune::YaspGrid<2>; }; -``` - -The problem class specifies initial and boundary conditions: - -```cpp -template<class TypeTag> -struct Problem<TypeTag, TTag::IncompressibleTest> { using type = OnePTestProblem<TypeTag>; }; -``` - -We define the spatial parameters for our simulation: - -```cpp -template<class TypeTag> -struct SpatialParams<TypeTag, TTag::IncompressibleTest> -{ -``` - -We define convenient shortcuts to the properties `GridGeometry` and `Scalar`: - -```cpp - using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; - using Scalar = GetPropType<TypeTag, Properties::Scalar>; -``` - -Finally, we set the spatial parameters: - -```cpp - using type = OnePTestSpatialParams<GridGeometry, Scalar>; -}; -``` - -We use the local residual that contains analytic derivative methods for incompressible flow: - -```cpp -template<class TypeTag> -struct LocalResidual<TypeTag, TTag::IncompressibleTest> { using type = OnePIncompressibleLocalResidual<TypeTag>; }; -``` - -In the following we define the fluid properties. - -```cpp -template<class TypeTag> -struct FluidSystem<TypeTag, TTag::IncompressibleTest> -{ -``` - -We define a convenient shortcut to the property Scalar: - -```cpp - using Scalar = GetPropType<TypeTag, Properties::Scalar>; -``` - -We create a fluid system that consists of one liquid water phase. We use the simple -description of water, which means we do not use tabulated values but more general equations of state. - -```cpp - using type = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >; -}; -``` - -We enable caching for the grid volume variables. - -```cpp -template<class TypeTag> -struct EnableGridVolumeVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; }; -``` - -We enable caching for the grid flux variables. - -```cpp -template<class TypeTag> -struct EnableGridFluxVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; }; -``` - -We enable caching for the FV grid geometry - -```cpp -template<class TypeTag> -struct EnableGridGeometryCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; }; -``` - -The cache stores values that were already calculated for later usage. This increases the memory demand but makes the simulation faster. -We leave the namespace Properties. - -```cpp -} -``` - ### The problem class We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation. -As this is a porous medium problem, we inherit from the base class `PorousMediumFlowProblem`. - +As this is a porous medium flow problem, we inherit from the base class `PorousMediumFlowProblem`. ```cpp +namespace Dumux { + template<class TypeTag> class OnePTestProblem : public PorousMediumFlowProblem<TypeTag> { @@ -496,6 +334,123 @@ We leave the namespace Dumux. +## The file `properties_1p.hh` + + +This file defines the `TypeTag` used for the single-phase simulation, for +which we then define the necessary properties. + +### Include files +The `TypeTag` defined for this simulation will inherit all properties from the +`OneP` type tag, a convenience type tag that predefines most of the required +properties for single-phase flow simulations in DuMuX. The properties that are +defined in this file are those that depend on user choices and no meaningful +default can be set. +```cpp +#include <dumux/porousmediumflow/1p/model.hh> +``` +We want to use `YaspGrid`, an implementation of the dune grid interface for structured grids: +```cpp +#include <dune/grid/yaspgrid.hh> +``` +In this example, we want to discretize the equations with the cell centered finite volume +scheme using two-point-flux approximation: +```cpp +#include <dumux/discretization/cctpfa.hh> +``` +The fluid properties are specified in the following headers (we use liquid water as the fluid phase): +```cpp +#include <dumux/material/components/simpleh2o.hh> +#include <dumux/material/fluidsystems/1pliquid.hh> +``` +The local residual for incompressible flow is included. +The one-phase flow model (included above) uses a default implementation of the +local residual for single-phase flow. However, in this example we are using an +incompressible fluid phase. Therefore, we are including the specialized local +residual which contains functionality to analytically compute the entries of +the Jacobian matrix. We will use this in the main file. +```cpp +#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh> +``` +We include the problem and spatial parameters headers used for this simulation. +```cpp +#include "problem_1p.hh" +#include "spatialparams_1p.hh" +``` +### Basic property definitions for the 1p problem +We enter the namespace Dumux::Properties in order to import the entire Dumux namespace for general use: +```cpp +namespace Dumux:: Properties { +``` +A `TypeTag` for our simulation is created which inherits from the one-phase flow model +and the cell centered finite volume scheme with two-point-flux discretization scheme: +```cpp +namespace TTag { +struct IncompressibleTest { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; }; +} +``` +We use a structured 2D grid: +```cpp +template<class TypeTag> +struct Grid<TypeTag, TTag::IncompressibleTest> { using type = Dune::YaspGrid<2>; }; +``` +The problem class specifies initial and boundary conditions: +```cpp +template<class TypeTag> +struct Problem<TypeTag, TTag::IncompressibleTest> { using type = OnePTestProblem<TypeTag>; }; +``` +We define the spatial parameters for our simulation: +```cpp +template<class TypeTag> +struct SpatialParams<TypeTag, TTag::IncompressibleTest> +{ +private: + using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; + using Scalar = GetPropType<TypeTag, Properties::Scalar>; +public: + using type = OnePTestSpatialParams<GridGeometry, Scalar>; +}; +``` +We use the local residual that contains analytic derivative methods for incompressible flow: +```cpp +template<class TypeTag> +struct LocalResidual<TypeTag, TTag::IncompressibleTest> { using type = OnePIncompressibleLocalResidual<TypeTag>; }; +``` +In the following we define the fluid system to be used: +```cpp +template<class TypeTag> +struct FluidSystem<TypeTag, TTag::IncompressibleTest> +{ +private: + using Scalar = GetPropType<TypeTag, Properties::Scalar>; +public: + using type = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >; +}; +``` +This enables grid-wide caching of the volume variables. +```cpp +template<class TypeTag> +struct EnableGridVolumeVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; }; +``` +This enables grid wide caching for the flux variables. +```cpp +template<class TypeTag> +struct EnableGridFluxVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; }; +``` +This enables grid-wide caching for the finite volume grid geometry +```cpp +template<class TypeTag> +struct EnableGridGeometryCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; }; +``` +The caches store values that were already calculated for later usage. This increases the memory demand but makes the simulation faster. +```cpp +``` +We leave the namespace Dumux::Properties. +```cpp +} // end namespace Dumux::Properties +``` + + ## The file `spatialparams_tracer.hh` @@ -626,220 +581,26 @@ private: ## The file `problem_tracer.hh` -Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties. -### Include files -Again, we use YaspGrid, the implementation of the dune grid interface for structured grids: - -```cpp -#include <dune/grid/yaspgrid.hh> -``` - -and the cell centered, two-point-flux discretization. - -```cpp -#include <dumux/discretization/cctpfa.hh> -``` - -Here, we include the tracer model: - +This file contains the __problem class__ which defines the initial and boundary +conditions for the tracer transport simulation. ```cpp -#include <dumux/porousmediumflow/tracer/model.hh> ``` - -We include again the porous medium problem class that this class is derived from: - +### Include files +This header contains the porous medium problem class that this class is derived from: ```cpp #include <dumux/porousmediumflow/problem.hh> ``` - -and the base fluidsystem. We will define a custom fluid system that inherits from that class. - -```cpp -#include <dumux/material/fluidsystems/base.hh> -``` - -We include the header that specifies all spatially variable parameters for the tracer problem: - +This header contains the class that specifies all spatially variable parameters +related to this problem. ```cpp #include "spatialparams_tracer.hh" ``` - -### Define basic properties for our simulation -We enter the namespace Dumux - -```cpp -namespace Dumux { -``` - -The problem class is forward declared: - -```cpp -template <class TypeTag> -class TracerTestProblem; -``` - -We enter the namespace Properties, - -```cpp -namespace Properties { -``` - -A `TypeTag` for our simulation is created which inherits from the tracer model and the -cell centered, two-point-flux discretization scheme. - -```cpp -namespace TTag { -struct TracerTest { using InheritsFrom = std::tuple<Tracer>; }; -struct TracerTestCC { using InheritsFrom = std::tuple<TracerTest, CCTpfaModel>; }; -} -``` - -We enable caching for the grid volume variables, the flux variables and the FV grid geometry. - -```cpp -template<class TypeTag> -struct EnableGridVolumeVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; }; -template<class TypeTag> -struct EnableGridFluxVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; }; -template<class TypeTag> -struct EnableGridGeometryCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; }; -``` - -We use the same grid as in the stationary one-phase model, a structured 2D grid: - -```cpp -template<class TypeTag> -struct Grid<TypeTag, TTag::TracerTest> { using type = Dune::YaspGrid<2>; }; -``` - -The problem class specifies initial and boundary conditions: - -```cpp -template<class TypeTag> -struct Problem<TypeTag, TTag::TracerTest> { using type = TracerTestProblem<TypeTag>; }; -``` - -We define the spatial parameters for our tracer simulation: - -```cpp -template<class TypeTag> -struct SpatialParams<TypeTag, TTag::TracerTest> -{ -private: - using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; - using Scalar = GetPropType<TypeTag, Properties::Scalar>; -public: - using type = TracerTestSpatialParams<GridGeometry, Scalar>; -}; -``` - -One can choose between a formulation in terms of mass or mole fractions. Here, we are using mass fractions. - -```cpp -template<class TypeTag> -struct UseMoles<TypeTag, TTag::TracerTest> { static constexpr bool value = false; }; -``` - -We use solution-independent molecular diffusion coefficients. Per default, solution-dependent -diffusion coefficients are assumed during the computation of the jacobian matrix entries. Specifying -solution-independent diffusion coefficients can speed up computations: - -```cpp -template<class TypeTag> -struct SolutionDependentMolecularDiffusion<TypeTag, TTag::TracerTestCC> -{ static constexpr bool value = false; }; -``` - -In the following, we create a new tracer fluid system and derive from the base fluid system. - -```cpp -template<class TypeTag> -class TracerFluidSystem : public FluidSystems::Base<GetPropType<TypeTag, Properties::Scalar>, - TracerFluidSystem<TypeTag>> -{ -``` - -We define some convenience aliases: - -```cpp - using Scalar = GetPropType<TypeTag, Properties::Scalar>; - using Problem = GetPropType<TypeTag, Properties::Problem>; - using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; - using Element = typename GridView::template Codim<0>::Entity; - using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; - using SubControlVolume = typename FVElementGeometry::SubControlVolume; - -public: -``` - -We specify that the fluid system only contains tracer components, - -```cpp - static constexpr bool isTracerFluidSystem() - { return true; } -``` - -We define the number of components of this fluid system (one single tracer component) - -```cpp - static constexpr int numComponents = 1; -``` - -This interface is designed to define the names of the components of the fluid system. -Here, we only have a single component, so `compIdx` should always be 0. -The component name is used for the vtk output. - -```cpp - static std::string componentName(int compIdx = 0) - { return "tracer_" + std::to_string(compIdx); } -``` - -We set the phase name for the phase index (`phaseIdx`) for velocity vtk output: -Here, we only have a single phase, so `phaseIdx` should always be zero. - -```cpp - static std::string phaseName(int phaseIdx = 0) - { return "Groundwater"; } -``` - -We set the molar mass of the tracer component with index `compIdx` (should again always be zero here). - -```cpp - static Scalar molarMass(unsigned int compIdx = 0) - { return 0.300; } -``` - -We set the value for the binary diffusion coefficient. This -might depend on spatial parameters like pressure / temperature. -But, in this case we neglect diffusion and return 0.0: - -```cpp - static Scalar binaryDiffusionCoefficient(unsigned int compIdx, - const Problem& problem, - const Element& element, - const SubControlVolume& scv) - { return 0.0; } -}; -``` - -We set the above created tracer fluid system: - -```cpp -template<class TypeTag> -struct FluidSystem<TypeTag, TTag::TracerTest> { using type = TracerFluidSystem<TypeTag>; }; -``` - -We leave the namespace Properties. - -```cpp -} -``` - ### The problem class We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation. -As this is a porous medium problem, we inherit from the base class `PorousMediumFlowProblem`. - +As this is a porous medium flow problem, we inherit from the base class `PorousMediumFlowProblem`. ```cpp +namespace Dumux { + template <class TypeTag> class TracerTestProblem : public PorousMediumFlowProblem<TypeTag> { @@ -994,6 +755,174 @@ This is everything the tracer problem class contains. We leave the namespace Dumux here. ```cpp +} + +``` + + + +## The file `properties_tracer.hh` + + +This file defines the `TypeTag` used for the tracer transport simulation, for +which we then define the necessary properties. + +### Include files +As for the single-phase problem, a`TypeTag` is defined for this simulation. +Here, we inherit all properties from the `Tracer` type tag, a convenience type tag +that predefines most of the required properties for tracer transport flow simulations in DuMuX. +```cpp +#include <dumux/porousmediumflow/tracer/model.hh> +``` +Again, we use YaspGrid, an implementation of the dune grid interface for structured grids: +```cpp +#include <dune/grid/yaspgrid.hh> +``` +and the cell centered, two-point-flux discretization. +```cpp +#include <dumux/discretization/cctpfa.hh> +``` +This includes the base class for fluid systems. We will define a custom fluid +system that inherits from that class. +```cpp +#include <dumux/material/fluidsystems/base.hh> +``` +We include the problem and spatial parameters headers used for this simulation. +```cpp +#include "problem_tracer.hh" +#include "spatialparams_tracer.hh" +``` +### Basic property definitions for the tracer transport problem +We enter the namespace Dumux +```cpp +namespace Dumux { +``` +In the following, we create a new tracer fluid system and derive from the base fluid system. +```cpp +template<class TypeTag> +class TracerFluidSystem : public FluidSystems::Base<GetPropType<TypeTag, Properties::Scalar>, + TracerFluidSystem<TypeTag>> +{ +``` +We define some convenience aliases to be used inside this class. +```cpp + using Scalar = GetPropType<TypeTag, Properties::Scalar>; + using Problem = GetPropType<TypeTag, Properties::Problem>; + using GridView = GetPropType<TypeTag, Properties::GridView>; + using Element = typename GridView::template Codim<0>::Entity; + using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + +public: +``` +We specify that the fluid system only contains tracer components, +```cpp + static constexpr bool isTracerFluidSystem() + { return true; } +``` +and that no component is the main component +```cpp + static constexpr int getMainComponent(int phaseIdx) + { return -1; } +``` +We define the number of components of this fluid system (one single tracer component) +```cpp + static constexpr int numComponents = 1; +``` +This interface is designed to define the names of the components of the fluid system. +Here, we only have a single component, so `compIdx` should always be 0. +The component name is used for the vtk output. +```cpp + static std::string componentName(int compIdx = 0) + { return "tracer_" + std::to_string(compIdx); } +``` +We set the phase name for the phase index (`phaseIdx`) for velocity vtk output: +Here, we only have a single phase, so `phaseIdx` should always be zero. +```cpp + static std::string phaseName(int phaseIdx = 0) + { return "Groundwater"; } +``` +We set the molar mass of the tracer component with index `compIdx` (should again always be zero here). +```cpp + static Scalar molarMass(unsigned int compIdx = 0) + { return 0.300; } +``` +We set the value for the binary diffusion coefficient. This +might depend on spatial parameters like pressure / temperature. +But, in this case we neglect diffusion and return 0.0: +```cpp + static Scalar binaryDiffusionCoefficient(unsigned int compIdx, + const Problem& problem, + const Element& element, + const SubControlVolume& scv) + { return 0.0; } +}; +``` +We enter the namespace Properties +```cpp +namespace Properties { +``` +A `TypeTag` for our simulation is created which inherits from the tracer model and the +cell centered discretization scheme using two-point flux approximation. +```cpp +namespace TTag { +struct TracerTest { using InheritsFrom = std::tuple<Tracer>; }; +struct TracerTestCC { using InheritsFrom = std::tuple<TracerTest, CCTpfaModel>; }; +} +``` +We enable caching for the grid volume variables, the flux variables and the FV grid geometry. +```cpp +template<class TypeTag> +struct EnableGridVolumeVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; }; +template<class TypeTag> +struct EnableGridFluxVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; }; +template<class TypeTag> +struct EnableGridGeometryCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; }; +``` +We use the same grid as in the stationary one-phase model, a structured 2D grid: +```cpp +template<class TypeTag> +struct Grid<TypeTag, TTag::TracerTest> { using type = Dune::YaspGrid<2>; }; +``` +The problem class that specifies initial and boundary conditions: +```cpp +template<class TypeTag> +struct Problem<TypeTag, TTag::TracerTest> { using type = TracerTestProblem<TypeTag>; }; +``` +We define the spatial parameters for our tracer simulation: +```cpp +template<class TypeTag> +struct SpatialParams<TypeTag, TTag::TracerTest> +{ +private: + using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; + using Scalar = GetPropType<TypeTag, Properties::Scalar>; +public: + using type = TracerTestSpatialParams<GridGeometry, Scalar>; +}; +``` +One can choose between a formulation in terms of mass or mole fractions. +Here, we are using mass fractions. +```cpp +template<class TypeTag> +struct UseMoles<TypeTag, TTag::TracerTest> { static constexpr bool value = false; }; +``` +We use solution-independent molecular diffusion coefficients. Per default, solution-dependent +diffusion coefficients are assumed during the computation of the jacobian matrix entries. Specifying +solution-independent diffusion coefficients can speed up computations: +```cpp +template<class TypeTag> +struct SolutionDependentMolecularDiffusion<TypeTag, TTag::TracerTestCC> +{ static constexpr bool value = false; }; +``` +We set the above created tracer fluid system: +```cpp +template<class TypeTag> +struct FluidSystem<TypeTag, TTag::TracerTest> { using type = TracerFluidSystem<TypeTag>; }; +``` +We leave the namespace Properties and Dumux. +```cpp +} // end namespace Properties } // end namespace Dumux ``` @@ -1012,12 +941,11 @@ the transport of an initial contamination through the model domain. ```cpp #include <config.h> ``` - -We include both problems in the main file, the `problem_1p.hh` and the `problem_tracer.hh`. - +This includes the `TypeTags` and properties to be used for the single-phase +and the tracer simulations. ```cpp -#include "problem_1p.hh" -#include "problem_tracer.hh" +#include "properties_1p.hh" +#include "properties_tracer.hh" ``` Further, we include a standard header file for C++, to get time and date information diff --git a/examples/1ptracer/doc/intro.md b/examples/1ptracer/doc/intro.md index faccdde14b..9ac7130fdc 100644 --- a/examples/1ptracer/doc/intro.md +++ b/examples/1ptracer/doc/intro.md @@ -21,22 +21,22 @@ with the darcy velocity $` \textbf v `$, the permeability $` \textbf K`$, the dy Darcy's law is inserted into the mass balance equation: ```math -\phi \frac{\partial \varrho}{\partial t} + \text{div} \textbf v = 0 +\phi \frac{\partial \varrho}{\partial t} + \text{div} \textbf v = 0, ``` where $`\phi`$ is the porosity. -The equation is discretized using a cell-centered finite volume scheme as spatial discretization for the pressure as primary variable. For details on the discretization scheme, have a look at the dumux [handbook](https://dumux.org/handbook). +The equation is discretized using cell-centered finite volumes with two-point flux approximation as spatial discretization scheme for the pressure as primary variable. For details on the discretization schemes available in DuMuX, have a look at the [handbook](https://dumux.org/handbook). ### Tracer Model -The transport of the contaminant component $`\kappa`$ is based on the previously evaluated velocity field $`\textbf v`$ with the help of the following mass balance equation: +The transport of the contaminant component $`\kappa`$ occurs with the velocity field $`\textbf v`$ +that is computed with the __1p model__ (see above): ```math \phi \frac{ \partial \varrho X^\kappa}{\partial t} - \text{div} \left\lbrace \varrho X^\kappa {\textbf v} + \varrho D^\kappa_\text{pm} \textbf{grad} X^\kappa \right\rbrace = 0, ``` where $`X^\kappa`$ is the mass fraction of the contaminant component $`\kappa`$ and $` D^\kappa_\text{pm} `$ is the effective diffusivity. - The effective diffusivity is a function of the diffusion coefficient of the component $`D^\kappa`$ and the porosity and tortuosity $`\tau`$ of the porous medium (see [dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh)): ```math diff --git a/examples/1ptracer/main.cc b/examples/1ptracer/main.cc index 55d2c5ae65..240d0736fe 100644 --- a/examples/1ptracer/main.cc +++ b/examples/1ptracer/main.cc @@ -25,91 +25,89 @@ // is computed from that pressure distribution, which is then passed to a tracer problem to solve // the transport of an initial contamination through the model domain. // ### Included header files +// Some generic includes. #include <config.h> - -// We include both problems in the main file, the `problem_1p.hh` and the `problem_tracer.hh`. -#include "problem_1p.hh" -#include "problem_tracer.hh" - -// Further, we include a standard header file for C++, to get time and date information #include <ctime> -// and another one for in- and output. #include <iostream> -// Dumux is based on DUNE, the Distributed and Unified Numerics Environment, which provides several grid managers and linear solvers. -// Here, we include classes related to parallel computations, time measurements and file I/O. +// These are DUNE helper classes related to parallel computations, time measurements and file I/O #include <dune/common/parallel/mpihelper.hh> #include <dune/common/timer.hh> #include <dune/grid/io/file/dgfparser/dgfexception.hh> #include <dune/grid/io/file/vtk.hh> -// In Dumux, the property system is used to specify classes and compile-time options to be used by the model. -// For this, different properties are defined containing type definitions, values and methods. -// All properties are declared in the file `properties.hh`. +// The following headers include functionality related to property definition or retrieval, as well as +// the retrieval of input parameters specified in the input file or via the command line. #include <dumux/common/properties.hh> -// The following file contains the parameter class, which manages the definition and retrieval of input -// parameters by a default value, the inputfile or the command line. #include <dumux/common/parameters.hh> -// The file `dumuxmessage.hh` contains the class defining the start and end message of the simulation. -#include <dumux/common/dumuxmessage.hh> -// The following file contains the class, which defines the sequential linear solver backends. +// The following files contains the available linear solver backends and the assembler for the linear +// systems arising from finite volume discretizations (box-scheme, tpfa-approximation, mpfa-approximation). #include <dumux/linear/seqsolverbackend.hh> -// Further we include the assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation). #include <dumux/assembly/fvassembler.hh> -// The containing class in the following file defines the different differentiation methods used to compute the derivatives of the residual. -#include <dumux/assembly/diffmethod.hh> +#include <dumux/assembly/diffmethod.hh> // analytic or numeric differentiation -// We need the following class to simplify the writing of dumux simulation data to VTK format. +// The following class provides a convenient way of writing of dumux simulation results to VTK format. #include <dumux/io/vtkoutputmodule.hh> -// The gridmanager constructs a grid from the information in the input or grid file. There is a specification for the different supported grid managers. +// The gridmanager constructs a grid from the information in the input or grid file. +// Many different Dune grid implementations are supported, of which a list can be found +// in `gridmanager.hh`. #include <dumux/io/grid/gridmanager.hh> -// ### Beginning of the main function +// For both the single-phase and the tracer problem, `TypeTags` are defined, which collect +// the properties that are required for the simulation. These type tags are defined +// in the headers that we include here. For detailed information, please have a look +// at the documentation provided therein. +#include "properties_1p.hh" +#include "properties_tracer.hh" + +// ### The main function +// We will now discuss the main program flow implemented within the `main` function. +// At the beginning of each program using Dune, an instance `Dune::MPIHelper` has to +// be created. Moreover, we parse the run-time arguments from the command line and the +// input file: int main(int argc, char** argv) try { using namespace Dumux; - // Convenience aliases for the type tags of the two problems, which are defined in the individual problem files. - using OnePTypeTag = Properties::TTag::IncompressibleTest; - using TracerTypeTag = Properties::TTag::TracerTestCC; - - // We initialize MPI. Finalization is done automatically on exit. + // The Dune MPIHelper must be instantiated for each program using Dune const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); - // We print the dumux start message. - if (mpiHelper.rank() == 0) - DumuxMessage::print(/*firstCall=*/true); - - // We parse the command line arguments. + // parse command line arguments and input file Parameters::init(argc, argv); - // ### Create the grid + // We define convenience aliases for the type tags of the two problems. The type + // tags contain all the properties that are needed to run the simulations. Throughout + // the main file, we will obtain types defined for these type tags using the property + // system, i.e. with `GetPropType`. + using OnePTypeTag = Properties::TTag::IncompressibleTest; + using TracerTypeTag = Properties::TTag::TracerTestCC; + + // ### Step 1: Create the grid // The `GridManager` class creates the grid from information given in the input file. - // This can either be a grid file, or in the case of structured grids, by specifying the coordinates + // This can either be a grid file, or in the case of structured grids, one can specify the coordinates // of the corners of the grid and the number of cells to be used to discretize each spatial direction. - // Here, we solve both the single-phase and the tracer problem on the same grid. - // Hence, the grid is only created once using the grid type defined by the type tag of the 1p problem. + // Here, we solve both the single-phase and the tracer problem on the same grid, and thus, + // the grid is only created once using the grid type defined by the type tag of the 1p problem. GridManager<GetPropType<OnePTypeTag, Properties::Grid>> gridManager; gridManager.init(); // We compute on the leaf grid view. const auto& leafGridView = gridManager.grid().leafGridView(); - // ### Set-up and solving of the 1p problem - // In the following section, we set up and solve the 1p problem. As the result of this problem, we obtain the pressure distribution in the domain. - // #### Set-up - // We create and initialize the finite volume grid geometry, the problem, the linear system, including the jacobian matrix, the residual and the solution vector and the gridvariables. - // We need the finite volume geometry to build up the subcontrolvolumes (scv) and subcontrolvolume faces (scvf) for each element of the grid partition. + // ### Step 2: Set-up and solving of the 1p problem + // First, a finite volume grid geometry is constructed from the grid that was created above. + // This builds the subcontrolvolumes (scv) and subcontrolvolume faces (scvf) for each element + // of the grid partition. using GridGeometry = GetPropType<OnePTypeTag, Properties::GridGeometry>; auto gridGeometry = std::make_shared<GridGeometry>(leafGridView); gridGeometry->update(); - // In the problem, we define the boundary and initial conditions. + // We now instantiate the problem, in which we define the boundary and initial conditions. using OnePProblem = GetPropType<OnePTypeTag, Properties::Problem>; auto problemOneP = std::make_shared<OnePProblem>(gridGeometry); - // The jacobian matrix (`A`), the solution vector (`p`) and the residual (`r`) are parts of the linear system. + // The jacobian matrix (`A`), the solution vector (`p`) and the residual (`r`) make up the linear system. using JacobianMatrix = GetPropType<OnePTypeTag, Properties::JacobianMatrix>; using SolutionVector = GetPropType<OnePTypeTag, Properties::SolutionVector>; SolutionVector p(leafGridView.size(0)); @@ -117,61 +115,59 @@ int main(int argc, char** argv) try auto A = std::make_shared<JacobianMatrix>(); auto r = std::make_shared<SolutionVector>(); - // The grid variables store variables (primary and secondary variables) on sub-control volumes and faces (volume and flux variables). + // The grid variables are used store variables (primary and secondary variables) on sub-control volumes and faces (volume and flux variables). using OnePGridVariables = GetPropType<OnePTypeTag, Properties::GridVariables>; auto onePGridVariables = std::make_shared<OnePGridVariables>(problemOneP, gridGeometry); onePGridVariables->init(p); - // #### Assembling the linear system - // We create and inizialize the assembler. + // We now instantiate the assembler class, assemble the linear system and solve it with the linear + // solver UMFPack. Besides that, the time needed for assembly and solve is measured and printed. using OnePAssembler = FVAssembler<OnePTypeTag, DiffMethod::analytic>; auto assemblerOneP = std::make_shared<OnePAssembler>(problemOneP, gridGeometry, onePGridVariables); - assemblerOneP->setLinearSystem(A, r); + assemblerOneP->setLinearSystem(A, r); // tell assembler to use our previously defined system - // We assemble the local jacobian and the residual and stop the time needed, which is displayed in the terminal output, using the `assemblyTimer`. Further, we start the timer to evaluate the total time of the assembly, solving and updating. Dune::Timer timer; Dune::Timer assemblyTimer; std::cout << "Assembling linear system ..." << std::flush; - assemblerOneP->assembleJacobianAndResidual(p); + assemblerOneP->assembleJacobianAndResidual(p); // assemble linear system around current solution assemblyTimer.stop(); std::cout << " took " << assemblyTimer.elapsed() << " seconds." << std::endl; - // We want to solve `Ax = -r`. - (*r) *= -1.0; + (*r) *= -1.0; // We want to solve `Ax = -r`. - // #### Solution - // We set the linear solver "UMFPack" as the linear solver. Afterwards we solve the linear system. The time needed to solve the system is recorded by the `solverTimer` and displayed in the terminal output. using LinearSolver = UMFPackBackend; Dune::Timer solverTimer; std::cout << "Solving linear system ..." << std::flush; auto linearSolver = std::make_shared<LinearSolver>(); linearSolver->solve(*A, p, *r); solverTimer.stop(); std::cout << " took " << solverTimer.elapsed() << " seconds." << std::endl; - // #### Update and output - // We update the grid variables with the new solution. Dune::Timer updateTimer; std::cout << "Updating variables ..." << std::flush; - onePGridVariables->update(p); + onePGridVariables->update(p); // update grid variables to new pressure distribution updateTimer.elapsed(); std::cout << " took " << updateTimer.elapsed() << std::endl; - - // We initialize the vtkoutput. Each model has a predefined model specific output with relevant parameters for that model. We add the pressure data from the solution vector (`p`) and the permeability field as output data. - using GridView = typename GetPropType<OnePTypeTag, Properties::GridGeometry>::GridView; + // The solution vector `p` now contains the pressure field that is the solution to the single-phase + // problem defined in `problem_1p.hh`. Let us now write this solution to a VTK file using the Dune + // `VTKWriter`. Moreover, we add the permeability distribution to the writer. + using GridView = GetPropType<OnePTypeTag, Properties::GridView>; Dune::VTKWriter<GridView> onepWriter(leafGridView); onepWriter.addCellData(p, "p"); - const auto& k = problemOneP->spatialParams().getKField(); - onepWriter.addCellData(k, "permeability"); - onepWriter.write("1p"); - // We stop the timer and display the total time of the simulation as well as the cumulative CPU time. - timer.stop(); + const auto& k = problemOneP->spatialParams().getKField(); // defined in spatialparams_1p.hh + onepWriter.addCellData(k, "permeability"); // add permeability to writer + onepWriter.write("1p"); // write the file "1p.vtk" + // print overall CPU time required for assembling and solving the 1p problem. + timer.stop(); const auto& comm = Dune::MPIHelper::getCollectiveCommunication(); std::cout << "Simulation took " << timer.elapsed() << " seconds on " << comm.size() << " processes.\n" << "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n"; - // ### Computation of the volume fluxes - // We use the results of the 1p problem to calculate the volume fluxes in the model domain. - using Scalar = GetPropType<OnePTypeTag, Properties::Scalar>; + // ### Step 3: Computation of the volume fluxes + // We use the results of the 1p problem to calculate the volume fluxes across all sub-control volume + // faces of the discretization and store them in the vector `volumeFlux`. In order to do so, we iterate + // over all elements of the grid, and in each element compute the volume fluxes for all sub-control volume + // faces embeded in that element. + using Scalar = GetPropType<OnePTypeTag, Properties::Scalar>; // type for scalar values std::vector<Scalar> volumeFlux(gridGeometry->numScvf(), 0.0); using FluxVariables = GetPropType<OnePTypeTag, Properties::FluxVariables>; @@ -181,14 +177,16 @@ int main(int argc, char** argv) try for (const auto& element : elements(leafGridView)) { // Compute the element-local views on geometry, primary and secondary variables - // as well as variables needed for flux computations - auto fvGeometry = localView(*gridGeometry); - fvGeometry.bind(element); + // as well as variables needed for flux computations. + // This creates instances of the local views + auto fvGeometry = localView(*gridGeometry); auto elemVolVars = localView(onePGridVariables->curGridVolVars()); - elemVolVars.bind(element, fvGeometry, p); - auto elemFluxVars = localView(onePGridVariables->gridFluxVarsCache()); + + // we now have to bind the views to the current element + fvGeometry.bind(element); + elemVolVars.bind(element, fvGeometry, p); elemFluxVars.bind(element, fvGeometry, elemVolVars); // We calculate the volume fluxes for all sub-control volume faces except for Neumann boundary faces @@ -198,34 +196,35 @@ int main(int argc, char** argv) try if (scvf.boundary() && problemOneP->boundaryTypes(element, scvf).hasNeumann()) continue; - // let the `FluxVariables` class do the flux computation. + // let the FluxVariables class do the flux computation. FluxVariables fluxVars; fluxVars.init(*problemOneP, element, fvGeometry, elemVolVars, scvf, elemFluxVars); volumeFlux[scvf.index()] = fluxVars.advectiveFlux(0, upwindTerm); } } - - // ### Set-up and solving of the tracer problem - // #### Set-up - // Similar to the 1p problem, we first create and initialize the problem. + // ### Step 4: Set-up and solving of the tracer problem + // First, we instantiate the tracer problem containing initial and boundary conditions, + // and pass to it the previously computed volume fluxes (see the documentation of the + // file `spatialparams_tracer.hh` for more details). using TracerProblem = GetPropType<TracerTypeTag, Properties::Problem>; auto tracerProblem = std::make_shared<TracerProblem>(gridGeometry); - - // We use the volume fluxes calculated in the previous section as input for the tracer model. tracerProblem->spatialParams().setVolumeFlux(volumeFlux); - // We create and initialize the solution vector. As the tracer problem is transient, the initial solution defined in the problem is applied to the solution vector. + // We create and initialize the solution vector. As, in contrast to the steady-state single-phase problem, + // the tracer problem is transient, the initial solution defined in the problem is applied to the solution vector. + // On the basis of this solution, we initialize then the grid variables. SolutionVector x(leafGridView.size(0)); tracerProblem->applyInitialSolution(x); auto xOld = x; - // We create and initialize the grid variables. using GridVariables = GetPropType<TracerTypeTag, Properties::GridVariables>; auto gridVariables = std::make_shared<GridVariables>(tracerProblem, gridGeometry); gridVariables->init(x); - // We read in some time loop parameters from the input file. The parameter `tEnd` defines the duration of the simulation, dt the initial time step size and `maxDt` the maximal time step size. + // Let us now instantiate the time loop. Therefore, we read in some time loop parameters from the input file. + // The parameter `tEnd` defines the duration of the simulation, `dt` the initial time step size and `maxDt` the maximal time step size. + // Moreover, we define 10 check points in the time loop at which we will write the solution to vtk files. const auto tEnd = getParam<Scalar>("TimeLoop.TEnd"); auto dt = getParam<Scalar>("TimeLoop.DtInitial"); const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize"); @@ -233,26 +232,35 @@ int main(int argc, char** argv) try // We instantiate the time loop. auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd); timeLoop->setMaxTimeStepSize(maxDt); + timeLoop->setPeriodicCheckPoint(tEnd/10.0); - // We create and inizialize the assembler with time loop for the instationary problem. + // We create and initialize the assembler with a time loop for the transient problem. + // Within the time loop, we will use this assembler in each time step to assemble the linear system. using TracerAssembler = FVAssembler<TracerTypeTag, DiffMethod::analytic, /*implicit=*/false>; auto assembler = std::make_shared<TracerAssembler>(tracerProblem, gridGeometry, gridVariables, timeLoop); assembler->setLinearSystem(A, r); - // We initialize the vtk output module and add a velocity output. + // The following lines of code initialize the vtk output module, add the velocity output facility + // and write out the initial solution. At each checkpoint, we will use the output module to write + // the solution of a time step into a corresponding vtk file. VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, tracerProblem->name()); + + // add model-specific output fields to the writer using IOFields = GetPropType<TracerTypeTag, Properties::IOFields>; IOFields::initOutputModule(vtkWriter); + + // add velocity output facility using VelocityOutput = GetPropType<TracerTypeTag, Properties::VelocityOutput>; vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables)); - vtkWriter.write(0.0); - - // We define 10 check points in the time loop at which we will write the solution to vtk files. - timeLoop->setPeriodicCheckPoint(tEnd/10.0); + // write initial solution + vtkWriter.write(0.0); // #### The time loop - // We start the time loop and calculate a new time step as long as `tEnd` is not reached. In every single time step, the problem is assembled and solved. + // We start the time loop and solve a new time step as long as `tEnd` is not reached. In every time step, + // the problem is assembled and solved, the solution is updated, and when a checkpoint is reached the solution + // is written to a new vtk file. In addition, statistics related to CPU time, the current simulation time + // and the time step sizes used is printed to the terminal. timeLoop->start(); do { // First we define the old solution as the solution of the previous time step for storage evaluations. @@ -269,13 +277,13 @@ int main(int argc, char** argv) try linearSolver->solve(*A, xDelta, *r); solveTimer.stop(); - // We calculate the actual solution and update it in the grid variables. + // update the solution vector and the grid variables. updateTimer.reset(); x -= xDelta; gridVariables->update(x); updateTimer.stop(); - // We display the statistics of the actual time step. + // display the statistics of the actual time step. const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed(); std::cout << "Assemble/solve/update time: " << assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/" @@ -283,43 +291,39 @@ int main(int argc, char** argv) try << updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)" << std::endl; - // The new solution is defined as the old solution. + // Update the old solution with the one computed in this time step and move to the next one xOld = x; gridVariables->advanceTimeStep(); - - // We advance the time loop to the next time step. timeLoop->advanceTimeStep(); - // We write the Vtk output on check points. + // Write the Vtk output on check points. if (timeLoop->isCheckPoint()) vtkWriter.write(timeLoop->time()); - // We report the statistics of this time step. + // report statistics of this time step timeLoop->reportTimeStep(); - // We set the time step size dt of the next time step. + // set the time step size for the next time step timeLoop->setTimeStepSize(dt); } while (!timeLoop->finished()); + // The following piece of code prints a final status report of the time loop + // before the program is terminated. timeLoop->finalize(leafGridView.comm()); - - // ### Final Output - if (mpiHelper.rank() == 0) - DumuxMessage::print(/*firstCall=*/false); - return 0; - } // ### Exception handling // In this part of the main file we catch and print possible exceptions that could // occur during the simulation. +// errors related to run-time parameters catch (Dumux::ParameterException &e) { std::cerr << std::endl << e << " ---> Abort!" << std::endl; return 1; } +// errors related to the parsing of Dune grid files catch (Dune::DGFException & e) { std::cerr << "DGF exception thrown (" << e << @@ -329,11 +333,13 @@ catch (Dune::DGFException & e) << " ---> Abort!" << std::endl; return 2; } +// generic error handling with Dune::Exception catch (Dune::Exception &e) { std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl; return 3; } +// other exceptions catch (...) { std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; diff --git a/examples/1ptracer/problem_1p.hh b/examples/1ptracer/problem_1p.hh index b3776c45b3..5e400ff712 100644 --- a/examples/1ptracer/problem_1p.hh +++ b/examples/1ptracer/problem_1p.hh @@ -22,94 +22,21 @@ // ## The file `problem_1p.hh` // +// This file contains the __problem class__ which defines the initial and boundary +// conditions for the single-phase flow simulation. // -// Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties. // ### Include files -// We use `YaspGrid`, an implementation of the dune grid interface for structured grids: -#include <dune/grid/yaspgrid.hh> - -// The cell centered, two-point-flux discretization scheme is included: -#include <dumux/discretization/cctpfa.hh> -// The one-phase flow model is included: -#include <dumux/porousmediumflow/1p/model.hh> -// This is the porous medium problem class that this class is derived from: +// This header contains the porous medium problem class that this class is derived from: #include <dumux/porousmediumflow/problem.hh> - -// The fluid properties are specified in the following headers (we use liquid water as the fluid phase): -#include <dumux/material/components/simpleh2o.hh> -#include <dumux/material/fluidsystems/1pliquid.hh> - -// The local residual for incompressible flow is included: -#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh> - -// We include the header that specifies all spatially variable parameters: +// This header contains the class that specifies all spatially variable parameters +// related to this problem. #include "spatialparams_1p.hh" -// ### Define basic properties for our simulation -// We enter the namespace Dumux in order to import the entire Dumux namespace for general use -namespace Dumux { - -// The problem class is forward declared: -template<class TypeTag> -class OnePTestProblem; - -// We enter the namespace Properties, which is a sub-namespace of the namespace Dumux: -namespace Properties { -// A `TypeTag` for our simulation is created which inherits from the one-phase flow model and the cell centered, two-point-flux discretization scheme. -namespace TTag { -struct IncompressibleTest { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; }; -} - -// We use a structured 2D grid: -template<class TypeTag> -struct Grid<TypeTag, TTag::IncompressibleTest> { using type = Dune::YaspGrid<2>; }; - -// The problem class specifies initial and boundary conditions: -template<class TypeTag> -struct Problem<TypeTag, TTag::IncompressibleTest> { using type = OnePTestProblem<TypeTag>; }; - -// We define the spatial parameters for our simulation: -template<class TypeTag> -struct SpatialParams<TypeTag, TTag::IncompressibleTest> -{ - // We define convenient shortcuts to the properties `GridGeometry` and `Scalar`: - using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; - using Scalar = GetPropType<TypeTag, Properties::Scalar>; - // Finally, we set the spatial parameters: - using type = OnePTestSpatialParams<GridGeometry, Scalar>; -}; - -// We use the local residual that contains analytic derivative methods for incompressible flow: -template<class TypeTag> -struct LocalResidual<TypeTag, TTag::IncompressibleTest> { using type = OnePIncompressibleLocalResidual<TypeTag>; }; - -// In the following we define the fluid properties. -template<class TypeTag> -struct FluidSystem<TypeTag, TTag::IncompressibleTest> -{ - // We define a convenient shortcut to the property Scalar: - using Scalar = GetPropType<TypeTag, Properties::Scalar>; - // We create a fluid system that consists of one liquid water phase. We use the simple - // description of water, which means we do not use tabulated values but more general equations of state. - using type = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >; -}; - -// We enable caching for the grid volume variables. -template<class TypeTag> -struct EnableGridVolumeVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; }; -// We enable caching for the grid flux variables. -template<class TypeTag> -struct EnableGridFluxVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; }; -// We enable caching for the FV grid geometry -template<class TypeTag> -struct EnableGridGeometryCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; }; -// The cache stores values that were already calculated for later usage. This increases the memory demand but makes the simulation faster. -// We leave the namespace Properties. -} - // ### The problem class // We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation. -// As this is a porous medium problem, we inherit from the base class `PorousMediumFlowProblem`. +// As this is a porous medium flow problem, we inherit from the base class `PorousMediumFlowProblem`. +namespace Dumux { + template<class TypeTag> class OnePTestProblem : public PorousMediumFlowProblem<TypeTag> { diff --git a/examples/1ptracer/problem_tracer.hh b/examples/1ptracer/problem_tracer.hh index 19390c7c60..880bfa4a8b 100644 --- a/examples/1ptracer/problem_tracer.hh +++ b/examples/1ptracer/problem_tracer.hh @@ -22,133 +22,21 @@ // ## The file `problem_tracer.hh` // +// This file contains the __problem class__ which defines the initial and boundary +// conditions for the tracer transport simulation. // -//Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties. // ### Include files -// Again, we use YaspGrid, the implementation of the dune grid interface for structured grids: -#include <dune/grid/yaspgrid.hh> -// and the cell centered, two-point-flux discretization. -#include <dumux/discretization/cctpfa.hh> -// Here, we include the tracer model: -#include <dumux/porousmediumflow/tracer/model.hh> -// We include again the porous medium problem class that this class is derived from: +// This header contains the porous medium problem class that this class is derived from: #include <dumux/porousmediumflow/problem.hh> -// and the base fluidsystem. We will define a custom fluid system that inherits from that class. -#include <dumux/material/fluidsystems/base.hh> -// We include the header that specifies all spatially variable parameters for the tracer problem: +// This header contains the class that specifies all spatially variable parameters +// related to this problem. #include "spatialparams_tracer.hh" -// ### Define basic properties for our simulation -// We enter the namespace Dumux -namespace Dumux { - -// The problem class is forward declared: -template <class TypeTag> -class TracerTestProblem; - -// We enter the namespace Properties, -namespace Properties { -// A `TypeTag` for our simulation is created which inherits from the tracer model and the -// cell centered, two-point-flux discretization scheme. -namespace TTag { -struct TracerTest { using InheritsFrom = std::tuple<Tracer>; }; -struct TracerTestCC { using InheritsFrom = std::tuple<TracerTest, CCTpfaModel>; }; -} - -// We enable caching for the grid volume variables, the flux variables and the FV grid geometry. -template<class TypeTag> -struct EnableGridVolumeVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; }; -template<class TypeTag> -struct EnableGridFluxVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; }; -template<class TypeTag> -struct EnableGridGeometryCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; }; - -// We use the same grid as in the stationary one-phase model, a structured 2D grid: -template<class TypeTag> -struct Grid<TypeTag, TTag::TracerTest> { using type = Dune::YaspGrid<2>; }; - -// The problem class specifies initial and boundary conditions: -template<class TypeTag> -struct Problem<TypeTag, TTag::TracerTest> { using type = TracerTestProblem<TypeTag>; }; - -// We define the spatial parameters for our tracer simulation: -template<class TypeTag> -struct SpatialParams<TypeTag, TTag::TracerTest> -{ -private: - using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; - using Scalar = GetPropType<TypeTag, Properties::Scalar>; -public: - using type = TracerTestSpatialParams<GridGeometry, Scalar>; -}; - -// One can choose between a formulation in terms of mass or mole fractions. Here, we are using mass fractions. -template<class TypeTag> -struct UseMoles<TypeTag, TTag::TracerTest> { static constexpr bool value = false; }; -// We use solution-independent molecular diffusion coefficients. Per default, solution-dependent -// diffusion coefficients are assumed during the computation of the jacobian matrix entries. Specifying -// solution-independent diffusion coefficients can speed up computations: -template<class TypeTag> -struct SolutionDependentMolecularDiffusion<TypeTag, TTag::TracerTestCC> -{ static constexpr bool value = false; }; - -// In the following, we create a new tracer fluid system and derive from the base fluid system. -template<class TypeTag> -class TracerFluidSystem : public FluidSystems::Base<GetPropType<TypeTag, Properties::Scalar>, - TracerFluidSystem<TypeTag>> -{ - // We define some convenience aliases: - using Scalar = GetPropType<TypeTag, Properties::Scalar>; - using Problem = GetPropType<TypeTag, Properties::Problem>; - using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; - using Element = typename GridView::template Codim<0>::Entity; - using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; - using SubControlVolume = typename FVElementGeometry::SubControlVolume; - -public: - // We specify that the fluid system only contains tracer components, - static constexpr bool isTracerFluidSystem() - { return true; } - - // We define the number of components of this fluid system (one single tracer component) - static constexpr int numComponents = 1; - - // This interface is designed to define the names of the components of the fluid system. - // Here, we only have a single component, so `compIdx` should always be 0. - // The component name is used for the vtk output. - static std::string componentName(int compIdx = 0) - { return "tracer_" + std::to_string(compIdx); } - - // We set the phase name for the phase index (`phaseIdx`) for velocity vtk output: - // Here, we only have a single phase, so `phaseIdx` should always be zero. - static std::string phaseName(int phaseIdx = 0) - { return "Groundwater"; } - - // We set the molar mass of the tracer component with index `compIdx` (should again always be zero here). - static Scalar molarMass(unsigned int compIdx = 0) - { return 0.300; } - - // We set the value for the binary diffusion coefficient. This - // might depend on spatial parameters like pressure / temperature. - // But, in this case we neglect diffusion and return 0.0: - static Scalar binaryDiffusionCoefficient(unsigned int compIdx, - const Problem& problem, - const Element& element, - const SubControlVolume& scv) - { return 0.0; } -}; - -// We set the above created tracer fluid system: -template<class TypeTag> -struct FluidSystem<TypeTag, TTag::TracerTest> { using type = TracerFluidSystem<TypeTag>; }; - -// We leave the namespace Properties. -} - - // ### The problem class // We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation. -// As this is a porous medium problem, we inherit from the base class `PorousMediumFlowProblem`. +// As this is a porous medium flow problem, we inherit from the base class `PorousMediumFlowProblem`. +namespace Dumux { + template <class TypeTag> class TracerTestProblem : public PorousMediumFlowProblem<TypeTag> { @@ -244,7 +132,7 @@ public: } private: -// We assign a private global variable for the epsilon: + // We assign a private global variable for the epsilon: static constexpr Scalar eps_ = 1e-6; // This is everything the tracer problem class contains. diff --git a/examples/1ptracer/properties_1p.hh b/examples/1ptracer/properties_1p.hh new file mode 100644 index 0000000000..9a782908f1 --- /dev/null +++ b/examples/1ptracer/properties_1p.hh @@ -0,0 +1,112 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ + +// ### Header guard +#ifndef DUMUX_ONEP_TEST_PROPERTIES_HH +#define DUMUX_ONEP_TEST_PROPERTIES_HH + +// This file defines the `TypeTag` used for the single-phase simulation, for +// which we then define the necessary properties. +// +// ### Include files +// The `TypeTag` defined for this simulation will inherit all properties from the +// `OneP` type tag, a convenience type tag that predefines most of the required +// properties for single-phase flow simulations in DuMuX. The properties that will be +// defined in this file are those that depend on user choices and no meaningful +// default can be set. +#include <dumux/porousmediumflow/1p/model.hh> + +// We want to use `YaspGrid`, an implementation of the dune grid interface for structured grids: +#include <dune/grid/yaspgrid.hh> +// In this example, we want to discretize the equations with the cell centered finite volume +// scheme using two-point-flux approximation: +#include <dumux/discretization/cctpfa.hh> +// The fluid properties are specified in the following headers (we use liquid water as the fluid phase): +#include <dumux/material/components/simpleh2o.hh> +#include <dumux/material/fluidsystems/1pliquid.hh> + +// The local residual for incompressible flow is included. +// The one-phase flow model (included above) uses a default implementation of the +// local residual for single-phase flow. However, in this example we are using an +// incompressible fluid phase. Therefore, we are including the specialized local +// residual which contains functionality to analytically compute the entries of +// the Jacobian matrix. We will use this in the main file. +#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh> + +// We include the problem and spatial parameters headers used for this simulation. +#include "problem_1p.hh" +#include "spatialparams_1p.hh" + +// ### Basic property definitions for the 1p problem +// We enter the namespace Dumux::Properties in order to import the entire Dumux namespace for general use: +namespace Dumux:: Properties { + +// A `TypeTag` for our simulation is created which inherits from the one-phase flow model +// and the cell centered finite volume scheme with two-point-flux discretization scheme: +namespace TTag { +struct IncompressibleTest { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; }; +} + +// We use a structured 2D grid: +template<class TypeTag> +struct Grid<TypeTag, TTag::IncompressibleTest> { using type = Dune::YaspGrid<2>; }; + +// The problem class specifies initial and boundary conditions: +template<class TypeTag> +struct Problem<TypeTag, TTag::IncompressibleTest> { using type = OnePTestProblem<TypeTag>; }; + +// We define the spatial parameters for our simulation: +template<class TypeTag> +struct SpatialParams<TypeTag, TTag::IncompressibleTest> +{ +private: + using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; + using Scalar = GetPropType<TypeTag, Properties::Scalar>; +public: + using type = OnePTestSpatialParams<GridGeometry, Scalar>; +}; + +// We use the local residual that contains analytic derivative methods for incompressible flow: +template<class TypeTag> +struct LocalResidual<TypeTag, TTag::IncompressibleTest> { using type = OnePIncompressibleLocalResidual<TypeTag>; }; + +// In the following we define the fluid system to be used: +template<class TypeTag> +struct FluidSystem<TypeTag, TTag::IncompressibleTest> +{ +private: + using Scalar = GetPropType<TypeTag, Properties::Scalar>; +public: + using type = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >; +}; + +// This enables grid-wide caching of the volume variables. +template<class TypeTag> +struct EnableGridVolumeVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; }; +//This enables grid wide caching for the flux variables. +template<class TypeTag> +struct EnableGridFluxVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; }; +// This enables grid-wide caching for the finite volume grid geometry +template<class TypeTag> +struct EnableGridGeometryCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; }; +// The caches store values that were already calculated for later usage. This increases the memory demand but makes the simulation faster. + +// We leave the namespace Dumux::Properties. +} // end namespace Dumux::Properties +#endif diff --git a/examples/1ptracer/properties_tracer.hh b/examples/1ptracer/properties_tracer.hh new file mode 100644 index 0000000000..d98cf4c6f2 --- /dev/null +++ b/examples/1ptracer/properties_tracer.hh @@ -0,0 +1,156 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ + +// ### Header guard +#ifndef DUMUX_TRACER_TEST_PROPERTIES_HH +#define DUMUX_TRACER_TEST_PROPERTIES_HH + +// This file defines the `TypeTag` used for the tracer transport simulation, for +// which we then define the necessary properties. +// +// ### Include files +// As for the single-phase problem, a`TypeTag` is defined for this simulation. +// Here, we inherit all properties from the `Tracer` type tag, a convenience type tag +// that predefines most of the required properties for tracer transport flow simulations in DuMuX. +#include <dumux/porousmediumflow/tracer/model.hh> + +// Again, we use YaspGrid, an implementation of the dune grid interface for structured grids: +#include <dune/grid/yaspgrid.hh> +// and the cell centered, two-point-flux discretization. +#include <dumux/discretization/cctpfa.hh> +// This includes the base class for fluid systems. We will define a custom fluid +// system that inherits from that class. +#include <dumux/material/fluidsystems/base.hh> + +// We include the problem and spatial parameters headers used for this simulation. +#include "problem_tracer.hh" +#include "spatialparams_tracer.hh" + +// ### Basic property definitions for the tracer transport problem +// We enter the namespace Dumux +namespace Dumux { + +// In the following, we create a new tracer fluid system and derive from the base fluid system. +template<class TypeTag> +class TracerFluidSystem : public FluidSystems::Base<GetPropType<TypeTag, Properties::Scalar>, + TracerFluidSystem<TypeTag>> +{ + // We define some convenience aliases to be used inside this class. + using Scalar = GetPropType<TypeTag, Properties::Scalar>; + using Problem = GetPropType<TypeTag, Properties::Problem>; + using GridView = GetPropType<TypeTag, Properties::GridView>; + using Element = typename GridView::template Codim<0>::Entity; + using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + +public: + // We specify that the fluid system only contains tracer components, + static constexpr bool isTracerFluidSystem() + { return true; } + + // and that no component is the main component + static constexpr int getMainComponent(int phaseIdx) + { return -1; } + + // We define the number of components of this fluid system (one single tracer component) + static constexpr int numComponents = 1; + + // This interface is designed to define the names of the components of the fluid system. + // Here, we only have a single component, so `compIdx` should always be 0. + // The component name is used for the vtk output. + static std::string componentName(int compIdx = 0) + { return "tracer_" + std::to_string(compIdx); } + + // We set the phase name for the phase index (`phaseIdx`) for velocity vtk output: + // Here, we only have a single phase, so `phaseIdx` should always be zero. + static std::string phaseName(int phaseIdx = 0) + { return "Groundwater"; } + + // We set the molar mass of the tracer component with index `compIdx` (should again always be zero here). + static Scalar molarMass(unsigned int compIdx = 0) + { return 0.300; } + + // We set the value for the binary diffusion coefficient. This + // might depend on spatial parameters like pressure / temperature. + // But, in this case we neglect diffusion and return 0.0: + static Scalar binaryDiffusionCoefficient(unsigned int compIdx, + const Problem& problem, + const Element& element, + const SubControlVolume& scv) + { return 0.0; } +}; + +// We enter the namespace Properties +namespace Properties { + +// A `TypeTag` for our simulation is created which inherits from the tracer model and the +// cell centered discretization scheme using two-point flux approximation. +namespace TTag { +struct TracerTest { using InheritsFrom = std::tuple<Tracer>; }; +struct TracerTestCC { using InheritsFrom = std::tuple<TracerTest, CCTpfaModel>; }; +} + +// We enable caching for the grid volume variables, the flux variables and the FV grid geometry. +template<class TypeTag> +struct EnableGridVolumeVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; }; +template<class TypeTag> +struct EnableGridFluxVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; }; +template<class TypeTag> +struct EnableGridGeometryCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; }; + +// We use the same grid as in the stationary one-phase model, a structured 2D grid: +template<class TypeTag> +struct Grid<TypeTag, TTag::TracerTest> { using type = Dune::YaspGrid<2>; }; + +// The problem class that specifies initial and boundary conditions: +template<class TypeTag> +struct Problem<TypeTag, TTag::TracerTest> { using type = TracerTestProblem<TypeTag>; }; + +// We define the spatial parameters for our tracer simulation: +template<class TypeTag> +struct SpatialParams<TypeTag, TTag::TracerTest> +{ +private: + using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; + using Scalar = GetPropType<TypeTag, Properties::Scalar>; +public: + using type = TracerTestSpatialParams<GridGeometry, Scalar>; +}; + +// One can choose between a formulation in terms of mass or mole fractions. +// Here, we are using mass fractions. +template<class TypeTag> +struct UseMoles<TypeTag, TTag::TracerTest> { static constexpr bool value = false; }; + +// We use solution-independent molecular diffusion coefficients. Per default, solution-dependent +// diffusion coefficients are assumed during the computation of the jacobian matrix entries. Specifying +// solution-independent diffusion coefficients can speed up computations: +template<class TypeTag> +struct SolutionDependentMolecularDiffusion<TypeTag, TTag::TracerTestCC> +{ static constexpr bool value = false; }; + +// We set the above created tracer fluid system: +template<class TypeTag> +struct FluidSystem<TypeTag, TTag::TracerTest> { using type = TracerFluidSystem<TypeTag>; }; + +// We leave the namespace Properties and Dumux. +} // end namespace Properties +} // end namespace Dumux + +#endif diff --git a/examples/1ptracer/spatialparams_1p.hh b/examples/1ptracer/spatialparams_1p.hh index 5cfb74f053..8befc0c6f1 100644 --- a/examples/1ptracer/spatialparams_1p.hh +++ b/examples/1ptracer/spatialparams_1p.hh @@ -23,26 +23,30 @@ // ## The file `spatialparams_1p.hh` // // -// In this file, we generate a random permeability field in the constructor of the `OnePTestSpatialParams` class. -// For this, we use the random number generation facilities provided by the C++ standard library. +// This file contains the __spatial parameters class__ which defines the +// distributions for the porous medium parameters permeability and porosity +// over the computational grid +// +// In this example, we use a randomly generated and element-wise distributed +// permeability field. For this, we use the random number generation facilitie +// provided by the C++ standard library. #include <random> -// We use the properties for porous medium flow models, declared in the file `properties.hh`. -#include <dumux/porousmediumflow/properties.hh> - -// We include the spatial parameters class for single-phase models discretized by finite volume schemes. -// The spatial parameters defined for this example will inherit from those. +// We include the spatial parameters class for single-phase models discretized +// by finite volume schemes, from which the spatial parameters defined for this +// example will inherit. #include <dumux/material/spatialparams/fv1p.hh> namespace Dumux { -// In the `OnePTestSpatialParams` class, we define all functions needed to describe the porous medium, e.g. porosity and permeability for the 1p_problem. +// In the `OnePTestSpatialParams` class, we define all functions needed to describe +// the porous medium, e.g. porosity and permeability, for the 1p_problem. template<class GridGeometry, class Scalar> class OnePTestSpatialParams : public FVSpatialParamsOneP<GridGeometry, Scalar, OnePTestSpatialParams<GridGeometry, Scalar>> { - // We declare aliases for types that we are going to need in this class. + // The following convenience aliases will be used throughout this class: using GridView = typename GridGeometry::GridView; using FVElementGeometry = typename GridGeometry::LocalView; using SubControlVolume = typename FVElementGeometry::SubControlVolume; @@ -58,15 +62,17 @@ public: // The spatial parameters must export the type used to define permeabilities. // Here, we are using scalar permeabilities, but tensors are also supported. using PermeabilityType = Scalar; + + // ### Generation of the random permeability field + // We generate the random permeability field upon construction of the spatial parameters class OnePTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry) : ParentType(gridGeometry), K_(gridGeometry->gridView().size(0), 0.0) { - // ### Generation of the random permeability field - // We get the permeability of the domain and the lens from the `params.input` file. + // The permeability of the domain and the lens are obtained from the `params.input` file. permeability_ = getParam<Scalar>("SpatialParams.Permeability"); permeabilityLens_ = getParam<Scalar>("SpatialParams.PermeabilityLens"); - // Furthermore, we get the position of the lens, which is defined by the position of the lower left and the upper right corner. + // Furthermore, the position of the lens, which is defined by the position of the lower left and the upper right corners, are obtained from the input file. lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft"); lensUpperRight_ =getParam<GlobalPosition>("SpatialParams.LensUpperRight"); @@ -92,12 +98,12 @@ public: const SubControlVolume& scv, const ElementSolution& elemSol) const { - return K_[scv.dofIndex()]; + return K_[scv.elementIndex()]; } // We set the porosity $`[-]`$ for the whole domain to a value of $`20 \%`$. - Scalar porosityAtPos(const GlobalPosition &globalPos) const + Scalar porosityAtPos(const GlobalPosition& globalPos) const { return 0.2; } // We reference to the permeability field. This is used in the main function to write an output for the permeability field. @@ -105,24 +111,23 @@ public: { return K_; } private: - // We have a convenient definition of the position of the lens. - bool isInLens_(const GlobalPosition &globalPos) const + // The following function returns true if a given position is inside the lens. + // We use an epsilon of 1.5e-7 here for floating point comparisons. + bool isInLens_(const GlobalPosition& globalPos) const { - for (int i = 0; i < dimWorld; ++i) { - if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_) + for (int i = 0; i < dimWorld; ++i) + { + if (globalPos[i] < lensLowerLeft_[i] + 1.5e-7 + || globalPos[i] > lensUpperRight_[i] - 1.5e-7) return false; } + return true; } - GlobalPosition lensLowerLeft_; - GlobalPosition lensUpperRight_; - + GlobalPosition lensLowerLeft_, lensUpperRight_; Scalar permeability_, permeabilityLens_; - std::vector<Scalar> K_; - - static constexpr Scalar eps_ = 1.5e-7; }; } // end namespace Dumux -- GitLab