main.md 16.7 KB
Newer Older
1
2
3
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->


4
| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 2](tracer.md) |
5
6
|---|---:|

7
# Part 3: Main program flow
8

9
10
11
12
13
14
15
16
We want to solve a single-phase flow problem
to obtain a pressure distribution in the domain. Subsequently, we compute a volume fluxes field
based on the obtained pressure distribution and pass it to the tracer problem.
Finally, we solve a transient transport problem for a tracer using the computed volume fluxes.
This main program flow is implemented in the `main()` function
of the program which is defined in the file `main.cc` described below.

The code documentation is structured as follows:
17
18
19
20
21

[[_TOC_]]


## The main file (`main.cc`)
22
23
24
25
26

<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>


27
### Included header files
28
<details><summary> Click to show includes</summary>
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
These are DUNE helper classes related to parallel computations, time measurements and file I/O

```cpp
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
```

The following headers include functionality related to property definition or retrieval, as well as
the retrieval of input parameters specified in the input file or via the command line.

```cpp
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
```

The following files contains the available linear solver backends and the assembler for the linear
systems arising from finite volume discretizations (box-scheme, tpfa-approximation, mpfa-approximation).

```cpp
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh> // analytic or numeric differentiation
```

The following class provides a convenient way of writing of dumux simulation results to VTK format.

```cpp
#include <dumux/io/vtkoutputmodule.hh>
```

The gridmanager constructs a grid from the information in the input or grid file.
Many different Dune grid implementations are supported, of which a list can be found
in `gridmanager.hh`.

```cpp
#include <dumux/io/grid/gridmanager.hh>
```

For both the single-phase and the tracer problem, `TypeTags` are defined, which collect
the properties that are required for the simulation. These type tags are defined
in the headers that we include here. For detailed information, please have a look
at the documentation provided therein.

```cpp
#include "properties_1p.hh"
#include "properties_tracer.hh"
```

</details>

79
### The main function
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
We will now discuss the main program flow implemented within the `main` function.
At the beginning of each program using Dune, an instance of `Dune::MPIHelper` has to
be created. Moreover, we parse the run-time arguments from the command line and the
input file:

```cpp
int main(int argc, char** argv) try
{
    using namespace Dumux;

    // The Dune MPIHelper must be instantiated for each program using Dune
    Dune::MPIHelper::instance(argc, argv);

    // parse command line arguments and input file
    Parameters::init(argc, argv);
```

We define convenience aliases for the type tags of the two problems. The type
tags contain all the properties that are needed to define the model and the problem
setup. Throughout the main file, we will obtain types defined for these type tags
using the property system, i.e. with `GetPropType`. A more detailed documentation
for the type tags and properties can be found in the documentation of the single-phase
and the tracer transport setups, respectively.

```cpp
    using OnePTypeTag = Properties::TTag::IncompressibleTest;
    using TracerTypeTag = Properties::TTag::TracerTest;
```

109
#### Step 1: Create the grid
110
111
112
113
114
115
116
117
118
119
120
121
122
123
The `GridManager` class creates the grid from information given in the input file.
This can either be a grid file, or in the case of structured grids, one can specify the coordinates
of the corners of the grid and the number of cells to be used to discretize each spatial direction.
Here, we solve both the single-phase and the tracer problem on the same grid, and thus,
the grid is only created once using the grid type defined by the `OnePTypeTag` of the single-phase problem.

```cpp
    GridManager<GetPropType<OnePTypeTag, Properties::Grid>> gridManager;
    gridManager.init();

    // We compute on the leaf grid view.
    const auto& leafGridView = gridManager.grid().leafGridView();
```

124
#### Step 2: Solving the single-phase problem
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
First, a finite volume grid geometry is constructed from the grid that was created above.
This builds the sub-control volumes (scv) and sub-control volume faces (scvf) for each element
of the grid partition.

```cpp
    using GridGeometry = GetPropType<OnePTypeTag, Properties::GridGeometry>;
    auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
    gridGeometry->update();
```

We now instantiate the problem, in which we define the boundary and initial conditions.

```cpp
    using OnePProblem = GetPropType<OnePTypeTag, Properties::Problem>;
    auto problemOneP = std::make_shared<OnePProblem>(gridGeometry);
```

The jacobian matrix (`A`), the solution vector (`p`) and the residual (`r`) make up the linear system.

```cpp
    using JacobianMatrix = GetPropType<OnePTypeTag, Properties::JacobianMatrix>;
    using SolutionVector = GetPropType<OnePTypeTag, Properties::SolutionVector>;
    SolutionVector p(leafGridView.size(0));

    auto A = std::make_shared<JacobianMatrix>();
    auto r = std::make_shared<SolutionVector>();
```

The grid variables are used store variables (primary and secondary variables) on sub-control volumes and faces (volume and flux variables).

