README.md

# Exercise Discrete Fractures (DuMuX Course)

In this exercise we are going to use the _DuMuX multidomain_ framework, and in particular the _facet coupling_ module. This allows you to define coupled problems, in which a (d-1)-dimensional domain lives on the element facets of a d-dimensional domain and the coupling occurs via the transfer fluxes between them. Here, we are going to make use of this in the context of fractured porous media, where the fractures are represented as lower-dimensional entities embedded in a porous matrix. In the _facet coupling_ framework it is required that the grid resolves the fractures by aligning the matrix element facets with the fracture geometries, for which we are using the open-source mesh generator _Gmsh_ ([see _Gmsh_ documentation](http://gmsh.info/doc/texinfo/gmsh.html)). Note that the Gmsh file format (_.msh_) is currently the only supported file format within the _facet coupling_ framework.

## Problem set-up

We consider a domain of 100 x 100 m with a set of ten fractures. This geometry has been taken from [Flemisch et al. (2018)](https://www.sciencedirect.com/science/article/pii/S0309170817300143), but was slightly modified to better fit to the purposes of this exercise. Apart from one fracture extending up to the top boundary, all other fractures are immersed, meaning that the fracture tips lie within the domain (see image below).

![](../extradoc/exercisefractures_setup.png)

In the initial setup of this exercise, we want to consider buoyancy-driven upwards migration of nitrogen (gas) in a fully water-saturated medium. Thus, as initial conditions we use hydrostatic pressure conditions and a zero nitrogen saturation. On the boundaries we apply Neumann no-flow boundary conditions on the left and right sides, while using Dirichlet boundary conditions on the top and the bottom of the domain. We set the initial conditions as boundary conditions on the Dirichlet segments, except for the middle part of the lower boundary ($`y = 0 \wedge x > 25 \wedge x < 75`$), where we set a non-zero nitrogen saturation (settable via the input file). Through this segment of the lower boundary, nitrogen can intrude the domain and flow upwards driven by buoyancy.

## Preparing the exercise

Navigate to the directory `dumux/exercises/exercise-fractures` and familiarize yourself with the files necessary for this exercise:

* The __main file__: `main.cc`
* The __input file__: `params.input`
* The __problem file__ for the matrix domain: `matrixproblem.hh`
* The __problem file__ for the fracture domain: `fractureproblem.hh`
* The __spatial parameters file__ for the matrix domain: `matrixspatialparams.hh`
* The __spatial parameters file__ for the matrix domain: `fracturespatialparams.hh`
* The __properties file__ for both domains: `properties.hh`

If you want to learn more about the program flow of coupled problems within the _DuMuX multidomain_ framework, take a closer look at the __main file__. Therein, you will find extensive comments above the key instructions that will help you to identify the main differences with respect to an uncoupled _DuMuX_ model. We will now briefly guide you through this, however, this is not necessary for this exercise. Feel free to continue with the section __Running the program__ and come back to this point later.

The basic idea of the _multidomain_ framework is as follows: you define the sub-problems individually, while all data and algorithms with respect to the coupling between the domains is contained in the `CouplingManager` object. This object is chosen depending on the coupling mode, thus, here we are using the implementation for models considering coupling across the element facets, i.e. the class `Dumux::FacetCouplingManager<...>` (see main file). As mentioned earlier, coupling occurs via the transfer fluxes between the d-dimensional matrix domain and the (d-1)-dimensional fracture domain. These fluxes act as additional source/sink terms for the fracture domain, which have to be implemented in the respective `source(...)` function (see `fractureproblem.hh`):

```cpp
//! Evaluate the source term at a given position
NumEqVector source(const Element& element,
                   const FVElementGeometry& fvGeometry,
                   const ElementVolumeVariables& elemVolVars,
                   const SubControlVolume& scv) const
{
    // evaluate sources from bulk domain using the function in the coupling manager
    auto source = couplingManagerPtr_->evalSourcesFromBulk(element, fvGeometry, elemVolVars, scv);

    // these sources are in kg/s, divide by volume and extrusion to have it in kg/s/m³
    source /= scv.volume()*elemVolVars[scv].extrusionFactor();
    return source;
}
```

As you can see, we call the convenience function `evalSourcesFromBulk(...)` available in the _facet coupling_ manager, which requires us to have access to it at this point. This access is provided in the __main file__ by handing over a shared pointer to the sub-problems, i.e.:

```cpp
// the coupling manager (needs the coupling mapper)
auto couplingManager = std::make_shared<TheCouplingManager>();
couplingManager->init(matrixProblem, fractureProblem, couplingMapper, x);

// we have to set coupling manager pointer in sub-problems
// they also have to be made accessible in them (see e.g. matrixproblem.hh)
matrixProblem->setCouplingManager(couplingManager);
fractureProblem->setCouplingManager(couplingManager);
```

From the matrix side, the coupling works a bit different. Since the fracture domain lives on the matrix domain's facets, the state inside the fracture (e.g. pressure, temperature etc.) will affect the fluxes computed across the respective element facets. However, flux assembly depends on the underlying finite-volume scheme and is not convenient to implement. We therefore created _TypeTag_ nodes for the matrix sub-problems in models considering facet coupling. By using these _TypeTags_, all necessary modifications for the flux computation on interior boundaries is done for you. For instance, in `properties.hh`, we define the _TypeTag_ for the matrix problem as follows (please read the provided comments):

```cpp
// We are using the framework for models that consider coupling
// across the element facets of the bulk domain. This has some
// properties defined, which we have to inherit here. In this
// exercise we want to use a cell-centered finite volume scheme
// with tpfa.
#include <dumux/multidomain/facet/cellcentered/tpfa/properties.hh>

namespace Dumux::Properties {

// create the type tag node for the matrix and fracture sub-problems
namespace TTag {
struct MatrixProblem { using InheritsFrom = std::tuple<CCTpfaFacetCouplingModel, TwoP>; };
struct FractureProblem { using InheritsFrom = std::tuple<TwoP, CCTpfaModel>; };
} // end namespace TTag

```

Additionally, we need to provide access to the coupling manager in the matrix problem, so that the flux assembly engine of the matrix domain has access to the state inside the fracture (which the fluxes depend on):

```cpp
//! returns reference to the coupling manager.
const CouplingManager& couplingManager() const
{ return *couplingManagerPtr_; }
```

At this point, we are done with the extra work involved in defining the sub-problems. There are two more differences with respect to the main program flow of uncoupled problems: we now have to use the multidomain assembler for the monolithic assembly of the coupled system of equations (see __main file__)

```cpp
// the assembler for the coupled problem
using Assembler = MultiDomainFVAssembler<TheMultiDomainTraits, TheCouplingManager, DiffMethod::numeric, /*implicit?*/true>;
auto assembler = std::make_shared<Assembler>( std::make_tuple(matrixProblem, fractureProblem),
                                              std::make_tuple(matrixFvGridGeometry, fractureFvGridGeometry),
                                              std::make_tuple(matrixGridVariables, fractureGridVariables),
                                              couplingManager,
                                              timeLoop);
```

and a specialized implementation for Newton's method as non-linear solver:

```cpp
// the non-linear solver
using NewtonSolver = Dumux::MultiDomainNewtonSolver<Assembler, LinearSolver, TheCouplingManager>;
auto newtonSolver = std::make_shared<NewtonSolver>(assembler, linearSolver, couplingManager);
```

## Running the program

Head to the build directory, compile the exercise and execute the program by typing

```bash
cd build-cmake/exercises/exercise-fractures
make exercise_fractures
./exercise_fractures params.input
```

Take a look at the results by opening the files `matrix.pvd` and `fractures.pvd` with _Paraview_. In order to increase the visibility of the fracture solution, you might want to apply the tube filter to it (e.g. with a tube radius of 0.2). The result should look like this:

![](../extradoc/exercisefractures_initsol.png)

You can see how the fractures act as preferential flowpaths in the upwards movement of nitrogen due to their higher permeabilities. Additionally, you can observe comparatively high nitrogen saturations in the fracture tips as a result of the fractures acting as capillary traps due to the lower capillary pressures inside them. Consider turning them into capillary barriers, e.g. by setting __SpatialParams.VGAlpha = 1e-5__ in the __Fracture__ group in the input file (`params.input`). The nitrogen saturations in the fractures (especially in the fracture tips) should now be lower than in the surrounding matrix. Do not forget to reset __SpatialParams.VGAlpha = 1e-1__ in the input file after you are done.

## Task A: Change the boundary conditions and disable gravity

In order for the influence of the fractures to be more visible in the resulting pressure fields, __switch off gravity__ in the input file and change the boundary conditions in the two sub-problems (i.e. modify the function `boundaryTypesAtPos(...)` in `matrixproblem.hh` and `fractureproblem.hh`). We want to define the boundary conditions such that Neumann no-flow boundaries are used everywhere in the fracture domain and everywhere in the matrix domain except for the lower left and the upper right part of the lateral sides (i.e. $`x = 0 \wedge y < 25`$ and $`x = 100 \wedge y > 75`$).

Furthermore, modify the function `dirichletAtPos(...)` in `matrixproblem.hh` such that an overpressure is added to the initial pressure and a non-zero saturation is applied on the right Dirichlet segment ($`x = 100 \wedge y > 75`$). For this you can use the private variables `boundaryOverPressure_` and `boundarySaturation_` which are stored in the matrix problem and are read from the input file. Compile and execute the program again. The resulting water pressure distribution should look like this:

![](../extradoc/exercisefractures_a.png)

## Task B: Turn the fractures into barriers

In this part, we want to change the parameterization of the fractures such that they act as both hydraulic and capillary barriers. In the `fracturespatialparams.hh` file you will see that the spatial parameters class already contains a set of parameters with the post-fix _Barrier_, which are read in from the input file in the class constructor.

Use these as return values for porosity, permeability and material law parameters. Take a look at the results.

You will see that only little nitrogen enters the domain, but if you display the pressure distribution in the matrix, you will notice that we cannot observe the pressure distribution we would expect. With the fracture permeabilities being this much lower than the matrix permeability, we would expect substantial pressure drops to be visible across them.

This has to do with the chosen coupling condition. The assembly of the transfer fluxes between fracture and matrix across the matrix element facets occurs on the basis of interface conditions, which type can be set in the function `interiorBoundaryTypes(...)`. Note that this function has to be implemented by the matrix problem. Take a look at this function in the `matrixproblem.hh` file and find the following comment:

```cpp
//! Specifies the type of interior boundary condition
BoundaryTypes interiorBoundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const
{
    BoundaryTypes values;

    // Here we set the type of condition to be used on faces that coincide
    // with a fracture. If Neumann is specified, a flux continuity condition
    // on the basis of the normal fracture permeability is evaluated. If this
    // permeability is lower than that of the matrix, this approach is able to
    // represent the resulting pressure jump across the fracture. If Dirichlet is set,
    // the pressure jump across the fracture is neglected and the pressure inside
    // the fracture is directly applied at the interface between fracture and matrix.
    // This assumption is justified for highly-permeable fractures, but lead to erroneous
    // results for low-permeable fractures.
    // Here, we consider "open" fractures for which we cannot define a normal permeability
    // and for which the pressure jump across the fracture is neglectable. Thus, we set
    // the interior boundary conditions to Dirichlet.
    // IMPORTANT: Note that you will never be asked to set any values at the interior boundaries!
    //            This simply chooses a different interface condition!
    // TODO dumux-course-task C
    // Change coupling conditions!
    values.setAllDirichlet();

    return values;
}
```

To summarize: when using Neumann-type interior boundary conditions, flux continuity conditions are evaluated on the two interfaces (the two sides of the fracture) between fracture and matrix ([see e.g. Martin et al. (2005)](https://link.springer.com/article/10.1007/s10596-012-9302-6)). Note that these interfaces geometrically coincide but allow for both different pressures and fluxes to prevail on the two sides of the fracture. When using Dirichlet-type interior boundary conditions, the pressure in the fracture is assumed to be invariant in normal direction of the fracture and is directly applied at the interfaces. While the fluxes still can be different on the two interfaces, the pressure is now the same. Thus, this is only valid for fractures that are highly permeable in normal direction (note that "open" fractures have infinite normal permeability).

This means we have to use Neumann-type interior boundary conditions here. Change this accordingly, compile and rerun the exercise again. You should now be able to see jumps in pressure in the matrix domain with an especially prominent one across the first vertical fracture (see image below).

![](../extradoc/exercisefractures_b.png)

### Task C: Define both open fractures and barriers

We have seen in the previous task what type of coupling or interior boundary conditions we have to use to describe open fractures and barriers. In this task, we want to define two fractures as barriers while considering the remaining fractures as "open". For this we want to make use of the domain markers available in gmsh. In the file `grids/complex.geo`, you will find that two different __physical indices__ were given to the different fracture segments:

```
// conductive fractures get physical index 1
Physical Line(1) = {6,7,9,11,12,13,14,15,16,17,18,19,20,21,23,25};

// blocking fractures get physical index 2
Physical Line(2) = {8,10,22,24,26};
```

We now want to give all fracture elements that are tagged with domain marker 2 the properties for barriers and assign the parameters for open fractures to all other elements. The domain markers are read from the grid file and can be obtained from the `GridManager` subsequent to grid creation. The main file of this exercise has been implemented such that these markers are already passed to the spatial parameters of the sub-domains (see `main.cc`, search for `// pass the model parameter group to the spatial params`). In the spatial parameters you will find a convenience function which allows you to obtain the domain marker for a given element (see `fracturespatialparams.hh`):

```cpp
//! returns the domain marker for an element
int getElementDomainMarker(const Element& element) const
{ return gridDataPtr_->getElementDomainMarker(element); }
```

The domain markers are also already added to the output (see `main.cc`, search for `// add domain markers to output`). To visualize them, open any of your previously produced results with _Paraview_ and take a look at them by selecting __domainMarker__.

Adjust the functions for permeability, porosity and material law parameters in the `fracturespatialparams.hh` file such that they are selected depending on the domain marker of the elements. You will see in the results that the pressure jump across the first vertical fracture is now lower than before, because there are highly permeable fractures crossing it, allowing for a pressure release into the other parts of the domain.

## Additional task:

We can now also set the interior boundary condition depending on the domain marker of the fracture domain elements. What we would like to do is using Dirichlet-type coupling conditions for open fractures and Neumann-type coupling conditions for barriers. This means we need to obtain the domain marker of the fracture domain element to which a matrix element facet couples. Consider using the following piece of code in your `interiorBoundaryTypes(...)` function in `matrixproblem.hh` in order to realize this:

```cpp
// we need to obtain the domain marker of the fracture element that is coupled to this face
// therefore we first get the fracture problem from the coupling manager. For this test, we
// know that the fracture domain id is 1 (see main file)
static constexpr auto fractureDomainId = Dune::index_constant<1>();
const auto& fractureProblem = couplingManager().problem(fractureDomainId);

// use helper function in coupling manager to obtain the element this face couples to
const auto fractureElement = couplingManager().getLowDimElement(element, scvf);

// obtain marker from the spatial params of the fracture problem
const auto fractureElementMarker = fractureProblem.spatialParams().getElementDomainMarker(fractureElement);
```

Note that this will not have a visible effect on the results, because the permeability chosen for the open fractures is very high, leading to identical results for the two approaches. However, as mentioned above, the use of interior Neumann-type boundary conditions involves the evaluation of flux continuity conditions at the interfaces between fracture and matrix on the basis of the fracture normal permeability. Here, we are defining scalar permeabilities on the fracture with the result that the normal and tangential permeabilities are the same (if you want them to be different, you have to define tensorial permeabilities). The high value for the "open" fractures leads to neglectable pressure jumps and seems to produce the right results. But, keep in mind that from a physical perspective we are using the wrong value as "open" fractures have infinite normal permeability. Using Dirichlet-type interior boundary conditions at the interfaces to "open" fractures is clearly the better choice.