README.md

# Exercise Mainfiles (DuMuX course)
<br>
## Problem set-up

This exercise will make you familiar the program sequence in DuMux and how different levels of complexity can be realized in the main file according to the complexity of your physical problem.

In order to do so, there are three examples of one phase flow problems. Two examples (a and b) are stationary problems and the third example (c) is an instationary problem.

The stationary examples differ in the fluidssystems they are using which means they differ in the fluid properties (e.g. density, thermal conductivity etc). The first problem (a) uses an incompressible fluid which means that the density does not change when pressure changes. This makes it possible to solve the system linearly. The second problem uses a compressible fluid, that means the density is a function of pressure and we need to use a nonlinear solver.

To summarize the problems differ in:
* exercise_1p_a: a one-phase incompressible, stationary problem
* exercise_1p_b: a one-phase compressible, stationary problem
* exercise_1p_c: a one-phase compressible, instationary problem

The problem set-up for all three examples is always the same: It is a two dimensional problem and the domain is $`1 m`$ by $`1 m`$. It is a heterogeneous set-up with a lens in the middle of the domain which has a lower permeability ($`1\cdot 10^{-12} m^2`$ compared to  $`1\cdot 10^{-10} m^2`$ in the rest of the domain).

<img src="https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/raw/master/exercises/extradoc/exercise1_1p_setup.png" width="1000">

In the beginning there is a uniform pressure of $`1\cdot 10^5 Pa`$ in the whole domain. On the top and the bottom border dirichlet boundary conditions are set with a pressure of  $`1\cdot 10^5 Pa`$ on top and  $`2 \cdot 10^5 Pa`$ on the bottom. At the sides there is no in- or outflow and there are no source terms.

## Preparing the exercise

* Navigate to the directory `dumux-course/exercises/exercise-mainfile`

<br><br>
### Task 1: Getting familiar with the code
<hr>

Locate all the files you will need for this exercise
* The __main file__ for the __1p incompressible, stationary__ problem : `exercise_1p_a.cc`
* The __main file__ for the __1p compressible, stationary__ problem : `exercise_1p_b.cc`
* The __main file__ for the __1p compressible, instationary__ problem : `exercise_1p_c.cc`
* The shared __problem file__: `1pproblem.hh`
* The shared __spatial parameters file__: `1pspatialparams.hh`
* The __input file__ for the __1p incompressible, stationary__ problem: `exercise_1p_a.input`
* The __input file__ for the __1p compressible, stationary__ problem: `exercise_1p_b.input`
* The __input file__ for the __1p compressible, instationary__ problem: `exercise_1p_c.input`

Please pay special attention to the similarities and differences in the three main files. The first main file is solved linearly and does not need a newton solver or any other nonlinear solver method. The second problem is a nonlinear problem and uses newton's method to solve the system. The third problem is nonlinear and additionally instationary. Therefore a time loop needs to be included in the main file.

The general structure of any main file in DuMux is:

* the specific problem TypeTag is defined for the problem. This example shows the TypeTag for the CompressibleProblem

```c++
// define the type tag for this problem
using TypeTag = TTAG(OnePCompressible);
```
The TypeTag is created in the `1pproblem.hh`. There you can see that it inherits from the __OneP__ and additionally from the __CCTpfaModel__ which defines the discretization method, which is in this case the cell-centered tpfa method.

* a gridmanager tries to create the grid either from a grid file or the input file

```c++
GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
gridManager.init();
```
* we create the finite volume grid geometry, the problem, solutionvector and the gridvariables and initialize them. Additionally we initialize the vtkoutput. Each model has a predefined model specific output with relevant parameters for that model.

```c++
// create the finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
fvGridGeometry->update();

// the problem (initial and boundary conditions)
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
auto problem = std::make_shared<Problem>(fvGridGeometry);

// the solution vector
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
SolutionVector x(fvGridGeometry->numDofs());

// the grid variables
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
gridVariables->init(x);

// initialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
vtkWriter.write(0.0);
```

* then we need to assemble and solve the system. Depending on the problem this can be done with a linear solver or with a nonlinear solver. If the problem is time dependent we additionally need a time loop. An example for that is given in exercise_1p_c:

```c++
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");

// instantiate time loop
auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);

// the assembler with time loop for instationary problem
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);

// the linear solver
using LinearSolver = ILU0BiCGSTABBackend;
auto linearSolver = std::make_shared<LinearSolver>();

// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);

// set some check points for the time loop
timeLoop->setPeriodicCheckPoint(tEnd/10.0);

// time loop
timeLoop->start(); do
{
    // set previous solution for storage evaluations
    assembler->setPreviousSolution(xOld);

    // linearize & solve
    nonLinearSolver.solve(x, *timeLoop);

    // make the new solution the old solution
    xOld = x;
    gridVariables->advanceTimeStep();

    // advance to the time loop to the next step
    timeLoop->advanceTimeStep();

    // write vtk output
    if (timeLoop->isCheckPoint())
        vtkWriter.write(timeLoop->time());

    // report statistics of this time step
    timeLoop->reportTimeStep();

    // set new dt as suggested by the newton solver
    timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));

} while (!timeLoop->finished());

timeLoop->finalize(leafGridView.comm());
```

<hr><br><br>
### Task 2: Compiling and running an executable
<hr>

* Change to the build-directory

```bash
cd ../../build-cmake/exercises/exercise-mainfile
```

* Compile all three executables `exercise_1p_a` and `exercise_1p_b` and `exercise_1p_c`

```bash
make exercise_1p_a exercise_1p_b exercise_1p_c
```

* Execute the three problems and inspect the result

```bash
./exercise_1p_a
./exercise_1p_b
./exercise_1p_c
```

* you can look at the results (e.g. for the first example) with paraview:

```bash
paraview 1p_incompressible_stationary.pvd
```

<hr><br><br>
### Task 3: Analytical differentiation
<hr>

In the input file `exercise_1p_a.input` you will see that there is a variable BaseEpsilon. This defines the base for the epsilon used in the numeric differentiation. If that value is too small, you will see that the solution of the numeric differentiation is not correct. Change that value to $`1 \cdot 10^{-15}`$ and have a look at the solution.

For the incompressible one phase problem it is possible to also have an analytic solution method. Then the epsilon does not play a role anymore since the derivatives are calculated analytically. To implement that follow the directions in the `exercise_1p_a.cc` and the `1pproblem.hh` marked by:

```c++
// TODO: dumux-course-task
```
For the analytic solution of your immiscible problem you need analytic solutions for the derivatives of the jacobian. For that we have a special local residual, the `OnePIncompressibleLocalResidual` which provides that. You just need to include that in your `1pproblem.hh` and use that instead of the `immisciblelocalresidual.hh` which is used as a standard for all immiscible models.

Additionally you need to set the differentiation method in the main file `exercise_1p_a.cc` to analytic.