```cpp
    using OnePGridVariables = GetPropType<OnePTypeTag, Properties::GridVariables>;
    auto onePGridVariables = std::make_shared<OnePGridVariables>(problemOneP, gridGeometry);
    onePGridVariables->init(p);
```

We now instantiate the assembler class, assemble the linear system and solve it with the linear
solver UMFPack. Besides that, the time needed for assembly and solve is measured and printed.

```cpp
    using OnePAssembler = FVAssembler<OnePTypeTag, DiffMethod::analytic>;
    auto assemblerOneP = std::make_shared<OnePAssembler>(problemOneP, gridGeometry, onePGridVariables);
    assemblerOneP->setLinearSystem(A, r); // tell assembler to use our previously defined system

    Dune::Timer timer;
    Dune::Timer assemblyTimer; std::cout << "Assembling linear system ..." << std::flush;
    assemblerOneP->assembleJacobianAndResidual(p); // assemble linear system around current solution
    assemblyTimer.stop(); std::cout << " took " << assemblyTimer.elapsed() << " seconds." << std::endl;

    (*r) *= -1.0; // We want to solve `Ax = -r`.

    using LinearSolver = UMFPackBackend;
    Dune::Timer solverTimer; std::cout << "Solving linear system ..." << std::flush;
    auto linearSolver = std::make_shared<LinearSolver>();
    linearSolver->solve(*A, p, *r);
    solverTimer.stop(); std::cout << " took " << solverTimer.elapsed() << " seconds." << std::endl;

    Dune::Timer updateTimer; std::cout << "Updating variables ..." << std::flush;
    onePGridVariables->update(p); // update grid variables to new pressure distribution
    updateTimer.elapsed(); std::cout << " took " << updateTimer.elapsed() << std::endl;
```

The solution vector `p` now contains the pressure field, i.e.the solution to the single-phase
problem defined in `problem_1p.hh`. Let us now write this solution to a VTK file using the Dune
`VTKWriter`. Moreover, we add the permeability distribution to the writer.

```cpp
192
    using GridView = typename GridGeometry::GridView;
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
    Dune::VTKWriter<GridView> onepWriter(leafGridView);
    onepWriter.addCellData(p, "p");

    const auto& k = problemOneP->spatialParams().getKField(); // defined in spatialparams_1p.hh
    onepWriter.addCellData(k, "permeability"); // add permeability to writer
    onepWriter.write("1p");  // write the file "1p.vtk"

    // print overall CPU time required for assembling and solving the 1p problem.
    timer.stop();
    const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
    std::cout << "Simulation took " << timer.elapsed() << " seconds on "
              << comm.size() << " processes.\n"
              << "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
```

208
#### Step 3: Computation of the volume fluxes
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
248
249
250
251
We use the results of the 1p problem to calculate the volume fluxes across all sub-control volume
faces of the discretization and store them in the vector `volumeFlux`. In order to do so, we iterate
over all elements of the grid, and in each element compute the volume fluxes for all sub-control volume
faces embeded in that element.

```cpp
    using Scalar =  GetPropType<OnePTypeTag, Properties::Scalar>; // type for scalar values
    std::vector<Scalar> volumeFlux(gridGeometry->numScvf(), 0.0);

    using FluxVariables =  GetPropType<OnePTypeTag, Properties::FluxVariables>;
    auto upwindTerm = [](const auto& volVars) { return volVars.mobility(0); };

    // We iterate over all elements
    for (const auto& element : elements(leafGridView))
    {
        // Compute the element-local views on geometry, primary and secondary variables
        // as well as variables needed for flux computations.

        // This creates instances of the local views
        auto fvGeometry = localView(*gridGeometry);
        auto elemVolVars = localView(onePGridVariables->curGridVolVars());
        auto elemFluxVars = localView(onePGridVariables->gridFluxVarsCache());

        // we now have to bind the views to the current element
        fvGeometry.bind(element);
        elemVolVars.bind(element, fvGeometry, p);
        elemFluxVars.bind(element, fvGeometry, elemVolVars);

        // We calculate the volume fluxes for all sub-control volume faces except for Neumann boundary faces
        for (const auto& scvf : scvfs(fvGeometry))
        {
            // skip Neumann boundary faces
            if (scvf.boundary() && problemOneP->boundaryTypes(element, scvf).hasNeumann())
                continue;

            // let the FluxVariables class do the flux computation.
            FluxVariables fluxVars;
            fluxVars.init(*problemOneP, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
            volumeFlux[scvf.index()] = fluxVars.advectiveFlux(0, upwindTerm);
        }
    }
```

252
#### Step 4: Solving the tracer transport problem
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
First, we instantiate the tracer problem containing initial and boundary conditions,
and pass to it the previously computed volume fluxes (see the documentation of the
file `spatialparams_tracer.hh` for more details).

