Commit 0a212b09 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[examples] Add PNM upscaling example

parent 033877be
...@@ -5,3 +5,4 @@ add_subdirectory(shallowwaterfriction) ...@@ -5,3 +5,4 @@ add_subdirectory(shallowwaterfriction)
add_subdirectory(freeflowchannel) add_subdirectory(freeflowchannel)
add_subdirectory(1protationsymmetry) add_subdirectory(1protationsymmetry)
add_subdirectory(liddrivencavity) add_subdirectory(liddrivencavity)
add_subdirectory(porenetwork_upscaling)
...@@ -129,3 +129,21 @@ You learn how to ...@@ -129,3 +129,21 @@ You learn how to
<figure><img src="liddrivencavity/img/setup.png" alt="liddriven result"/></figure></td> <figure><img src="liddrivencavity/img/setup.png" alt="liddriven result"/></figure></td>
</a></td> </a></td>
</tr></table> </tr></table>
### [:open_file_folder: Example 7: Permeability estimation using a pore-network model](porenetwork_upscaling/README.md)
<table><tr><td>
In this example, we use a single-phase flow pore-network model to estimate the upscaled Darcy permeability of a randomly
generated grid.
You learn how to
* solve a single-phase-flow pore-network problem
* use the total mass flow rate to estimate $`K_{xx}`$, $`K_{yy}`$, $`K_{zz}`$
</td>
<td width="20%"><a href="liddrivencavity/README.md">
<figure><img src="porenetwork_upscaling/img/result.png" alt="pnm result"/></figure></td>
</a></td>
</tr></table>
{
"README.md" : [
"doc/_intro.md"
],
"doc/problem.md" : [
"doc/problem_intro.md",
"properties.hh",
"problem.hh"
],
"doc/main.md" : [
"doc/main_intro.md",
"main.cc"
],
"doc/upscalinghelper.md" : [
"doc/upscalinghelper_intro.md",
"upscalinghelper.hh"
],
"navigation" : {
"mainpage" : "README.md",
"subpages" : [
"doc/problem.md",
"doc/main.md",
"doc/upscalinghelper.md"
],
"subtitles" : [
"Permeability upscaling with a pore-network model",
"Main program flow",
"Upscaling"
]
}
}
dune_symlink_to_source_files(FILES "params.input")
dumux_add_test(NAME example_pnm1p_permeabilityupscaling
LABELS porenetwork example
SOURCES main.cc
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_CURRENT_BINARY_DIR}/example_pnm1p_permeabilityupscaling
CMD_ARGS -Problem.ReferenceData "3.083710e-13 2.870186e-13 2.671460e-13"
-Problem.TestEpsilon 1e-7)
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
# Determining the Darcy permeability of a pore network
__In this example, you will learn how to__
* simulate flow on a pore network by applying a pressure gradient in a given direction
* perform an upscaling in order to determine the intrinsic single-phase permeability $`\mathbf{K}`$ [m$`^2`$]
__Result__.
As a result of the simulation of this example, you will get the intrinsic single-phase permeabilities for each spatial direction $`K_{xx}`$, $`K_{yy}`$, $`K_{zz}`$ [m$`^2`$]
as direct output on your terminal.
Additionally, the output contains the values of the auxiliary parameters: cross-sectional area, mass flux integrated over the cross-sectional area,
and the resulting Darcy velocity as in the example below for the x-direction permeability $`K_{xx}`$:
```
x-direction: Area = 1.600000e-07 m^2; Massflux = 4.509338e-13 kg/s; v_Darcy = 2.818337e-09 m/s; K = 2.818337e-13 m^2
```
The pressure distribution throughout the pore network as well as pore-network characteristics will also be written to a vtp output file that can be viewed with ParaView.
Figure 1 shows the pressure distribution within the pore network for the case of flow in x-direction.
<figure>
<center>
<img src="img/result.png" alt="Pore-network pressure distribution" width="60%"/>
<figcaption> <b> Fig.1 </b> - Pressure distribution within the pore network for flow in x-direction. </figcaption>
</center>
</figure>
__Table of contents__. This description is structured as follows:
[[_TOC_]]
## Problem setup
We consider a single-phase problem within a randomly generated pore network of 20x20x20 pores cubed from which some of the pore throats were [deleted randomly](https://doi.org/10.1029/2010WR010180).
The inscribed pore body radii follow a truncated log-normal distribution.
To calculate the upscaled permeability, a pressure difference of 4000 Pa is applied sequentially in every direction, while all lateral sides are closed.
The resulting mass flow rate is then used to determine $`\mathbf{K}`$, as described [later](upscalinghelper.md).
## Mathematical and numerical model
In this example we are using the single-phase pore-network model of DuMu<sup>x</sup>, which considers a Hagen-Poiseuille-type law to relate the volume flow from on pore body to another to discrete pressure drops $`\Delta p = p_i - p_j`$ between the pore bodies.
We require mass conservation at each pore body $`i`$:
```math
\sum_j Q_{ij} = 0,
```
where $`Q_{ij}`$ is the discrete volume flow rate in a throat connecting pore bodies $`i`$ and $`j`$:
```math
Q_{ij} = g_{ij} (p_i - p_j).
```
# Implementation & Post processing
In the following, we take a closer look at the source files for this example.
We will discuss the different parts of the code in detail subsequently.
```
└── pnmpermeabilityupscaling/
├── CMakeLists.txt -> build system file
├── main.cc -> main program flow
├── params.input -> runtime parameters
├── properties.hh -> compile time settings for the simulation
├── problem.hh -> boundary & initial conditions for the simulation
├── spatialparams.hh -> parameter distributions for the simulation
└── upscalinghelper.hh -> helper for the upscaling calculations and writing the results
```
## Part 1: Permeability upscaling with a pore-network model
| [:arrow_right: Click to continue with part 1 of the documentation](doc/problem.md) |
|---:|
## Part 2: Main program flow
| [:arrow_right: Click to continue with part 2 of the documentation](doc/main.md) |
|---:|
## Part 3: Upscaling
| [:arrow_right: Click to continue with part 3 of the documentation](doc/upscalinghelper.md) |
|---:|
# Determining the Darcy permeability of a pore network
__In this example, you will learn how to__
* simulate flow on a pore network by applying a pressure gradient in a given direction
* perform an upscaling in order to determine the intrinsic single-phase permeability $`\mathbf{K}`$ [m$`^2`$]
__Result__.
As a result of the simulation of this example, you will get the intrinsic single-phase permeabilities for each spatial direction $`K_{xx}`$, $`K_{yy}`$, $`K_{zz}`$ [m$`^2`$]
as direct output on your terminal.
Additionally, the output contains the values of the auxiliary parameters: cross-sectional area, mass flux integrated over the cross-sectional area,
and the resulting Darcy velocity as in the example below for the x-direction permeability $`K_{xx}`$:
```
x-direction: Area = 1.600000e-07 m^2; Massflux = 4.509338e-13 kg/s; v_Darcy = 2.818337e-09 m/s; K = 2.818337e-13 m^2
```
The pressure distribution throughout the pore network as well as pore-network characteristics will also be written to a vtp output file that can be viewed with ParaView.
Figure 1 shows the pressure distribution within the pore network for the case of flow in x-direction.
<figure>
<center>
<img src="img/result.png" alt="Pore-network pressure distribution" width="60%"/>
<figcaption> <b> Fig.1 </b> - Pressure distribution within the pore network for flow in x-direction. </figcaption>
</center>
</figure>
__Table of contents__. This description is structured as follows:
[[_TOC_]]
## Problem setup
We consider a single-phase problem within a randomly generated pore network of 20x20x20 pores cubed from which some of the pore throats were [deleted randomly](https://doi.org/10.1029/2010WR010180).
The inscribed pore body radii follow a truncated log-normal distribution.
To calculate the upscaled permeability, a pressure difference of 4000 Pa is applied sequentially in every direction, while all lateral sides are closed.
The resulting mass flow rate is then used to determine $`\mathbf{K}`$, as described [later](upscalinghelper.md).
## Mathematical and numerical model
In this example we are using the single-phase pore-network model of DuMu<sup>x</sup>, which considers a Hagen-Poiseuille-type law to relate the volume flow from on pore body to another to discrete pressure drops $`\Delta p = p_i - p_j`$ between the pore bodies.
We require mass conservation at each pore body $`i`$:
```math
\sum_j Q_{ij} = 0,
```
where $`Q_{ij}`$ is the discrete volume flow rate in a throat connecting pore bodies $`i`$ and $`j`$:
```math
Q_{ij} = g_{ij} (p_i - p_j).
```
# Implementation & Post processing
In the following, we take a closer look at the source files for this example.
We will discuss the different parts of the code in detail subsequently.
```
└── pnmpermeabilityupscaling/
├── CMakeLists.txt -> build system file
├── main.cc -> main program flow
├── params.input -> runtime parameters
├── properties.hh -> compile time settings for the simulation
├── problem.hh -> boundary & initial conditions for the simulation
├── spatialparams.hh -> parameter distributions for the simulation
└── upscalinghelper.hh -> helper for the upscaling calculations and writing the results
```
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 1](problem.md) | [:arrow_right: Continue with part 3](upscalinghelper.md) |
|---|---|---:|
# Part 2: Main program flow
The main program flow is implemented in file `main.cc` described below.
For each spatial direction x, y and z, flow through the network is simulated and the resulting mass flow rate
is used to determine the permeability.
The code documentation is structured as follows:
[[_TOC_]]
## The main program (`main.cc`)
This file contains the main program flow. In this example, we use a single-phase
Pore-Network-Model to evaluate the upscaled Darcy permeability of a given network.
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>
### Includes
<details><summary> Click to show includes</summary>
```cpp
#include <config.h>
#include <iostream>
#include <dune/common/float_cmp.hh> // for floating point comparison
#include <dumux/common/properties.hh> // for GetPropType
#include <dumux/common/parameters.hh> // for getParam
#include <dumux/linear/seqsolverbackend.hh> // for ILU0BiCGSTABBackend
#include <dumux/linear/pdesolver.hh> // for LinearPDESolver
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
#include <dumux/io/grid/porenetwork/gridmanager.hh> // for pore-network grid
#include <dumux/porenetworkflow/common/pnmvtkoutputmodule.hh>
#include <dumux/porenetworkflow/common/boundaryflux.hh> // for getting the total mass flux leaving the network
#include "upscalinghelper.hh"
#include "properties.hh"
```
</details>
### Beginning of the main function
```cpp
int main(int argc, char** argv) try
{
using namespace Dumux;
// We parse the command line arguments.
Parameters::init(argc, argv);
// Convenience alias for the type tag of the problem.
using TypeTag = Properties::TTag::PNMUpscaling;
```
### Create the grid and the grid geometry
```cpp
// The grid manager can be used to create a grid from the input file
using GridManager = PoreNetwork::GridManager<3>;
GridManager gridManager;
gridManager.init();
// We compute on the leaf grid view.
const auto& leafGridView = gridManager.grid().leafGridView();
// instantiate the grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
auto gridData = gridManager.getGridData();
gridGeometry->update(*gridData);
```
### Initialize the problem and grid variables
```cpp
using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
auto spatialParams = std::make_shared<SpatialParams>(gridGeometry);
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry, spatialParams);
// instantiate and initialize the discrete and exact solution vectors
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(gridGeometry->numDofs());
// instantiate and initialize the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
```
### Initialize VTK output
```cpp
using VtkOutputFields = GetPropType<TypeTag, Properties::IOFields>;
using VtkWriter = PoreNetwork::VtkOutputModule<GridVariables, GetPropType<TypeTag, Properties::FluxVariables>, SolutionVector>;
VtkWriter vtkWriter(*gridVariables, x, problem->name());
VtkOutputFields::initOutputModule(vtkWriter);
vtkWriter.addField(gridGeometry->poreVolume(), "poreVolume", VtkWriter::FieldType::vertex);
vtkWriter.addField(gridGeometry->throatShapeFactor(), "throatShapeFactor", VtkWriter::FieldType::element);
vtkWriter.addField(gridGeometry->throatCrossSectionalArea(), "throatCrossSectionalArea", VtkWriter::FieldType::element);
```
### Instantiate the solver
We use the `LinearPDESolver` class, which is instantiated on the basis
of an assembler and a linear solver. When the `solve` function of the
`LinearPDESolver` is called, it uses the assembler and linear
solver classes to assemble and solve the linear system around the provided
solution and stores the result therein.
```cpp
using Assembler = FVAssembler<TypeTag, DiffMethod::analytic>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
using LinearSolver = UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
LinearPDESolver<Assembler, LinearSolver> solver(assembler, linearSolver);
solver.setVerbose(false); // suppress output during solve()
```
### Prepare the upscaling procedure.
Specify the directions for which the permeability shall be determined (default: x, y, z for 3D).
```cpp
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
const auto defaultDirections = GridView::dimensionworld == 3 ? std::vector<int>{0, 1, 2}
: std::vector<int>{0, 1};
const auto directions = getParam<std::vector<int>>("Problem.Directions", defaultDirections);
```
Set up a helper class to determine the total mass flux leaving the network
```cpp
const auto boundaryFlux = PoreNetwork::BoundaryFlux(*gridVariables, assembler->localResidual(), x);
```
Set the side lengths used for applying the pressure gradient and calculating the REV outflow area.
One can either specify these values manually (usually more accurate) or let the UpscalingHelper struct
determine it automatically based on the network's bounding box.
```cpp
const auto sideLengths = [&]()
{
using GlobalPosition = typename GridGeometry::GlobalCoordinate;
if (hasParam("Problem.SideLength"))
return getParam<GlobalPosition>("Problem.SideLength");
else
return UpscalingHelper::getSideLengths(*gridGeometry);
}();
// pass the side lengths to the problem
problem->setSideLengths(sideLengths);
```
### The actual upscaling procedure
Iterate over all directions specified before, apply the pressure gradient, calculated the mass flux
and finally determine the permeability.
```cpp
for (int dimIdx : directions)
{
// reset the solution
x = 0;
// set the direction in which the pressure gradient will be applied
problem->setDirection(dimIdx);
// solve problem
solver.solve(x);
// write the vtu file for the given direction
vtkWriter.write(dimIdx);
// get the Scalar type
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
// calculate the permeability
const Scalar totalFluidMassFlux = boundaryFlux.getFlux(std::vector<int>{problem->outletPoreLabel()})[0];
const Scalar K = UpscalingHelper::getDarcyPermeability(*problem, totalFluidMassFlux);
// optionally compare with reference data
static const auto referenceData = getParam<std::vector<Scalar>>("Problem.ReferenceData", std::vector<Scalar>{});
if (!referenceData.empty())
{
static const Scalar eps = getParam<Scalar>("Problem.TestEpsilon");
if (Dune::FloatCmp::ne<Scalar>(K, referenceData[dimIdx], eps))
{
std::cerr << "Calculated permeability of " << K << " does not match with reference value of " << referenceData[dimIdx] << std::endl;
return 1;
}
}
}
// program end, return with 0 exit code (success)
return 0;
}
```
### Exception handling
In this part of the main file we catch and print possible exceptions that could
occur during the simulation.
<details><summary> Click to show error handler</summary>
```cpp
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;
return 3;
}
```
</details>
</details>
| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 1](problem.md) | [:arrow_right: Continue with part 3](upscalinghelper.md) |
|---|---|---:|
# Part 2: Main program flow
The main program flow is implemented in file `main.cc` described below.
For each spatial direction x, y and z, flow through the network is simulated and the resulting mass flow rate
is used to determine the permeability.
The code documentation is structured as follows:
[[_TOC_]]
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_right: Continue with part 2](main.md) |
|---|---:|
# Part 1: Simulation setup
The code documentation is structured as follows:
[[_TOC_]]
## Compile-time settings (`properties.hh`)
This file defines the `TypeTag` used for the simulation in this example, for
which we specialize a number of compile-time `properties`.
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../properties.hh))</summary>
### Includes
<details><summary> Click to show includes</summary>
```cpp
#include <dune/grid/yaspgrid.hh> // for `Dune::YaspGrid`
```
The `OneP` type tag specializes most of the `properties` required for single-
phase flow simulations in DuMu<sup>x</sup>. We will use this in the following to inherit the
respective properties, and subsequently specialize those properties for our
type tag, which we want to modify or for which no meaningful default can be set.
```cpp
#include <dumux/porenetworkflow/1p/model.hh>// for `TTag::PNMOneP`
```
The local residual for incompressible flow is included.
The one-phase flow model (included above) uses a default implementation of the
local residual for single-phase flow. However, in this example we are using an
incompressible fluid phase. Therefore, we are including the specialized local
residual which contains functionality to analytically compute the entries of
the Jacobian matrix. We will use this in the main file.
```cpp
#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh>
```
We will use a single liquid phase consisting of a component with constant fluid properties.
```cpp
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
```
The classes that define the problem and parameters used in this simulation
```cpp
#include "problem.hh"
#include "spatialparams.hh"
```
</details>
### `TypeTag` definition
A `TypeTag` for our simulation is defined, which inherits properties from the
single-phase flow model and the box scheme.
```cpp
namespace Dumux::Properties {
namespace TTag {
struct PNMUpscaling { using InheritsFrom = std::tuple<PNMOneP>; };
}
```
### Property specializations
In the following piece of code, mandatory `properties` for which no meaningful
default can be set, are specialized for our type tag `PNMUpscaling`.
```cpp
// We use `dune-foamgrid`, which is especially tailored for 1D networks.
template<class TypeTag>
struct Grid<TypeTag, TTag::PNMUpscaling>
{ using type = Dune::FoamGrid<1, 3>; };
// The problem class specifying initial and boundary conditions:
template<class TypeTag>
struct Problem<TypeTag, TTag::PNMUpscaling>
{ using type = UpscalingProblem<TypeTag>; };
//! The spatial parameters to be employed.
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::PNMUpscaling>
{
private:
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = Dumux::PoreNetwork::UpscalingSpatialParams<GridGeometry, Scalar>;
};
//! The advection type.
template<class TypeTag>
struct AdvectionType<TypeTag, TTag::PNMUpscaling>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using TransmissibilityLaw = Dumux::PoreNetwork::TransmissibilityPatzekSilin<Scalar, true/*considerPoreBodyResistance*/>;
public:
using type = Dumux::PoreNetwork::CreepingFlow<Scalar, TransmissibilityLaw>;
};
// We use a single liquid phase consisting of a component with constant fluid properties.
template<class TypeTag>