GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/multidomain/embedded/1d3d/root_soil_benchmark/couplingreconstruction.hh
Date: 2024-09-21 20:52:54
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 1 KirchhoffTransformation()
40 1 {
41 1 alpha = getParam<double>("Soil.SpatialParams.VanGenuchtenAlpha");
42 1 n = getParam<double>("Soil.SpatialParams.VanGenuchtenN");
43 1 l = getParam<double>("Soil.SpatialParams.VanGenuchtenL", 0.5);
44 1 m = 1.0 - 1.0/n;
45 1 }
46
47 38472785 inline double kr(const double pw) const
48 {
49 38472785 using std::pow;
50 38472785 const auto pc = 1.0e5-pw;
51 38472785 const auto swe = pow(pow(alpha*pc, n) + 1.0, -m);
52 38472785 const auto r = 1.0 - pow(1.0 - pow(swe, 1.0/m), m);
53 38472785 return pow(swe, l)*r*r;
54 }
55
56 //! Kirchhoff transform
57 inline double computeTransformed(const double p) const
58 19286393 { 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 1 KirchhoffTransformationCached()
81 1 {
82
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 pMin = getParam<double>("MixedDimension.Reconstruction.PressureMin", -2.6e7);
83
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 pMax = getParam<double>("MixedDimension.Reconstruction.PressureMax", 1.0e5);
84
1/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 1 samples = getParam<std::size_t>("MixedDimension.Reconstruction.Samples", 10000);
85
86
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 p_.resize(samples);
87
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 psi_.resize(samples);
88 1 const double pStep = (pMax-pMin)/samples;
89
2/2
✓ Branch 0 taken 100001 times.
✓ Branch 1 taken 1 times.
 100002 for (std::size_t i = 0; i <= samples; ++i)
90 {
91 100001 const auto pValue = pMin + i*pStep;
92 100001 p_[i] = pValue;
93 100001 psi_[i] = transformation_.computeTransformed(pValue);
94 }
95
96 1 std::cout << "[coupling] Presampled Kirchhoff transformation with "
97
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 2 << samples << " samples (step size: " << pStep << " Pa)."
98
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 << std::endl;
99 1 }
100
101 //! Kirchhoff transform
102 inline double computeTransformed(const double p) const
103 102574630 { return interpolate<InterpolationPolicy::LinearTable>(p, p_, psi_); }
104
105 //! Inverse Kirchhoff transform
106 inline double computeInverseTransformed(const double psi) const
107 14 { 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 1329 : q_{}
119 {
120 1329 p0_.fill(1e20);
121 1329 pRoot_.fill(1e20);
122 }
123
124 template<class F>
125 154753281 double operator() (const double p0, const double pRoot, const F& f)
126 {
127
2/2
✓ Branch 0 taken 323218140 times.
✓ Branch 1 taken 2438772 times.
 325656912 for (std::size_t i = 0; i < size; ++i)
128
12/12
✓ Branch 0 taken 244076898 times.
✓ Branch 1 taken 79141242 times.
✓ Branch 2 taken 244076898 times.
✓ Branch 3 taken 79141242 times.
✓ Branch 4 taken 244076898 times.
✓ Branch 5 taken 79141242 times.
✓ Branch 6 taken 152314509 times.
✓ Branch 7 taken 91762389 times.
✓ Branch 8 taken 152314509 times.
✓ Branch 9 taken 91762389 times.
✓ Branch 10 taken 152314509 times.
✓ Branch 11 taken 91762389 times.
 969654420 if (std::abs(p0-p0_[i]) < 1e-14 && std::abs(pRoot-pRoot_[i]) < 1e-14)
129 304629018 return q_[i];
130
131 2438772 const auto q = f(p0, pRoot);
132
2/2
✓ Branch 0 taken 1626696 times.
✓ Branch 1 taken 812075 times.
 2438771 p0_[writePos] = p0;
133
2/2
✓ Branch 0 taken 1626696 times.
✓ Branch 1 taken 812075 times.
 2438771 pRoot_[writePos] = pRoot;
134
2/2
✓ Branch 0 taken 1626696 times.
✓ Branch 1 taken 812075 times.
 2438771 q_[writePos] = q;
135
2/2
✓ Branch 0 taken 1626696 times.
✓ Branch 1 taken 812075 times.
 2438771 writePos = (writePos >= size-1) ? 0 : writePos+1;
136 2438771 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 1 CouplingReconstruction()
160 2 {
161
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 K_ = getParam<double>("Soil.SpatialParams.Permeability");
162
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 porosity_ = getParam<double>("Soil.SpatialParams.Porosity");
163 // relative tolerance in term of the unknown approximation (psBar)
164
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 nonlinearSolverTolerance_ = getParam<double>("MixedDimension.Reconstruction.RelativeTolerance", 1e-4);
165
166
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 2 std::cout << "\n" << "[coupling] Test Kirchhoff transform:\n";
167 const auto& tf = transformation_;
168
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 1 times.
 8 for (const auto p : {0.9e5, 0.5e5, 0.1e5, 0.0, -1.0e5, -5.0e5, -1.5e6})
169 {
170 7 const auto psi = tf.computeTransformed(p);
171 7 const auto pNum = tf.computeInverseTransformed(psi);
172
4/8
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 7 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 7 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 7 times.
✗ Branch 14 not taken.
 28 std::cout << "pw: " << p << ", psi: " << psi << ", pw (inverse): " << pNum << '\n';
173
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7 times.
 14 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 1 times.
✗ Branch 2 not taken.
 1 std::cout << std::endl;
177 1 }
178
179 2438772 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 2438772 const auto& tf = transformation_;
186 2438772 const auto psi0 = tf.computeTransformed(p0);
187 2438772 const auto sourceFactor = (std::log(kernelRadius/rootRadius) - 0.5 + 0.5*delta*delta/kernelRadius/kernelRadius)*viscosity/(2.0*M_PI*K_);
188
189 295529988 const auto residual = [&](const double psBar){
190 195394144 const auto sourceTerm = sourceFactor*2.0*M_PI*rootRadius*Kr*(psBar - pRoot);
191 48848536 return psi0 - tf.computeTransformed(psBar) - sourceTerm;
192 };
193
194 2438772 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 154753281 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 154753281 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 154753281 times.
 309506562 if (!sourceCache_)
216 ✗ return computeSource(p0, pRoot, kernelRadius, rootRadius, Kr, density, viscosity, sourceId, delta);
217
218 309506562 auto& cacheLine = (*sourceCache_)[lEIdx];
219
5/8
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 1388802 times.
✓ Branch 3 taken 2722467 times.
✓ Branch 4 taken 21466530 times.
✓ Branch 5 taken 21489217 times.
✓ Branch 6 taken 131897949 times.
✗ Branch 7 not taken.
 643224808 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 154753281 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 154753281 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 154753281 times.
 464259843 if (bIt == cacheLine.end())
221 ✗ DUNE_THROW(Dune::InvalidStateException, "Source cache does not exist for index pair (l,b): " << lEIdx << "," << bEIdx);
222
223 157192053 return bIt->second(p0, pRoot, [&](const auto p0_, const auto pRoot_)
224 {
225 19510175 const auto psBar = computePsBar(p0_, pRoot_, kernelRadius, rootRadius, Kr, density, viscosity, sourceId, delta);
226 7316313 return -2.0*M_PI*rootRadius*Kr*(psBar - pRoot_)*density;
227 157192053 });
228 }
229
230 const Dumux::Detail::RootSoil::KirchhoffTransformationCached& transformation() const
231 { return transformation_; }
232
233 template<class CouplingManager>
234 1 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 1 times.
 1 const auto& lowDimGG = couplingManager.problem(CouplingManager::lowDimIdx).gridGeometry();
239
2/4
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 1 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 1170 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1170 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1170 times.
✓ Branch 5 taken 1 times.
 1171 for (std::size_t lowDimElementIdx = 0; lowDimElementIdx < sourceCache_->size(); ++lowDimElementIdx)
242 {
243 2340 const auto& couplingStencil = couplingManager.couplingStencils(CouplingManager::lowDimIdx).at(lowDimElementIdx);
244 4680 (*sourceCache_)[lowDimElementIdx].reserve(couplingStencil.size());
245
4/4
✓ Branch 0 taken 1329 times.
✓ Branch 1 taken 1170 times.
✓ Branch 2 taken 1329 times.
✓ Branch 3 taken 1170 times.
 6168 for (std::size_t bulkElementIdx : couplingStencil)
246 5316 (*sourceCache_)[lowDimElementIdx].emplace_back(bulkElementIdx, Detail::RootSoil::SourceCache{});
247 }
248 1 }
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