diff --git a/examples/2pinfiltration/README.md b/examples/2pinfiltration/README.md
index 9d179cafb3a3ae13165150ec0fb8cd9c545cb6cb..20da8e72f56b837e8e9f1fa88265b23e32d5f64a 100644
--- a/examples/2pinfiltration/README.md
+++ b/examples/2pinfiltration/README.md
@@ -25,14 +25,14 @@ __Table of contents__. This description is structured as follows:
## Scenario and mathematical model
We model a soil contamination problem where DNAPL infiltrates a porous medium. The initial distribution of DNAPL is known and we can read it from a txt-file.
-To describe that problem we use a two phase model of two immiscible fluids with the multiphase Darcy's law as the description of momentum, i.e.:
+We describe the problem using a two phase model with two immiscible fluid phases (subscripts $`w`$ and $`n`$). We use multiphase Darcy's law for the momentum balance equations of the fluid phases, i.e.:
```math
v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K}
\left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right)
```
-If we insert this into the conservation equations for each phase $`\alpha`$ that leads to:
+Inserting this into the conservation equations for each phase $`\alpha`$ yields:
```math
\phi \frac{\partial \varrho_\alpha S_\alpha}{\partial t}
@@ -40,10 +40,16 @@ If we insert this into the conservation equations for each phase $`\alpha`$ that
\right\} - q_\alpha = 0
```
-To reduce the number of unknowns and close the system we need closure relations for this equations. For that, we make use of a $`p_c - S_w`$ as well as a $`k_r - S_w`$ - relationship. In this problem we use a Van-Genuchten parameterization. The parameters for that relationship are specified in the `spatialparams.hh` file.
+To reduce the number of unknowns and close the system we need closure relations for these equations. In this example, we use the
+Van Genuchten-Mualem relationships (see
+[Van Genuchten (1980)](https://doi.org/10.2136/sssaj1980.03615995004400050002x)
+and
+[Mualem (1976)](https://doi.org/10.1029/WR012i003p00513))
+for the capillary pressure $`pc = p_n - p_w`$ and the relative permeabilities $`k_r\alpha`$.
+The parameters for these relationships are specified in the `spatialparams.hh` file.
-With the additional constraint that $`S_w + S_n = 1`$ we reduce the number of primary variables to two.
-In this example we use the wetting phase pressure $`p_0`$ and the saturation of the nonwetting phase $`S_1`$ as primary variables. It is also possible to switch that formulation to the nonwetting pressure and the wetting saturation.
+With the additional constraint that $`S_w + S_n = 1`$, the number of unknowns is reduced to two.
+In this example we use the wetting phase pressure $`p_w`$ and the saturation of the nonwetting phase $`S_n`$ as primary variables. It is also possible to switch that formulation to the nonwetting pressure and the wetting saturation.
The two-dimensional model domain is 6m x 4m and contains a lens with a lower permeability and porosity. We read the initial values for the DNAPL saturation and the water pressure from a file.
The lens and the initial saturation can be seen in Figures 1 and 2.
@@ -55,7 +61,7 @@ DNAPL enters the model domain at the upper boundary between 1.75m ≤ x ≤ 2m w
In addition, the DNAPL is injected at a point source at x = 0.502m and y = 3.02m with a rate of 0.1 kg/s.
We discretize the equations with a cell-centered finite volume TPFA scheme in space and an implicit Euler scheme in time. We use Newton's method to solve the system of nonlinear equations.
-The grid is adapitvely refined around the injection. The adaptive behaviour can be changed with input parameters in the `params.input` file.
+The grid is adaptively refined around the injection. The adaptive behaviour can be changed with input parameters in the `params.input` file.
# Implementation
diff --git a/examples/2pinfiltration/doc/2p.md b/examples/2pinfiltration/doc/2p.md
index c7e5e94a6179092ad2d23ff0b75f21e8d1c58f22..5f522154db4a36ed0d315bd66b6a0e2292d192bc 100644
--- a/examples/2pinfiltration/doc/2p.md
+++ b/examples/2pinfiltration/doc/2p.md
@@ -128,7 +128,9 @@ struct FluidSystem<TypeTag, TTag::PointSourceExample>
The two-phase model implements different primary variable formulations.
Here we choose the pressure of the first phase and the saturation of the second phase.
The order of phases is specified by the fluid system.
-In this case that means that the primary variables are water pressure and DNAPL saturation.
+With our choice of FluidSystem this means that the first phase is water and the second
+phase is DNAPL. We choose the `p0s1` formulation, which means a $`p_w - S_n`$ formulation,
+as the first phase (i.e. water) is the wetting phase in this example.
```cpp
template<class TypeTag>
@@ -217,10 +219,11 @@ public:
```
#### Boundary types
-We define the type of boundary conditions depending on location. Two types of boundary conditions
-can be specified: Dirichlet or Neumann boundary condition. On a Dirichlet boundary, the values of the
+We define the type of boundary conditions depending on the location. Two types of boundary conditions
+can be specified: Dirichlet or Neumann boundary conditions. On a Dirichlet boundary, the values of the
primary variables need to be fixed. On a Neumann boundary condition, values for derivatives need to be fixed.
-Mixed boundary conditions (different types for different equations on the same boundary) are not accepted.
+Depending on the chosen discretization scheme, mixed boundary conditions (different types for different equations
+on the same boundary) may not be accepted.
```cpp
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
@@ -462,7 +465,8 @@ We set the parameter object for the Van Genuchten material law.
}
```
-Here we can define which phase is to be considered as the wetting phase. Here the wetting phase is the the first phase of the fluidsystem. In this case that is water.
+Here we can define which phase is to be considered as the wetting phase. We choose the first phase
+in the fluid system, which is water in our case (see definition of `FluidSystem` property).
```cpp
template<class FluidSystem>
diff --git a/examples/2pinfiltration/doc/_intro.md b/examples/2pinfiltration/doc/_intro.md
index 2a5af295469cdfb328100983d8b71be2a1c99177..75982c3351463060846f5ebe096eb0a17cae11d1 100644
--- a/examples/2pinfiltration/doc/_intro.md
+++ b/examples/2pinfiltration/doc/_intro.md
@@ -23,14 +23,14 @@ __Table of contents__. This description is structured as follows:
## Scenario and mathematical model
We model a soil contamination problem where DNAPL infiltrates a porous medium. The initial distribution of DNAPL is known and we can read it from a txt-file.
-To describe that problem we use a two phase model of two immiscible fluids with the multiphase Darcy's law as the description of momentum, i.e.:
+We describe the problem using a two phase model with two immiscible fluid phases (subscripts $`w`$ and $`n`$). We use multiphase Darcy's law for the momentum balance equations of the fluid phases, i.e.:
```math
v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K}
\left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right)
```
-If we insert this into the conservation equations for each phase $`\alpha`$ that leads to:
+Inserting this into the conservation equations for each phase $`\alpha`$ yields:
```math
\phi \frac{\partial \varrho_\alpha S_\alpha}{\partial t}
@@ -38,10 +38,16 @@ If we insert this into the conservation equations for each phase $`\alpha`$ that
\right\} - q_\alpha = 0
```
-To reduce the number of unknowns and close the system we need closure relations for this equations. For that, we make use of a $`p_c - S_w`$ as well as a $`k_r - S_w`$ - relationship. In this problem we use a Van-Genuchten parameterization. The parameters for that relationship are specified in the `spatialparams.hh` file.
+To reduce the number of unknowns and close the system we need closure relations for these equations. In this example, we use the
+Van Genuchten-Mualem relationships (see
+[Van Genuchten (1980)](https://doi.org/10.2136/sssaj1980.03615995004400050002x)
+and
+[Mualem (1976)](https://doi.org/10.1029/WR012i003p00513))
+for the capillary pressure $`pc = p_n - p_w`$ and the relative permeabilities $`k_r\alpha`$.
+The parameters for these relationships are specified in the `spatialparams.hh` file.
-With the additional constraint that $`S_w + S_n = 1`$ we reduce the number of primary variables to two.
-In this example we use the wetting phase pressure $`p_0`$ and the saturation of the nonwetting phase $`S_1`$ as primary variables. It is also possible to switch that formulation to the nonwetting pressure and the wetting saturation.
+With the additional constraint that $`S_w + S_n = 1`$, the number of unknowns is reduced to two.
+In this example we use the wetting phase pressure $`p_w`$ and the saturation of the nonwetting phase $`S_n`$ as primary variables. It is also possible to switch that formulation to the nonwetting pressure and the wetting saturation.
The two-dimensional model domain is 6m x 4m and contains a lens with a lower permeability and porosity. We read the initial values for the DNAPL saturation and the water pressure from a file.
The lens and the initial saturation can be seen in Figures 1 and 2.
@@ -53,7 +59,7 @@ DNAPL enters the model domain at the upper boundary between 1.75m ≤ x ≤ 2m w
In addition, the DNAPL is injected at a point source at x = 0.502m and y = 3.02m with a rate of 0.1 kg/s.
We discretize the equations with a cell-centered finite volume TPFA scheme in space and an implicit Euler scheme in time. We use Newton's method to solve the system of nonlinear equations.
-The grid is adapitvely refined around the injection. The adaptive behaviour can be changed with input parameters in the `params.input` file.
+The grid is adaptively refined around the injection. The adaptive behaviour can be changed with input parameters in the `params.input` file.
# Implementation
diff --git a/examples/2pinfiltration/doc/main.md b/examples/2pinfiltration/doc/main.md
index ce6354fc082225c3aed8c4316c727c6352e7e5d3..468bbb09c7ed72c0f3080e796fa07f222da21d03 100644
--- a/examples/2pinfiltration/doc/main.md
+++ b/examples/2pinfiltration/doc/main.md
@@ -6,8 +6,8 @@
# Part 2: Main program flow
-This main program flow is implemented in the `main()` function
-of the program which is defined in the file `main.cc` described below. Additionally, we can see how to implement an adaptive grid with a saturation dependent indicator.
+This main program flow is implemented in the `main` function
+of the program which is defined in the file `main.cc` described below. Additionally, we can see how to implement an adaptive grid with a saturation-dependent indicator.
The code documentation is structured as follows:
@@ -92,7 +92,9 @@ Finally, we include the properties which configure the simulation
</details>
### The main function
-At the beginning of the simulation, we create a type tag with a suitable alias for our problem. // This will contain all the properties that are needed to define the problem set-up and the model we use. Additionally, we have to create an instance of `Dune::MPIHelper` and parse the run-time arguments:
+At the beginning of the simulation, we create an alias for our problem type tag, for which we have
+defined the properties of our simulation (see `properties.hh`). Additionally, we have to create an
+instance of `Dune::MPIHelper` and parse the run-time arguments:
```cpp
int main(int argc, char** argv) try
@@ -102,16 +104,17 @@ int main(int argc, char** argv) try
// we define the type tag for this problem
using TypeTag = Properties::TTag::PointSourceExample;
- // maybe initialize MPI and/or multithreading backend
+ // (maybe) initialize MPI and/or multithreading backend
Dumux::initialize(argc, argv);
const auto& mpiHelper = Dune::MPIHelper::instance();
- //We parse command line arguments and input file
+ // Construct parameter tree from command line arguments (and input file if given in the arguments)
Parameters::init(argc, argv);
```
#### Create the grid
-A gridmanager tries to create the grid either from a grid file or the input file.
+The `gridManager` creates a grid from our parameter choices in the input file (which could be a grid
+file or specified domain dimensions and number of cells, as done in this example).
```cpp
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
@@ -124,26 +127,27 @@ A gridmanager tries to create the grid either from a grid file or the input file
```
#### Set-up of the problem
-
-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.
+We build the finite volume geometry, which allows us to iterate over subcontrolvolumes (scv) and
+subcontrolvolume faces (scvf) embedded in the elements of the grid partition.
```cpp
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
```
-In the problem, we define the boundary and initial conditions and compute the point sources. The `computePointSourceMap` method is inherited from the fvproblem and therefore specified in the `dumux/common/fvproblem.hh`. It calls the `addPointSources` method specified in the `problem.hh` file.
+In the problem, we define the boundary and initial conditions and compute the point sources.
+The `computePointSourceMap` function is inherited from `FVProblem` (see `dumux/common/fvproblem.hh`).
+It calls the `addPointSources` method, which is empty per default, but can be overloaded in user problems.
+As we have seen, in this example we do specify a point source (see `problem.hh`).
```cpp
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
- // We call the `computePointSourceMap` method to compute the point sources.
problem->computePointSourceMap();
```
-We initialize the solution vector and then use the solution vector to initialize the `gridVariables`.
+We initialize the solution vector (all primary variables on the grid) and then use it to initialize the
+`gridVariables`, which provide access to both primary & secondary variables (per sub-control volume).
```cpp
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
@@ -157,7 +161,8 @@ We initialize the solution vector and then use the solution vector to initialize
```
##### Grid adaption
-We instantiate the indicator for grid adaption & the data transfer, we read some parameters for indicator from the input file.
+The following piece of code shows how we instantiate an indicator class which determines where to adapt the grid,
+and a data transfer class that maps a solution to an adapted grid.
```cpp
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -173,63 +178,45 @@ We instantiate the indicator for grid adaption & the data transfer, we read some
TwoPGridDataTransfer<TypeTag> dataTransfer(problem, gridGeometry, gridVariables, x);
```
-We do an initial refinement around sources/BCs. We use the `GridAdaptInitializationIndicator` defined in `dumux/adaptive/initializationindicator.hh` for that.
-
-```cpp
- GridAdaptInitializationIndicator<TypeTag> initIndicator(problem, gridGeometry, gridVariables);
-```
-
-We refine up to the maximum level. For every level, the indicator used for the refinement/coarsening is calculated. If any grid cells have to be adapted, the gridvariables and the pointsourcemap are updated.
+Convenience functions to perform a grid adaptation given a pre-calculated indicator.
```cpp
- const auto maxLevel = getParam<std::size_t>("Adaptive.MaxLevel", 0);
- for (std::size_t i = 0; i < maxLevel; ++i)
- {
- //we calculate the initial indicator for adaption for each grid cell using the initial solution x
- initIndicator.calculate(x);
-
- //and then we mark the elements that were adapted.
+ const auto adaptGridAndVariables = [&] (auto& indicator) {
+ // We mark and adapt the elements according to the indicator.
bool wasAdapted = false;
- if (markElements(gridManager.grid(), initIndicator))
+ if (markElements(gridManager.grid(), indicator))
wasAdapted = adapt(gridManager.grid(), dataTransfer);
- // In case of a grid adaptation, the gridvariables and the pointsourcemap are updated.
+ // In case of any adapted elements, the gridvariables and the pointsourcemap are updated.
if (wasAdapted)
{
// We overwrite the old solution with the new (resized & interpolated) one
xOld = x;
- //We initialize the secondary variables to the new (and "new old") solution
+ // We initialize the secondary variables to the new (and "new old") solution
gridVariables->updateAfterGridAdaption(x);
// we update the point source map after adaption
problem->computePointSourceMap();
}
- }
+ };
```
-Depending on the initial conditions, another grid adaptation might be necessary.
-The gridadaptindicator uses the input parameters `Adaptive.RefineTolerance` and `Adaptive.CoarsenTolerance` for this step.
-Again, if elements were adapted, we mark them and update gridvariables and the pointsourcemap accordingly.
+We do initial refinements around sources/BCs, for which we use the `GridAdaptInitializationIndicator` defined in `dumux/adaptive/initializationindicator.hh`.
+Afterwards, depending on the initial conditions, another grid adaptation using our `indicator` above might be necessary, which uses
+`Adaptive.RefineTolerance` and `Adaptive.CoarsenTolerance` for this step.
```cpp
- indicator.calculate(x, refineTol, coarsenTol);
-
- //we mark the elements that were adapted
- bool wasAdapted = false;
- if (markElements(gridManager.grid(), indicator))
- wasAdapted = adapt(gridManager.grid(), dataTransfer);
-
- // In case of an additional grid adaptation, the gridvariables and the pointsourcemap are updated again.
- if (wasAdapted)
+ GridAdaptInitializationIndicator<TypeTag> initIndicator(problem, gridGeometry, gridVariables);
+ const auto maxLevel = getParam<std::size_t>("Adaptive.MaxLevel", 0);
+ for (std::size_t i = 0; i < maxLevel; ++i)
{
- // Overwrite the old solution with the new (resized & interpolated) one
- xOld = x;
- // Initialize the secondary variables to the new (and "new old") solution
- gridVariables->updateAfterGridAdaption(x);
- // Update the point source map
- problem->computePointSourceMap();
+ initIndicator.calculate(x);
+ adaptGridAndVariables(initIndicator);
}
+ indicator.calculate(x, refineTol, coarsenTol);
+ adaptGridAndVariables(indicator);
```
+[[/codeblock]]
#### Solving the problem
We get some time loop parameters from the input file params.input
@@ -283,50 +270,33 @@ We start the time loop. In each time step before we start calculating a new solu
timeLoop->start(); do
{
// We only want to refine/coarsen after first time step is finished, not before.
- //The initial refinement was already done before the start of the time loop.
- //This means we only refine when the time is greater than 0.
+ // The initial refinement was already done before the start of the time loop.
+ // This means we only refine when the time is greater than 0. Note that we now also
+ // have to update the assembler, since the sizes of the residual vector and jacobian matrix change.
if (timeLoop->time() > 0)
{
- // again we compute the refinement indicator with the `TwoPGridAdaptIndicator`
indicator.calculate(x, refineTol, coarsenTol);
-
- //we mark elements and adapt grid if necessary
- wasAdapted = false;
- if (markElements(gridManager.grid(), indicator))
- wasAdapted = adapt(gridManager.grid(), dataTransfer);
-
- //In case of a grid adaptation, the gridvariables and the pointsourcemap are updated again.
- if (wasAdapted)
- {
- // We overwrite the old solution with the new (resized & interpolated) one
- xOld = x;
- // We initialize the secondary variables to the new (and "new old") solution
- gridVariables->updateAfterGridAdaption(x);
- // We also need to update the assembler (sizes of residual and Jacobian change)
- assembler->updateAfterGridAdaption();
- // We update the point source map
- problem->computePointSourceMap();
- }
- // we leave the refinement step
+ adaptGridAndVariables(indicator);
+ assembler->updateAfterGridAdaption();
}
// We solve the non-linear system with time step control.
nonLinearSolver.solve(x, *timeLoop);
- //We make the new solution the old solution.
+ // We make the new solution the old solution.
xOld = x;
gridVariables->advanceTimeStep();
- //We advance to the time loop to the next step.
+ // We advance to the time loop to the next step.
timeLoop->advanceTimeStep();
- //We write vtk output for each time step
+ // We write vtk output for each time step
vtkWriter.write(timeLoop->time());
- //We report statistics of this time step
+ // We report statistics of this time step
timeLoop->reportTimeStep();
- //We set a new dt as suggested by the newton solver for the next time step
+ // We set a new dt as suggested by the newton solver for the next time step
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
@@ -358,15 +328,6 @@ catch (const Dumux::ParameterException &e)
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
-catch (const Dune::DGFException & e)
-{
- std::cerr << "DGF exception thrown (" << e <<
- "). Most likely, the DGF file name is wrong "
- "or the DGF file is corrupted, "
- "e.g. missing hash at end of file or wrong number (dimensions) of entries."
- << " ---> Abort!" << std::endl;
- return 2;
-}
catch (const Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
diff --git a/examples/2pinfiltration/doc/main_intro.md b/examples/2pinfiltration/doc/main_intro.md
index cb7b9fef0fe04ee6f7155cb444a1fc79c1b4f5bf..c3fc1452fb00a77c432df157f14323a5cb2658ef 100644
--- a/examples/2pinfiltration/doc/main_intro.md
+++ b/examples/2pinfiltration/doc/main_intro.md
@@ -1,7 +1,7 @@
# Part 2: Main program flow
-This main program flow is implemented in the `main()` function
-of the program which is defined in the file `main.cc` described below. Additionally, we can see how to implement an adaptive grid with a saturation dependent indicator.
+This main program flow is implemented in the `main` function
+of the program which is defined in the file `main.cc` described below. Additionally, we can see how to implement an adaptive grid with a saturation-dependent indicator.
The code documentation is structured as follows:
diff --git a/examples/2pinfiltration/main.cc b/examples/2pinfiltration/main.cc
index 980d17a2ab767089a8248e9d20b2b719499bd7df..41db9c2cb599a2995180989e66f8e364bb166b11 100644
--- a/examples/2pinfiltration/main.cc
+++ b/examples/2pinfiltration/main.cc
@@ -53,7 +53,9 @@
//
// ### The main function
-// At the beginning of the simulation, we create a type tag with a suitable alias for our problem. // This will contain all the properties that are needed to define the problem set-up and the model we use. Additionally, we have to create an instance of `Dune::MPIHelper` and parse the run-time arguments:
+// At the beginning of the simulation, we create an alias for our problem type tag, for which we have
+// defined the properties of our simulation (see `properties.hh`). Additionally, we have to create an
+// instance of `Dune::MPIHelper` and parse the run-time arguments:
// [[codeblock]]
int main(int argc, char** argv) try
{
@@ -62,16 +64,17 @@ int main(int argc, char** argv) try
// we define the type tag for this problem
using TypeTag = Properties::TTag::PointSourceExample;
- // maybe initialize MPI and/or multithreading backend
+ // (maybe) initialize MPI and/or multithreading backend
Dumux::initialize(argc, argv);
const auto& mpiHelper = Dune::MPIHelper::instance();
- //We parse command line arguments and input file
+ // Construct parameter tree from command line arguments (and input file if given in the arguments)
Parameters::init(argc, argv);
// [[/codeblock]]
// #### Create the grid
- // A gridmanager tries to create the grid either from a grid file or the input file.
+ // The `gridManager` creates a grid from our parameter choices in the input file (which could be a grid
+ // file or specified domain dimensions and number of cells, as done in this example).
// [[codeblock]]
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();
@@ -83,22 +86,23 @@ int main(int argc, char** argv) try
// [[/codeblock]]
// #### Set-up of the problem
- //
- // 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.
+ // We build the finite volume geometry, which allows us to iterate over subcontrolvolumes (scv) and
+ // subcontrolvolume faces (scvf) embedded in the elements of the grid partition.
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
- // In the problem, we define the boundary and initial conditions and compute the point sources. The `computePointSourceMap` method is inherited from the fvproblem and therefore specified in the `dumux/common/fvproblem.hh`. It calls the `addPointSources` method specified in the `problem.hh` file.
+ // In the problem, we define the boundary and initial conditions and compute the point sources.
+ // The `computePointSourceMap` function is inherited from `FVProblem` (see `dumux/common/fvproblem.hh`).
+ // It calls the `addPointSources` method, which is empty per default, but can be overloaded in user problems.
+ // As we have seen, in this example we do specify a point source (see `problem.hh`).
// [[codeblock]]
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
- // We call the `computePointSourceMap` method to compute the point sources.
problem->computePointSourceMap();
// [[/codeblock]]
- // We initialize the solution vector and then use the solution vector to initialize the `gridVariables`.
+ // We initialize the solution vector (all primary variables on the grid) and then use it to initialize the
+ // `gridVariables`, which provide access to both primary & secondary variables (per sub-control volume).
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x;
problem->applyInitialSolution(x);
@@ -109,7 +113,8 @@ int main(int argc, char** argv) try
gridVariables->init(x);
// ##### Grid adaption
- // We instantiate the indicator for grid adaption & the data transfer, we read some parameters for indicator from the input file.
+ // The following piece of code shows how we instantiate an indicator class which determines where to adapt the grid,
+ // and a data transfer class that maps a solution to an adapted grid.
// [[codeblock]]
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const Scalar refineTol = getParam<Scalar>("Adaptive.RefineTolerance");
@@ -124,56 +129,39 @@ int main(int argc, char** argv) try
TwoPGridDataTransfer<TypeTag> dataTransfer(problem, gridGeometry, gridVariables, x);
// [[/codeblock]]
- // We do an initial refinement around sources/BCs. We use the `GridAdaptInitializationIndicator` defined in `dumux/adaptive/initializationindicator.hh` for that.
- GridAdaptInitializationIndicator<TypeTag> initIndicator(problem, gridGeometry, gridVariables);
-
- //We refine up to the maximum level. For every level, the indicator used for the refinement/coarsening is calculated. If any grid cells have to be adapted, the gridvariables and the pointsourcemap are updated.
- // [[codeblock]]
- const auto maxLevel = getParam<std::size_t>("Adaptive.MaxLevel", 0);
- for (std::size_t i = 0; i < maxLevel; ++i)
- {
- //we calculate the initial indicator for adaption for each grid cell using the initial solution x
- initIndicator.calculate(x);
-
- //and then we mark the elements that were adapted.
+ // Convenience functions to perform a grid adaptation given a pre-calculated indicator.
+ // [[codeblock]]
+ const auto adaptGridAndVariables = [&] (auto& indicator) {
+ // We mark and adapt the elements according to the indicator.
bool wasAdapted = false;
- if (markElements(gridManager.grid(), initIndicator))
+ if (markElements(gridManager.grid(), indicator))
wasAdapted = adapt(gridManager.grid(), dataTransfer);
- // In case of a grid adaptation, the gridvariables and the pointsourcemap are updated.
+ // In case of any adapted elements, the gridvariables and the pointsourcemap are updated.
if (wasAdapted)
{
// We overwrite the old solution with the new (resized & interpolated) one
xOld = x;
- //We initialize the secondary variables to the new (and "new old") solution
+ // We initialize the secondary variables to the new (and "new old") solution
gridVariables->updateAfterGridAdaption(x);
// we update the point source map after adaption
problem->computePointSourceMap();
}
- }
+ };
// [[/codeblock]]
- // Depending on the initial conditions, another grid adaptation might be necessary.
- //The gridadaptindicator uses the input parameters `Adaptive.RefineTolerance` and `Adaptive.CoarsenTolerance` for this step.
- //Again, if elements were adapted, we mark them and update gridvariables and the pointsourcemap accordingly.
- // [[codeblock]]
- indicator.calculate(x, refineTol, coarsenTol);
-
- //we mark the elements that were adapted
- bool wasAdapted = false;
- if (markElements(gridManager.grid(), indicator))
- wasAdapted = adapt(gridManager.grid(), dataTransfer);
-
- // In case of an additional grid adaptation, the gridvariables and the pointsourcemap are updated again.
- if (wasAdapted)
+ // We do initial refinements around sources/BCs, for which we use the `GridAdaptInitializationIndicator` defined in `dumux/adaptive/initializationindicator.hh`.
+ // Afterwards, depending on the initial conditions, another grid adaptation using our `indicator` above might be necessary, which uses
+ // `Adaptive.RefineTolerance` and `Adaptive.CoarsenTolerance` for this step.
+ GridAdaptInitializationIndicator<TypeTag> initIndicator(problem, gridGeometry, gridVariables);
+ const auto maxLevel = getParam<std::size_t>("Adaptive.MaxLevel", 0);
+ for (std::size_t i = 0; i < maxLevel; ++i)
{
- // Overwrite the old solution with the new (resized & interpolated) one
- xOld = x;
- // Initialize the secondary variables to the new (and "new old") solution
- gridVariables->updateAfterGridAdaption(x);
- // Update the point source map
- problem->computePointSourceMap();
+ initIndicator.calculate(x);
+ adaptGridAndVariables(initIndicator);
}
+ indicator.calculate(x, refineTol, coarsenTol);
+ adaptGridAndVariables(indicator);
// [[/codeblock]]
// #### Solving the problem
@@ -214,50 +202,33 @@ int main(int argc, char** argv) try
timeLoop->start(); do
{
// We only want to refine/coarsen after first time step is finished, not before.
- //The initial refinement was already done before the start of the time loop.
- //This means we only refine when the time is greater than 0.
+ // The initial refinement was already done before the start of the time loop.
+ // This means we only refine when the time is greater than 0. Note that we now also
+ // have to update the assembler, since the sizes of the residual vector and jacobian matrix change.
if (timeLoop->time() > 0)
{
- // again we compute the refinement indicator with the `TwoPGridAdaptIndicator`
indicator.calculate(x, refineTol, coarsenTol);
-
- //we mark elements and adapt grid if necessary
- wasAdapted = false;
- if (markElements(gridManager.grid(), indicator))
- wasAdapted = adapt(gridManager.grid(), dataTransfer);
-
- //In case of a grid adaptation, the gridvariables and the pointsourcemap are updated again.
- if (wasAdapted)
- {
- // We overwrite the old solution with the new (resized & interpolated) one
- xOld = x;
- // We initialize the secondary variables to the new (and "new old") solution
- gridVariables->updateAfterGridAdaption(x);
- // We also need to update the assembler (sizes of residual and Jacobian change)
- assembler->updateAfterGridAdaption();
- // We update the point source map
- problem->computePointSourceMap();
- }
- // we leave the refinement step
+ adaptGridAndVariables(indicator);
+ assembler->updateAfterGridAdaption();
}
// We solve the non-linear system with time step control.
nonLinearSolver.solve(x, *timeLoop);
- //We make the new solution the old solution.
+ // We make the new solution the old solution.
xOld = x;
gridVariables->advanceTimeStep();
- //We advance to the time loop to the next step.
+ // We advance to the time loop to the next step.
timeLoop->advanceTimeStep();
- //We write vtk output for each time step
+ // We write vtk output for each time step
vtkWriter.write(timeLoop->time());
- //We report statistics of this time step
+ // We report statistics of this time step
timeLoop->reportTimeStep();
- //We set a new dt as suggested by the newton solver for the next time step
+ // We set a new dt as suggested by the newton solver for the next time step
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
@@ -265,7 +236,7 @@ int main(int argc, char** argv) try
// The following piece of code prints a final status report of the time loop
// before the program is terminated and we print he dumux end message
- // [[codeblock]]
+ // [[codeblock]]
timeLoop->finalize(leafGridView.comm());
// print dumux end message
@@ -286,15 +257,6 @@ catch (const Dumux::ParameterException &e)
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
-catch (const Dune::DGFException & e)
-{
- std::cerr << "DGF exception thrown (" << e <<
- "). Most likely, the DGF file name is wrong "
- "or the DGF file is corrupted, "
- "e.g. missing hash at end of file or wrong number (dimensions) of entries."
- << " ---> Abort!" << std::endl;
- return 2;
-}
catch (const Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
diff --git a/examples/2pinfiltration/problem.hh b/examples/2pinfiltration/problem.hh
index ac67777942e2a162151bd48c709b4694cc8d1101..3c47952304ec385382119b246729f323999026e1 100644
--- a/examples/2pinfiltration/problem.hh
+++ b/examples/2pinfiltration/problem.hh
@@ -64,10 +64,11 @@ public:
}
// #### Boundary types
- // We define the type of boundary conditions depending on location. Two types of boundary conditions
- // can be specified: Dirichlet or Neumann boundary condition. On a Dirichlet boundary, the values of the
+ // We define the type of boundary conditions depending on the location. Two types of boundary conditions
+ // can be specified: Dirichlet or Neumann boundary conditions. On a Dirichlet boundary, the values of the
// primary variables need to be fixed. On a Neumann boundary condition, values for derivatives need to be fixed.
- // Mixed boundary conditions (different types for different equations on the same boundary) are not accepted.
+ // Depending on the chosen discretization scheme, mixed boundary conditions (different types for different equations
+ // on the same boundary) may not be accepted.
// [[codeblock]]
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
diff --git a/examples/2pinfiltration/properties.hh b/examples/2pinfiltration/properties.hh
index 05198d8fc9940359f02af5887ef92e5325f6d271..00d2f93d646c50997fcc5d9f2178d434bbb8609e 100644
--- a/examples/2pinfiltration/properties.hh
+++ b/examples/2pinfiltration/properties.hh
@@ -101,7 +101,9 @@ struct FluidSystem<TypeTag, TTag::PointSourceExample>
// The two-phase model implements different primary variable formulations.
// Here we choose the pressure of the first phase and the saturation of the second phase.
// The order of phases is specified by the fluid system.
-// In this case that means that the primary variables are water pressure and DNAPL saturation.
+// With our choice of FluidSystem this means that the first phase is water and the second
+// phase is DNAPL. We choose the `p0s1` formulation, which means a $`p_w - S_n`$ formulation,
+// as the first phase (i.e. water) is the wetting phase in this example.
template<class TypeTag>
struct Formulation<TypeTag, TTag::PointSourceExample>
{ static constexpr auto value = TwoPFormulation::p0s1; };
diff --git a/examples/2pinfiltration/spatialparams.hh b/examples/2pinfiltration/spatialparams.hh
index 965ccc84488236806521c0466a054f79111d7448..cef70502a226255c0dc4076e5dcec4543674ef3f 100644
--- a/examples/2pinfiltration/spatialparams.hh
+++ b/examples/2pinfiltration/spatialparams.hh
@@ -94,7 +94,8 @@ public:
}
- // Here we can define which phase is to be considered as the wetting phase. Here the wetting phase is the the first phase of the fluidsystem. In this case that is water.
+ // Here we can define which phase is to be considered as the wetting phase. We choose the first phase
+ // in the fluid system, which is water in our case (see definition of `FluidSystem` property).
template<class FluidSystem>
int wettingPhaseAtPos(const GlobalPosition& globalPos) const
{ return FluidSystem::phase0Idx; }