From f2574821820d997828c59e60bb8e4434b3c733e4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= <dennis.glaeser@iws.uni-stuttgart.de>
Date: Fri, 20 Oct 2023 11:47:36 +0200
Subject: [PATCH] [examples] regenerate docs
---
examples/2pinfiltration/doc/2p.md | 14 +--
examples/2pinfiltration/doc/main.md | 133 ++++++++++------------------
2 files changed, 56 insertions(+), 91 deletions(-)
diff --git a/examples/2pinfiltration/doc/2p.md b/examples/2pinfiltration/doc/2p.md
index c7e5e94a61..5f522154db 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/main.md b/examples/2pinfiltration/doc/main.md
index ce6354fc08..468bbb09c7 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;
--
GitLab