main.md 8.74 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
<!-- 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) |
|---|---|---:|