Dennis Gläser committed Jul 19, 2018 1 ``````# Exercise Discrete Fractures (DuMuX Course) `````` Dennis Gläser committed Jul 19, 2018 2 `````` `````` Dennis Gläser committed Jul 19, 2018 3 ``````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. `````` Dennis Gläser committed Jul 19, 2018 4 5 6 `````` ## Problem set-up `````` Dennis Gläser committed Jul 19, 2018 7 ``````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). `````` Dennis Gläser committed Jul 19, 2018 8 `````` `````` Dennis Gläser committed Jul 19, 2018 9 ``````![](../extradoc/exercisefractures_setup.png) `````` Dennis Gläser committed Jul 19, 2018 10 `````` `````` Dennis Gläser committed Jul 19, 2018 11 ``````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. `````` Dennis Gläser committed Jul 19, 2018 12 13 14 `````` ## Preparing the exercise `````` Dennis Gläser committed Jul 19, 2018 15 ``````Navigate to the directory `dumux/exercises/exercise-fractures` and familiarize yourself with the files necessary for this exercise: `````` Dennis Gläser committed Jul 19, 2018 16 `````` `````` Katharina Heck committed Oct 02, 2019 17 18 ``````* The __main file__: `main.cc` * The __input file__: `params.input` `````` Dennis Gläser committed Jul 19, 2018 19 20 21 22 ``````* 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` `````` Maziar Veyskarami committed Apr 09, 2020 23 ``````* The __properties file__ for both domains: `properties.hh` `````` Dennis Gläser committed Jul 19, 2018 24 `````` `````` Dennis Gläser committed Jul 19, 2018 25 ``````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. `````` Dennis Gläser committed Jul 19, 2018 26 `````` `````` Dennis Gläser committed Jul 19, 2018 27 ``````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`): `````` Dennis Gläser committed Jul 19, 2018 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 `````` ```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; } ``` `````` Dennis Gläser committed Jul 19, 2018 45 ``````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.: `````` Dennis Gläser committed Jul 19, 2018 46 47 48 49 50 51 52 53 54 55 56 57 `````` ```cpp // the coupling manager (needs the coupling mapper) auto couplingManager = std::make_shared(); 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); ``` `````` Maziar Veyskarami committed Apr 09, 2020 58 ``````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): `````` Dennis Gläser committed Jul 19, 2018 59 60 61 62 63 64 65 66 67 `````` ```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 `````` Maziar Veyskarami committed Apr 09, 2020 68 ``````namespace Dumux::Properties { `````` Dennis Gläser committed Jul 19, 2018 69 `````` `````` Maziar Veyskarami committed Apr 09, 2020 70 ``````// create the type tag node for the matrix and fracture sub-problems `````` Theresa Schollenberger committed Dec 17, 2018 71 ``````namespace TTag { `````` Felix Weinhardt committed Dec 20, 2018 72 ``````struct MatrixProblem { using InheritsFrom = std::tuple; }; `````` Maziar Veyskarami committed Apr 09, 2020 73 ``````struct FractureProblem { using InheritsFrom = std::tuple; }; `````` Theresa Schollenberger committed Dec 17, 2018 74 ``````} // end namespace TTag `````` Maziar Veyskarami committed Apr 09, 2020 75 `````` `````` Dennis Gläser committed Jul 19, 2018 76 77 ````````` `````` Dennis Gläser committed Jul 19, 2018 78 ``````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): `````` Dennis Gläser committed Jul 19, 2018 79 80 81 82 83 84 85 `````` ```cpp //! returns reference to the coupling manager. const CouplingManager& couplingManager() const { return *couplingManagerPtr_; } ``` `````` Dennis Gläser committed Jul 19, 2018 86 ``````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__) `````` Dennis Gläser committed Jul 19, 2018 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 `````` ```cpp // the assembler for the coupled problem using Assembler = MultiDomainFVAssembler; auto assembler = std::make_shared( 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; auto newtonSolver = std::make_shared(assembler, linearSolver, couplingManager); ``` `````` Dennis Gläser committed Jul 19, 2018 106 107 108 109 110 111 112 ``````## 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 `````` Katharina Heck committed Oct 02, 2019 113 ``````./exercise_fractures params.input `````` Dennis Gläser committed Jul 19, 2018 114 115 ````````` `````` Dennis Gläser committed Jul 19, 2018 116 117 ``````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: `````` Dennis Gläser committed Jul 19, 2018 118 ``````![](../extradoc/exercisefractures_initsol.png) `````` Dennis Gläser committed Jul 19, 2018 119 `````` `````` Katharina Heck committed Oct 02, 2019 120 ``````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. `````` Dennis Gläser committed Jul 19, 2018 121 122 123 `````` ## Task A: Change the boundary conditions and disable gravity `````` Simon Scholz committed Jul 19, 2018 124 ``````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`\$). `````` Dennis Gläser committed Jul 19, 2018 125 `````` `````` Dennis Gläser committed Jul 19, 2018 126 ``````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: `````` Dennis Gläser committed Jul 19, 2018 127 `````` `````` Dennis Gläser committed Jul 19, 2018 128 ``````![](../extradoc/exercisefractures_a.png) `````` Dennis Gläser committed Jul 19, 2018 129 130 131 `````` ## Task B: Turn the fractures into barriers `````` Simon Scholz committed Jul 19, 2018 132 ``````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. `````` Dennis Gläser committed Jul 19, 2018 133 `````` `````` Simon Scholz committed Jul 19, 2018 134 135 136 ``````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. `````` Dennis Gläser committed Jul 19, 2018 137 138 139 140 141 142 143 144 145 146 147 148 149 `````` 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 `````` Dennis Gläser committed Jul 19, 2018 150 `````` // represent the resulting pressure jump across the fracture. If Dirichlet is set, `````` Dennis Gläser committed Jul 19, 2018 151 152 153 154 155 156 157 158 159 `````` // 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! `````` Farid Mohammadi committed Mar 31, 2020 160 161 `````` // TODO dumux-course-task C // Change coupling conditions! `````` Dennis Gläser committed Jul 19, 2018 162 163 164 165 166 167 `````` values.setAllDirichlet(); return values; } ``` `````` Dennis Gläser committed Jul 19, 2018 168 169 170 ``````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). `````` Dennis Gläser committed Jul 19, 2018 171 `````` `````` Dennis Gläser committed Jul 19, 2018 172 ``````![](../extradoc/exercisefractures_b.png) `````` Dennis Gläser committed Jul 19, 2018 173 174 175 `````` ### Task C: Define both open fractures and barriers `````` Dennis Gläser committed Jul 19, 2018 176 ``````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: `````` Dennis Gläser committed Jul 19, 2018 177 178 179 180 181 182 183 184 185 `````` ``` // 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}; ``` `````` Farid Mohammadi committed Mar 31, 2020 186 ``````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`): `````` Dennis Gläser committed Jul 19, 2018 187 188 189 190 191 192 193 `````` ```cpp //! returns the domain marker for an element int getElementDomainMarker(const Element& element) const { return gridDataPtr_->getElementDomainMarker(element); } ``` `````` Farid Mohammadi committed Mar 31, 2020 194 ``````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__. `````` Dennis Gläser committed Jul 19, 2018 195 `````` `````` Dennis Gläser committed Jul 19, 2018 196 ``````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. `````` Dennis Gläser committed Jul 19, 2018 197 `````` `````` Simon Scholz committed Jul 19, 2018 198 199 200 ``````## 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: `````` Dennis Gläser committed Jul 19, 2018 201 202 203 204 205 `````` ```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) `````` Theresa Schollenberger committed Dec 17, 2018 206 207 ``````static constexpr auto fractureDomainId = Dune::index_constant<1>(); const auto& fractureProblem = couplingManager().problem(fractureDomainId); `````` Dennis Gläser committed Jul 19, 2018 208 209 210 211 212 213 214 `````` // 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); ``` `````` Dennis Gläser committed Jul 19, 2018 215 `````` `````` Dennis Gläser committed Jul 19, 2018 216 ``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.``