GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/multidomain/embedded/1d3d/root_soil_benchmark/couplingreconstruction.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 84 89 94.4%
Functions: 10 12 83.3%
Branches: 68 120 56.7%
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 * \file
9 * \ingroup EmbeddedCoupling
10 * \brief Coupling reconstruction for 1d-3d mixed-dimension method with distributed sources
11 */
12 #ifndef DUMUX_MULTIDOMAIN_EMBEDDED_COUPLING_RECONSTRUCTION_1D3D_KERNEL_ROOTSOIL_HH
13 #define DUMUX_MULTIDOMAIN_EMBEDDED_COUPLING_RECONSTRUCTION_1D3D_KERNEL_ROOTSOIL_HH
14
15 #include <cmath>
16 #include <algorithm>
17 #include <memory>
18 #include <vector>
19
20 #include <dune/common/exceptions.hh>
21
22 #include <dumux/common/integrate.hh>
23 #include <dumux/common/parameters.hh>
24 #include <dumux/common/math.hh>
25
26 #include <dumux/nonlinear/findscalarroot.hh>
27
28 namespace Dumux::Detail::RootSoil {
29
30 /*!
31 * \file
32 * \ingroup EmbeddedCoupling
33 * \brief Kirchhoff transformation for the Richards equation
34 */
35 class KirchhoffTransformation
36 {
37 double alpha, n, m, l;
38 public:
39 2 KirchhoffTransformation()
40 2 {
41 2 alpha = getParam<double>("Soil.SpatialParams.VanGenuchtenAlpha");
42 2 n = getParam<double>("Soil.SpatialParams.VanGenuchtenN");
43 2 l = getParam<double>("Soil.SpatialParams.VanGenuchtenL", 0.5);
44 2 m = 1.0 - 1.0/n;
45 2 }
46
47 76945570 inline double kr(const double pw) const
48 {
49 76945570 using std::pow;
50 76945570 const auto pc = 1.0e5-pw;
51 76945570 const auto swe = pow(pow(alpha*pc, n) + 1.0, -m);
52 76945570 const auto r = 1.0 - pow(1.0 - pow(swe, 1.0/m), m);
53 76945570 return pow(swe, l)*r*r;
54 }
55
56 //! Kirchhoff transform
57 inline double computeTransformed(const double p) const
58 38572786 { return integrateScalarFunction([this](const auto& p){ return kr(p); }, 0.0, p, 1e-15); }
59
60 //! Inverse Kirchhoff transform
61 inline double computeInverseTransformed(const double psi) const
62 {
63 const auto residual = [this, psi](const auto& p){ return psi - computeTransformed(p); };
64 return findScalarRootBrent(-2.6e7, 1.0e5, residual, 1e-15, 200000);
65 }
66 };
67
68 /*!
69 * \file
70 * \ingroup EmbeddedCoupling
71 * \brief Kirchhoff transformation for the Richards equation with values cached in a lookup table
72 */
73 class KirchhoffTransformationCached
74 {
75 double pMin, pMax;
76 std::size_t samples;
77 KirchhoffTransformation transformation_;
78 std::vector<double> p_, psi_;
79 public:
80 2 KirchhoffTransformationCached()
81 2 {
82
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 pMin = getParam<double>("MixedDimension.Reconstruction.PressureMin", -2.6e7);
83
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 pMax = getParam<double>("MixedDimension.Reconstruction.PressureMax", 1.0e5);
84
1/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 2 samples = getParam<std::size_t>("MixedDimension.Reconstruction.Samples", 10000);
85
86
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 p_.resize(samples);
87
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 psi_.resize(samples);
88 2 const double pStep = (pMax-pMin)/samples;
89
2/2
✓ Branch 0 taken 200002 times.
✓ Branch 1 taken 2 times.
 200004 for (std::size_t i = 0; i <= samples; ++i)
90 {
91 200002 const auto pValue = pMin + i*pStep;
92 200002 p_[i] = pValue;
93 200002 psi_[i] = transformation_.computeTransformed(pValue);
94 }
95
96 2 std::cout << "[coupling] Presampled Kirchhoff transformation with "
97
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 4 << samples << " samples (step size: " << pStep << " Pa)."
98
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << std::endl;
99 2 }
100
101 //! Kirchhoff transform
102 inline double computeTransformed(const double p) const
103 198680122 { return interpolate<InterpolationPolicy::LinearTable>(p, p_, psi_); }
104
105 //! Inverse Kirchhoff transform
106 inline double computeInverseTransformed(const double psi) const
107 28 { return interpolate<InterpolationPolicy::LinearTable>(psi, psi_, p_); }
108 };
109
110 /*!
111 * \file
112 * \ingroup EmbeddedCoupling
113 * \brief Caching the source computation to speed up Jacobian approximation during assembly
114 */
115 struct SourceCache
116 {
117 SourceCache()
118 2658 : q_{}
119 {
120 2658 p0_.fill(1e20);
121 2658 pRoot_.fill(1e20);
122 }
123
124 template<class F>
125 298723402 double operator() (const double p0, const double pRoot, const F& f)
126 {
127
2/2
✓ Branch 0 taken 515558308 times.
✓ Branch 1 taken 4716678 times.
 520274986 for (std::size_t i = 0; i < size; ++i)
128
12/12
✓ Branch 0 taken 417376205 times.
✓ Branch 1 taken 98182103 times.
✓ Branch 2 taken 417376205 times.
✓ Branch 3 taken 98182103 times.
✓ Branch 4 taken 417376205 times.
✓ Branch 5 taken 98182103 times.
✓ Branch 6 taken 294006724 times.
✓ Branch 7 taken 123369481 times.
✓ Branch 8 taken 294006724 times.
✓ Branch 9 taken 123369481 times.
✓ Branch 10 taken 294006724 times.
✓ Branch 11 taken 123369481 times.
 1546674924 if (std::abs(p0-p0_[i]) < 1e-14 && std::abs(pRoot-pRoot_[i]) < 1e-14)
129 588013448 return q_[i];
130
131 4716678 const auto q = f(p0, pRoot);
132
2/2
✓ Branch 0 taken 3145743 times.
✓ Branch 1 taken 1570934 times.
 4716677 p0_[writePos] = p0;
133
2/2
✓ Branch 0 taken 3145743 times.
✓ Branch 1 taken 1570934 times.
 4716677 pRoot_[writePos] = pRoot;
134
2/2
✓ Branch 0 taken 3145743 times.
✓ Branch 1 taken 1570934 times.
 4716677 q_[writePos] = q;
135
2/2
✓ Branch 0 taken 3145743 times.
✓ Branch 1 taken 1570934 times.
 4716677 writePos = (writePos >= size-1) ? 0 : writePos+1;
136 4716677 return q;
137 }
138
139 private:
140 static constexpr std::size_t size = 3; // original and each variable deflected
141 std::array<double, size> q_, p0_, pRoot_;
142 std::size_t writePos = 0;
143 };
144
145 } // end namespace Dumux::Detail::RootSoil
146
147 namespace Dumux::RootSoil {
148
149 /*!
150 * \file
151 * \ingroup EmbeddedCoupling
152 * \brief Coupling reconstruction as derived in Koch et al (2022), https://doi.org/10.1016/j.jcp.2021.110823
153 */
154 template<class PcKrSwCurve>
155 class CouplingReconstruction
156 {
157 using CacheStorage = std::vector<std::vector<std::pair<std::size_t, Detail::RootSoil::SourceCache>>>;
158 public:
159 2 CouplingReconstruction()
160 4 {
161
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 K_ = getParam<double>("Soil.SpatialParams.Permeability");
162
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 porosity_ = getParam<double>("Soil.SpatialParams.Porosity");
163 // relative tolerance in term of the unknown approximation (psBar)
164
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 nonlinearSolverTolerance_ = getParam<double>("MixedDimension.Reconstruction.RelativeTolerance", 1e-4);
165
166
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 4 std::cout << "\n" << "[coupling] Test Kirchhoff transform:\n";
167 const auto& tf = transformation_;
168
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 2 times.
 16 for (const auto p : {0.9e5, 0.5e5, 0.1e5, 0.0, -1.0e5, -5.0e5, -1.5e6})
169 {
170 14 const auto psi = tf.computeTransformed(p);
171 14 const auto pNum = tf.computeInverseTransformed(psi);
172
4/8
✓ Branch 2 taken 14 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 14 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 14 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 14 times.
✗ Branch 14 not taken.
 56 std::cout << "pw: " << p << ", psi: " << psi << ", pw (inverse): " << pNum << '\n';
173
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 14 times.
 28 if (std::abs(p-pNum) > 100)
174 ✗ std::cerr << "Transformation is not precise enough. Maybe decrease the tolerances!" << std::endl;
175 }
176
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 std::cout << std::endl;
177 2 }
178
179 4716678 double computePsBar(const double p0, const double pRoot,
180 const double kernelRadius, const double rootRadius,
181 const double Kr, const double density, const double viscosity,
182 const std::size_t sourceId,
183 const double delta) const
184 {
185 4716678 const auto& tf = transformation_;
186 4716678 const auto psi0 = tf.computeTransformed(p0);
187 4716678 const auto sourceFactor = (std::log(kernelRadius/rootRadius) - 0.5 + 0.5*delta*delta/kernelRadius/kernelRadius)*viscosity/(2.0*M_PI*K_);
188
189 572456892 const auto residual = [&](const double psBar){
190 378493476 const auto sourceTerm = sourceFactor*2.0*M_PI*rootRadius*Kr*(psBar - pRoot);
191 94623369 return psi0 - tf.computeTransformed(psBar) - sourceTerm;
192 };
193
194 4716678 return findScalarRootBrent(-2.5e7, 1.0e5, residual, nonlinearSolverTolerance_, 20000);
195 }
196
197 double computeSource(const double p0, const double pRoot,
198 const double kernelRadius, const double rootRadius,
199 const double Kr, const double density, const double viscosity,
200 const std::size_t sourceId,
201 const double delta) const
202 {
203 ✗ const auto psBar = computePsBar(p0, pRoot, kernelRadius, rootRadius, Kr, density, viscosity, sourceId, delta);
204 ✗ return -2.0*M_PI*rootRadius*Kr*(psBar - pRoot)*density;
205 }
206
207 // try to be smart and use a cache to reduce nonlinear solves
208 298723402 double computeSource(std::size_t bEIdx, std::size_t lEIdx,
209 const double p0, const double pRoot,
210 const double kernelRadius, const double rootRadius,
211 const double Kr, const double density, const double viscosity,
212 const std::size_t sourceId,
213 const double delta) const
214 {
215
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 298723402 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 298723402 times.
 597446804 if (!sourceCache_)
216 ✗ return computeSource(p0, pRoot, kernelRadius, rootRadius, Kr, density, viscosity, sourceId, delta);
217
218 597446804 auto& cacheLine = (*sourceCache_)[lEIdx];
219
5/8
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 2681038 times.
✓ Branch 3 taken 5255630 times.
✓ Branch 4 taken 41436788 times.
✓ Branch 5 taken 41480565 times.
✓ Branch 6 taken 254605576 times.
✗ Branch 7 not taken.
 1241629803 auto bIt = std::find_if(cacheLine.begin(), cacheLine.end(), [&](const auto& p){ return p.first == bEIdx; });
220
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 298723402 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 298723402 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 298723402 times.
 896170206 if (bIt == cacheLine.end())
221 ✗ DUNE_THROW(Dune::InvalidStateException, "Source cache does not exist for index pair (l,b): " << lEIdx << "," << bEIdx);
222
223 303440080 return bIt->second(p0, pRoot, [&](const auto p0_, const auto pRoot_)
224 {
225 37733423 const auto psBar = computePsBar(p0_, pRoot_, kernelRadius, rootRadius, Kr, density, viscosity, sourceId, delta);
226 14150031 return -2.0*M_PI*rootRadius*Kr*(psBar - pRoot_)*density;
227 303440080 });
228 }
229
230 const Dumux::Detail::RootSoil::KirchhoffTransformationCached& transformation() const
231 { return transformation_; }
232
233 template<class CouplingManager>
234 2 void enableCache(const CouplingManager& couplingManager)
235 {
236 // create a cache vector for each low dim element
237 // the vector shall contain one element for each coupled bulk element (where p0 is measured)
238
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 const auto& lowDimGG = couplingManager.problem(CouplingManager::lowDimIdx).gridGeometry();
239
2/4
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 2 sourceCache_ = std::make_unique<CacheStorage>(lowDimGG.numDofs());
240 // lets use the stencil to find these coupled elements (probably only works for tpfa)
241
6/6
✓ Branch 0 taken 2340 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2340 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2340 times.
✓ Branch 5 taken 2 times.
 2342 for (std::size_t lowDimElementIdx = 0; lowDimElementIdx < sourceCache_->size(); ++lowDimElementIdx)
242 {
243 4680 const auto& couplingStencil = couplingManager.couplingStencils(CouplingManager::lowDimIdx).at(lowDimElementIdx);
244 9360 (*sourceCache_)[lowDimElementIdx].reserve(couplingStencil.size());
245
4/4
✓ Branch 0 taken 2658 times.
✓ Branch 1 taken 2340 times.
✓ Branch 2 taken 2658 times.
✓ Branch 3 taken 2340 times.
 12336 for (std::size_t bulkElementIdx : couplingStencil)
246 10632 (*sourceCache_)[lowDimElementIdx].emplace_back(bulkElementIdx, Detail::RootSoil::SourceCache{});
247 }
248 2 }
249 private:
250 std::unique_ptr<CacheStorage> sourceCache_;
251
252 Detail::RootSoil::KirchhoffTransformationCached transformation_;
253 double K_, porosity_;
254 double nonlinearSolverTolerance_;
255 };
256
257 } // end namespace Dumux::RootSoil
258
259 #endif
260