| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | |||
| 8 | #ifndef DUMUX_INCOMPRESSIBLE_ONEP_TEST_SPATIAL_PARAMS_HH | ||
| 9 | #define DUMUX_INCOMPRESSIBLE_ONEP_TEST_SPATIAL_PARAMS_HH | ||
| 10 | |||
| 11 | // ## Parameter distributions (`spatialparams_1p.hh`) | ||
| 12 | // | ||
| 13 | // This file contains the __spatial parameters class__ which defines the | ||
| 14 | // distributions for the porous medium parameters permeability and porosity | ||
| 15 | // over the computational grid. | ||
| 16 | // | ||
| 17 | // [[content]] | ||
| 18 | // | ||
| 19 | // ### Include files | ||
| 20 | // In this example, we use a randomly generated and element-wise distributed | ||
| 21 | // permeability field. For this, we use the random number generation facilities | ||
| 22 | // provided by the C++ standard library. | ||
| 23 | // DuMu<sup>x</sup> additionally implements some simpler but | ||
| 24 | // standard-library-implementation-independent random distributions | ||
| 25 | // that are accurate enough for our purposes | ||
| 26 | // but allow to generate reproducible and portable random fields (when using a constant seed) | ||
| 27 | // for testing purposes. | ||
| 28 | #include <random> | ||
| 29 | #include <dumux/common/random.hh> | ||
| 30 | |||
| 31 | // We include the spatial parameters class for single-phase models discretized | ||
| 32 | // by finite volume schemes, from which the spatial parameters defined for this | ||
| 33 | // example will inherit. | ||
| 34 | #include <dumux/porousmediumflow/fvspatialparams1p.hh> | ||
| 35 | // | ||
| 36 | // ### The spatial parameters class | ||
| 37 | // | ||
| 38 | // In the `OnePTestSpatialParams` class, we define all functions needed to describe | ||
| 39 | // the porous medium, e.g. porosity and permeability, for the simulated single-phase flow scenario. | ||
| 40 | // We inherit from the `FVPorousMediumFlowSpatialParamsOneP` class here, which is the base class | ||
| 41 | // for spatial parameters in the context of single-phase porous medium flow | ||
| 42 | // applications using finite volume discretization schemes. | ||
| 43 | // [[codeblock]] | ||
| 44 | namespace Dumux { | ||
| 45 | |||
| 46 | template<class GridGeometry, class Scalar> | ||
| 47 | class OnePTestSpatialParams | ||
| 48 | : public FVPorousMediumFlowSpatialParamsOneP<GridGeometry, Scalar, | ||
| 49 | OnePTestSpatialParams<GridGeometry, Scalar>> | ||
| 50 | { | ||
| 51 | // The following convenience aliases will be used throughout this class. | ||
| 52 | using GridView = typename GridGeometry::GridView; | ||
| 53 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 54 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 55 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 56 | using ParentType = FVPorousMediumFlowSpatialParamsOneP<GridGeometry, Scalar, | ||
| 57 | OnePTestSpatialParams<GridGeometry, Scalar>>; | ||
| 58 | |||
| 59 | static constexpr int dimWorld = GridView::dimensionworld; | ||
| 60 | using GlobalPosition = typename SubControlVolume::GlobalPosition; | ||
| 61 | |||
| 62 | public: | ||
| 63 | // The spatial parameters must export the type used to define permeabilities. | ||
| 64 | // Here, we are using scalar permeabilities, but tensors are also supported. | ||
| 65 | using PermeabilityType = Scalar; | ||
| 66 | // [[/codeblock]] | ||
| 67 | // | ||
| 68 | // #### Generation of the random permeability field | ||
| 69 | // We generate a random permeability field upon construction of the spatial parameters class | ||
| 70 | // using lognormal distributions. The mean values for the permeability inside and outside of a | ||
| 71 | // low-permeable lens (given by the coordinates `lensLowerLeft_` and `lensUpperRight_`) are defined | ||
| 72 | // in the variables `permeabilityLens_` and `permeability_`. The respective values are obtained | ||
| 73 | // from the input file (`params.input`) making use of the free function `getParam`. We use a standard deviarion | ||
| 74 | // of 10% here and compute permeabily values for all elements of the computational grid. | ||
| 75 | // [[codeblock]] | ||
| 76 | 1 | OnePTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry) | |
| 77 |
10/26✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 1 times.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
|
7 | : ParentType(gridGeometry), K_(gridGeometry->gridView().size(0), 0.0) |
| 78 | { | ||
| 79 | // The permeability of the domain and the lens are obtained from the `params.input` file. | ||
| 80 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | permeability_ = getParam<Scalar>("SpatialParams.Permeability"); |
| 81 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | permeabilityLens_ = getParam<Scalar>("SpatialParams.PermeabilityLens"); |
| 82 | |||
| 83 | // Furthermore, the position of the lens, which is defined by the position of the lower left and the upper right corners, are obtained from the input file. | ||
| 84 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft"); |
| 85 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight"); |
| 86 | |||
| 87 | // We generate random fields for the permeability using lognormal distributions, | ||
| 88 | // with `permeability_` as mean value and 10 % of it as standard deviation. | ||
| 89 | // A separate distribution is used for the lens using `permeabilityLens_`. | ||
| 90 | // A permeability value is created for each element of the grid and is stored in the vector `K_`. | ||
| 91 | // For a "truly" random field (different every time we execute) the seed `0` for the | ||
| 92 | // Mersenne-Twister pseudo-random-number generator should be replaced | ||
| 93 | // by `std::random_device{}()`. For unbiased high-quality distributions (usually not needed) | ||
| 94 | // replace `Dumux::SimpleLogNormalDistribution` by `std::lognormal_distribution`. | ||
| 95 | 1 | std::mt19937 rand(0); | |
| 96 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | Dumux::SimpleLogNormalDistribution<Scalar> K(std::log(permeability_), std::log(permeability_)*0.1); |
| 97 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | Dumux::SimpleLogNormalDistribution<Scalar> KLens(std::log(permeabilityLens_), std::log(permeabilityLens_)*0.1); |
| 98 | |||
| 99 | // loop over all elements and compute a permeability value | ||
| 100 |
4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2501 times.
✗ Branch 10 not taken.
|
5003 | for (const auto& element : elements(gridGeometry->gridView())) |
| 101 | { | ||
| 102 | 7500 | const auto eIdx = gridGeometry->elementMapper().index(element); | |
| 103 | 2500 | const auto globalPos = element.geometry().center(); | |
| 104 |
2/2✓ Branch 0 taken 900 times.
✓ Branch 1 taken 1600 times.
|
5000 | K_[eIdx] = isInLens_(globalPos) ? KLens(rand) : K(rand); |
| 105 | } | ||
| 106 | 1 | } | |
| 107 | // [[/codeblock]] | ||
| 108 | |||
| 109 | // #### Properties of the porous matrix | ||
| 110 | // This function returns the permeability $`[m^2]`$ to be used within a sub-control volume (`scv`) inside the element `element`. | ||
| 111 | // One can define the permeability as function of the primary variables on the element, which are given in the provided | ||
| 112 | // instance of `ElementSolution`. Here, we use element-wise distributed permeabilities that were randomly generated in | ||
| 113 | // the constructor and stored in the vector `K_`(see above). | ||
| 114 | template<class ElementSolution> | ||
| 115 | ✗ | const PermeabilityType& permeability(const Element& element, | |
| 116 | const SubControlVolume& scv, | ||
| 117 | const ElementSolution& elemSol) const | ||
| 118 | { | ||
| 119 | ✗ | return K_[scv.elementIndex()]; | |
| 120 | } | ||
| 121 | |||
| 122 | // We set the porosity $`[-]`$ for the whole domain to a value of $`20 \%`$. | ||
| 123 | // Note that in case you want to use solution-dependent porosities, you can | ||
| 124 | // use the overload | ||
| 125 | // `porosity(const Element&, const SubControlVolume&, const ElementSolution&)` | ||
| 126 | // that is defined in the base class `FVPorousMediumFlowSpatialParamsOneP`. Per default, this | ||
| 127 | // forwards to the `porosityAtPos` function per default, which we overload here. | ||
| 128 | ✗ | Scalar porosityAtPos(const GlobalPosition& globalPos) const | |
| 129 | ✗ | { return 0.2; } | |
| 130 | |||
| 131 | // We reference to the permeability field. This is used in the main function to write an output for the permeability field. | ||
| 132 | const std::vector<Scalar>& getKField() const | ||
| 133 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | { return K_; } |
| 134 | |||
| 135 | // The remainder of this class contains a convenient function to determine if | ||
| 136 | // a position is inside the lens and defines the data members. | ||
| 137 | // [[details]] private data members and member functions | ||
| 138 | // [[codeblock]] | ||
| 139 | private: | ||
| 140 | // The following function returns true if a given position is inside the lens. | ||
| 141 | // We use an epsilon of 1.5e-7 here for floating point comparisons. | ||
| 142 | bool isInLens_(const GlobalPosition& globalPos) const | ||
| 143 | { | ||
| 144 |
2/2✓ Branch 0 taken 4000 times.
✓ Branch 1 taken 900 times.
|
4900 | for (int i = 0; i < dimWorld; ++i) |
| 145 | { | ||
| 146 |
4/4✓ Branch 0 taken 3200 times.
✓ Branch 1 taken 800 times.
✓ Branch 2 taken 3200 times.
✓ Branch 3 taken 800 times.
|
8000 | if (globalPos[i] < lensLowerLeft_[i] + 1.5e-7 |
| 147 |
8/8✓ Branch 0 taken 3200 times.
✓ Branch 1 taken 800 times.
✓ Branch 2 taken 2400 times.
✓ Branch 3 taken 800 times.
✓ Branch 4 taken 2400 times.
✓ Branch 5 taken 800 times.
✓ Branch 6 taken 2400 times.
✓ Branch 7 taken 800 times.
|
4000 | || globalPos[i] > lensUpperRight_[i] - 1.5e-7) |
| 148 | return false; | ||
| 149 | } | ||
| 150 | |||
| 151 | return true; | ||
| 152 | } | ||
| 153 | |||
| 154 | GlobalPosition lensLowerLeft_, lensUpperRight_; | ||
| 155 | Scalar permeability_, permeabilityLens_; | ||
| 156 | std::vector<Scalar> K_; | ||
| 157 | }; // end class definition of OnePTestSpatialParams | ||
| 158 | } // end namespace Dumux | ||
| 159 | // [[/codeblock]] | ||
| 160 | // [[/details]] | ||
| 161 | // [[/content]] | ||
| 162 | #endif | ||
| 163 |