```cpp
    using TracerProblem = GetPropType<TracerTypeTag, Properties::Problem>;
    auto tracerProblem = std::make_shared<TracerProblem>(gridGeometry);
    tracerProblem->spatialParams().setVolumeFlux(volumeFlux);
```

We create and initialize the solution vector. In contrast to the steady-state single-phase problem,
the tracer problem is transient, and the initial solution defined in the problem is applied to the
solution vector. On the basis of this solution, we initialize then the grid variables.

```cpp
    SolutionVector x(leafGridView.size(0));
    tracerProblem->applyInitialSolution(x);
    auto xOld = x;

    using GridVariables = GetPropType<TracerTypeTag, Properties::GridVariables>;
    auto gridVariables = std::make_shared<GridVariables>(tracerProblem, gridGeometry);
    gridVariables->init(x);
```

Let us now instantiate the time loop. Therefore, we read in some time loop parameters from the input file.
The parameter `tEnd` defines the duration of the simulation, `dt` the initial time step size and `maxDt` the maximal time step size.
Moreover, we define 10 check points in the time loop at which we will write the solution to vtk files.

```cpp
    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
    auto dt = getParam<Scalar>("TimeLoop.DtInitial");
    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");

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

We create and initialize the assembler with a time loop for the transient problem.
Within the time loop, we will use this assembler in each time step to assemble the linear system.

```cpp
    using TracerAssembler = FVAssembler<TracerTypeTag, DiffMethod::analytic, /*implicit=*/false>;
297
    auto assembler = std::make_shared<TracerAssembler>(tracerProblem, gridGeometry, gridVariables, timeLoop, xOld);
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
    assembler->setLinearSystem(A, r);
```

The following lines of code initialize the vtk output module, add the velocity output facility
and write out the initial solution. At each checkpoint, we will use the output module to write
the solution of a time step into a corresponding vtk file.

```cpp
    VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, tracerProblem->name());

    // add model-specific output fields to the writer
    using IOFields = GetPropType<TracerTypeTag, Properties::IOFields>;
    IOFields::initOutputModule(vtkWriter);

    // add velocity output facility
    using VelocityOutput = GetPropType<TracerTypeTag, Properties::VelocityOutput>;
    vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));

    // write initial solution
    vtkWriter.write(0.0);
```

320
##### The time loop
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
We start the time loop and solve a new time step as long as `tEnd` is not reached. In every time step,
the problem is assembled and solved, the solution is updated, and when a checkpoint is reached the solution
is written to a new vtk file. In addition, statistics related to CPU time, the current simulation time
and the time step sizes used is printed to the terminal.

```cpp
    timeLoop->start(); do
    {
        // Then the linear system is assembled.
        Dune::Timer assembleTimer;
        assembler->assembleJacobianAndResidual(x);
        assembleTimer.stop();

        // We solve the linear system `A(xOld-xNew) = r`.
        Dune::Timer solveTimer;
        SolutionVector xDelta(x);
        linearSolver->solve(*A, xDelta, *r);
        solveTimer.stop();

        // update the solution vector and the grid variables.
        updateTimer.reset();
        x -= xDelta;
        gridVariables->update(x);
        updateTimer.stop();

        // display the statistics of the actual time step.
        const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
        std::cout << "Assemble/solve/update time: "
                  <<  assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/"
                  <<  solveTimer.elapsed() << "(" << 100*solveTimer.elapsed()/elapsedTot << "%)/"
                  <<  updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
                  <<  std::endl;

        // Update the old solution with the one computed in this time step and move to the next one
        xOld = x;
        gridVariables->advanceTimeStep();
        timeLoop->advanceTimeStep();

        // Write the Vtk output on check points.
        if (timeLoop->isCheckPoint())
            vtkWriter.write(timeLoop->time());

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

        // set the time step size for the next time step
        timeLoop->setTimeStepSize(dt);

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

The following piece of code prints a final status report of the time loop
before the program is terminated.

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

    return 0;
}
```

382
#### Exception handling
383
384
In this part of the main file we catch and print possible exceptions that could
occur during the simulation.
385
<details><summary> Click to show exception handler</summary>
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419

```cpp
// errors related to run-time parameters
catch (Dumux::ParameterException &e)
{
    std::cerr << std::endl << e << " ---> Abort!" << std::endl;
    return 1;
}
// errors related to the parsing of Dune grid files
catch (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;
}
// generic error handling with Dune::Exception
catch (Dune::Exception &e)
{
    std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
    return 3;
}
// other exceptions
catch (...)
{
    std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
    return 4;
}
```

</details>

420
421
</details>

422

423
| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 2](tracer.md) |
424
425
|---|---:|