diff --git a/dumux/common/doubleexpintegrationconstants.hh b/dumux/common/doubleexpintegrationconstants.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9c5fdaee4d6795bfe8ffd1b473970bff24fd3969
--- /dev/null
+++ b/dumux/common/doubleexpintegrationconstants.hh
@@ -0,0 +1,498 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ * *
+ * This version is modified after the original version by John D. Cook *
+ * see https://www.codeproject.com/ *
+ * Articles/31550/Fast-Numerical-Integration *
+ * which is licensed under BSD-2-clause, which reads as follows: *
+ * Copyright John D. Cook *
+ * *
+ * Redistribution and use in source and binary forms, with or without *
+ * modification, are permitted provided that the following *
+ * conditions are met: *
+ * 1. Redistributions of source code must retain the above *
+ * copyright notice, this list of conditions *
+ and the following disclaimer. *
+ * 2. Redistributions in binary form must reproduce the above *
+ * copyright notice, this list of conditions *
+ * and the following disclaimer *
+ * in the documentation and/or other materials *
+ * provided with the distribution. *
+ * *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT *
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS *
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE *
+ * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, *
+ * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL *
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE *
+ * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS *
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, *
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE *
+ * OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, *
+ * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
+ * *
+ *****************************************************************************/
+
+#ifndef DUMUX_COMMON_DOUBLEEXP_INTEGRATION_CONSTANTS_HH
+#define DUMUX_COMMON_DOUBLEEXP_INTEGRATION_CONSTANTS_HH
+
+/*
+The double exponential rule is based on the observation that the trapezoid rule converges
+very rapidly for functions on the entire real line that go to zero like exp( - exp(t) ).
+The change of variables x = tanh( pi sinh(t) /2) transforms an integral over [-1, 1]
+into an integral with integrand suited to the double exponential rule.
+
+The transformed integral is infinite, but we truncate the domain of integration to [-3, 3].
+The limit '3' was chosen for two reasons: for t = 3, the transformed x values
+are nearly equal to 1 for 12 or more significant figures. Also, for t = 3, the
+smallest weights are 12 orders of magnitude smaller than the largest weights; setting
+the cutoff larger than 3 would not have a significant impact on the integral value
+unless there is a strong singularity at one of the end points.
+
+The change of variables x(t) is an odd function, i.e. x(-t) = -x(t), and so we need only
+store the positive x values. Also, the derivative w(t) = x'(t) is even, i.e. w(-t) = w(t),
+and so we need only store the weights corresponding to positive values of x.
+
+The integration first applies the trapezoid rule to [-3, 3] in steps of size 1.
+Then it subsequently cuts the step size in half each time, comparing the results.
+Integration stops when subsequent iterations are close enough together or the maximum
+integration points have been used.
+By cutting h in half, the previous integral can be reused; we only need evaluate the
+integrand at the newly added points.
+
+Finally, note that we're not strictly using the trapezoid rule: we don't treat the
+end points differently. This is because we assume the values at the ends of the interval
+hardly matter due to the rapid decay of the integrand.
+
+All values below were calculated with Mathematica.
+*/
+
+namespace Dumux {
+
+constexpr double doubleExponentialIntegrationAbcissas[] =
+{
+ // 1st layer abcissas: transformed 0, 1, 2, 3
+ 0.00000000000000000000,
+ 0.95136796407274694573,
+ 0.99997747719246159286,
+ 0.99999999999995705839,
+ // 2nd layer abcissas: transformed 1/2, 3/2, 5/2
+ 0.67427149224843582608,
+ 0.99751485645722438683,
+ 0.99999998887566488198,
+ // 3rd layer abcissas: transformed 1/4, 3/4, ...
+ 0.37720973816403417379,
+ 0.85956905868989663517,
+ 0.98704056050737689169,
+ 0.99968826402835320905,
+ 0.99999920473711471266,
+ 0.99999999995285644818,
+ // 4th layer abcissas: transformed 1/8, 3/8, ...
+ 0.19435700332493543161,
+ 0.53914670538796776905,
+ 0.78060743898320029925,
+ 0.91487926326457461091,
+ 0.97396686819567744856,
+ 0.99405550663140214329,
+ 0.99906519645578584642,
+ 0.99990938469514399984,
+ 0.99999531604122052843,
+ 0.99999989278161241838,
+ 0.99999999914270509218,
+ 0.99999999999823216531,
+ // 5th layer abcissa: transformed 1/16, 3/16, ...
+ 0.097923885287832333262,
+ 0.28787993274271591456,
+ 0.46125354393958570440,
+ 0.61027365750063894488,
+ 0.73101803479256151149,
+ 0.82331700550640237006,
+ 0.88989140278426019808,
+ 0.93516085752198468323,
+ 0.96411216422354729193,
+ 0.98145482667733517003,
+ 0.99112699244169880223,
+ 0.99610866543750854254,
+ 0.99845420876769773751,
+ 0.99945143443527460584,
+ 0.99982882207287494166,
+ 0.99995387100562796075,
+ 0.99998948201481850361,
+ 0.99999801714059543208,
+ 0.99999969889415261122,
+ 0.99999996423908091534,
+ 0.99999999678719909830,
+ 0.99999999978973286224,
+ 0.99999999999039393352,
+ 0.99999999999970809734,
+ // 6th layer abcissas: transformed 1/32, 3/32, ...
+ 0.049055967305077886315,
+ 0.14641798429058794053,
+ 0.24156631953888365838,
+ 0.33314226457763809244,
+ 0.41995211127844715849,
+ 0.50101338937930910152,
+ 0.57558449063515165995,
+ 0.64317675898520470128,
+ 0.70355000514714201566,
+ 0.75669390863372994941,
+ 0.80279874134324126576,
+ 0.84221924635075686382,
+ 0.87543539763040867837,
+ 0.90301328151357387064,
+ 0.92556863406861266645,
+ 0.94373478605275715685,
+ 0.95813602271021369012,
+ 0.96936673289691733517,
+ 0.97797623518666497298,
+ 0.98445883116743083087,
+ 0.98924843109013389601,
+ 0.99271699719682728538,
+ 0.99517602615532735426,
+ 0.99688031812819187372,
+ 0.99803333631543375402,
+ 0.99879353429880589929,
+ 0.99928111192179195541,
+ 0.99958475035151758732,
+ 0.99976797159956083506,
+ 0.99987486504878034648,
+ 0.99993501992508242369,
+ 0.99996759306794345976,
+ 0.99998451990227082442,
+ 0.99999293787666288565,
+ 0.99999693244919035751,
+ 0.99999873547186590954,
+ 0.99999950700571943689,
+ 0.99999981889371276701,
+ 0.99999993755407837378,
+ 0.99999997987450320175,
+ 0.99999999396413420165,
+ 0.99999999832336194826,
+ 0.99999999957078777261,
+ 0.99999999989927772326,
+ 0.99999999997845533741,
+ 0.99999999999582460688,
+ 0.99999999999927152627,
+ 0.99999999999988636130,
+ // 7th layer abcissas: transformed 1/64, 3/64, ...
+ 0.024539763574649160379,
+ 0.073525122985671294475,
+ 0.12222912220155764235,
+ 0.17046797238201051811,
+ 0.21806347346971200463,
+ 0.26484507658344795046,
+ 0.31065178055284596083,
+ 0.35533382516507453330,
+ 0.39875415046723775644,
+ 0.44078959903390086627,
+ 0.48133184611690504422,
+ 0.52028805069123015958,
+ 0.55758122826077823080,
+ 0.59315035359195315880,
+ 0.62695020805104287950,
+ 0.65895099174335012438,
+ 0.68913772506166767176,
+ 0.71750946748732412721,
+ 0.74407838354734739913,
+ 0.76886868676824658459,
+ 0.79191549237614211447,
+ 0.81326360850297385168,
+ 0.83296629391941087564,
+ 0.85108400798784873261,
+ 0.86768317577564598669,
+ 0.88283498824466895513,
+ 0.89661425428007602579,
+ 0.90909831816302043511,
+ 0.92036605303195280235,
+ 0.93049693799715340631,
+ 0.93957022393327475539,
+ 0.94766419061515309734,
+ 0.95485549580502268541,
+ 0.96121861515111640753,
+ 0.96682537031235585284,
+ 0.97174454156548730892,
+ 0.97604156025657673933,
+ 0.97977827580061576265,
+ 0.98301279148110110558,
+ 0.98579936302528343597,
+ 0.98818835380074264243,
+ 0.99022624046752774694,
+ 0.99195566300267761562,
+ 0.99341551316926403900,
+ 0.99464105571251119672,
+ 0.99566407681695316965,
+ 0.99651305464025377317,
+ 0.99721334704346870224,
+ 0.99778739195890653083,
+ 0.99825491617199629344,
+ 0.99863314864067747762,
+ 0.99893703483351217373,
+ 0.99917944893488591716,
+ 0.99937140114093768690,
+ 0.99952223765121720422,
+ 0.99963983134560036519,
+ 0.99973076151980848263,
+ 0.99980048143113838630,
+ 0.99985347277311141171,
+ 0.99989338654759256426,
+ 0.99992317012928932869,
+ 0.99994518061445869309,
+ 0.99996128480785666613,
+ 0.99997294642523223656,
+ 0.99998130127012072679,
+ 0.99998722128200062811,
+ 0.99999136844834487344,
+ 0.99999423962761663478,
+ 0.99999620334716617675,
+ 0.99999752962380516793,
+ 0.99999841381096473542,
+ 0.99999899541068996962,
+ 0.99999937270733536947,
+ 0.99999961398855024275,
+ 0.99999976602333243312,
+ 0.99999986037121459941,
+ 0.99999991800479471056,
+ 0.99999995264266446185,
+ 0.99999997311323594362,
+ 0.99999998500307631173,
+ 0.99999999178645609907,
+ 0.99999999558563361584,
+ 0.99999999767323673790,
+ 0.99999999879798350040,
+ 0.99999999939177687583,
+ 0.99999999969875436925,
+ 0.99999999985405611550,
+ 0.99999999993088839501,
+ 0.99999999996803321674,
+ 0.99999999998556879008,
+ 0.99999999999364632387,
+ 0.99999999999727404948,
+ 0.99999999999886126543,
+ 0.99999999999953722654,
+ 0.99999999999981720098,
+ 0.99999999999992987953
+}; // end abcissas
+
+constexpr double doubleExponentialIntegrationWeights[] =
+{
+ // First layer weights
+ 1.5707963267948966192,
+ 0.230022394514788685,
+ 0.00026620051375271690866,
+ 1.3581784274539090834e-12,
+ // 2nd layer weights
+ 0.96597657941230114801,
+ 0.018343166989927842087,
+ 2.1431204556943039358e-7,
+ // 3rd layer weights
+ 1.3896147592472563229,
+ 0.53107827542805397476,
+ 0.076385743570832304188,
+ 0.0029025177479013135936,
+ 0.000011983701363170720047,
+ 1.1631165814255782766e-9,
+ // 4th layer weights
+ 1.5232837186347052132,
+ 1.1934630258491569639,
+ 0.73743784836154784136,
+ 0.36046141846934367417,
+ 0.13742210773316772341,
+ 0.039175005493600779072,
+ 0.0077426010260642407123,
+ 0.00094994680428346871691,
+ 0.000062482559240744082891,
+ 1.8263320593710659699e-6,
+ 1.8687282268736410132e-8,
+ 4.9378538776631926964e-11,
+ // 5th layer weights
+ 1.5587733555333301451,
+ 1.466014426716965781,
+ 1.297475750424977998,
+ 1.0816349854900704074,
+ 0.85017285645662006895,
+ 0.63040513516474369106,
+ 0.44083323627385823707,
+ 0.290240679312454185,
+ 0.17932441211072829296,
+ 0.10343215422333290062,
+ 0.055289683742240583845,
+ 0.027133510013712003219,
+ 0.012083543599157953493,
+ 0.0048162981439284630173,
+ 0.0016908739981426396472,
+ 0.00051339382406790336017,
+ 0.00013205234125609974879,
+ 0.000028110164327940134749,
+ 4.8237182032615502124e-6,
+ 6.4777566035929719908e-7,
+ 6.5835185127183396672e-8,
+ 4.8760060974240625869e-9,
+ 2.5216347918530148572e-10,
+ 8.6759314149796046502e-12,
+ // 6th layer weights
+ 1.5677814313072218572,
+ 1.5438811161769592204,
+ 1.4972262225410362896,
+ 1.4300083548722996676,
+ 1.3452788847662516615,
+ 1.2467012074518577048,
+ 1.1382722433763053734,
+ 1.0240449331118114483,
+ 0.90787937915489531693,
+ 0.79324270082051671787,
+ 0.68306851634426375464,
+ 0.57967810308778764708,
+ 0.48475809121475539287,
+ 0.39938474152571713515,
+ 0.32408253961152890402,
+ 0.258904639514053516,
+ 0.20352399885860174519,
+ 0.15732620348436615027,
+ 0.11949741128869592428,
+ 0.089103139240941462841,
+ 0.065155533432536205042,
+ 0.046668208054846613644,
+ 0.032698732726609031113,
+ 0.022379471063648476483,
+ 0.014937835096050129696,
+ 0.0097072237393916892692,
+ 0.0061300376320830301252,
+ 0.0037542509774318343023,
+ 0.0022250827064786427022,
+ 0.0012733279447082382027,
+ 0.0007018595156842422708,
+ 0.00037166693621677760301,
+ 0.00018856442976700318572,
+ 0.000091390817490710122732,
+ 0.000042183183841757600604,
+ 0.000018481813599879217116,
+ 7.6595758525203162562e-6,
+ 2.9916615878138787094e-6,
+ 1.0968835125901264732e-6,
+ 3.7595411862360630091e-7,
+ 1.1992442782902770219e-7,
+ 3.5434777171421953043e-8,
+ 9.6498888961089633609e-9,
+ 2.4091773256475940779e-9,
+ 5.482835779709497755e-10,
+ 1.1306055347494680536e-10,
+ 2.0989335404511469109e-11,
+ 3.4841937670261059685e-12,
+ // 7th layer weights
+ 1.5700420292795931467,
+ 1.5640214037732320999,
+ 1.5520531698454121192,
+ 1.5342817381543034316,
+ 1.5109197230741697127,
+ 1.48224329788553807,
+ 1.4485862549613225916,
+ 1.4103329714462590129,
+ 1.3679105116808964881,
+ 1.3217801174437728579,
+ 1.2724283455378627082,
+ 1.2203581095793582207,
+ 1.1660798699324345766,
+ 1.1101031939653403796,
+ 1.0529288799552666556,
+ 0.99504180404613271514,
+ 0.93690461274566793366,
+ 0.87895234555278212039,
+ 0.82158803526696470334,
+ 0.7651792989089561367,
+ 0.71005590120546898385,
+ 0.65650824613162753076,
+ 0.60478673057840362158,
+ 0.55510187800363350959,
+ 0.5076251588319080997,
+ 0.4624903980553677613,
+ 0.41979566844501548066,
+ 0.37960556938665160999,
+ 0.3419537959230168323,
+ 0.30684590941791694932,
+ 0.27426222968906810637,
+ 0.24416077786983990868,
+ 0.21648020911729617038,
+ 0.19114268413342749532,
+ 0.16805663794826916233,
+ 0.14711941325785693248,
+ 0.12821973363120098675,
+ 0.11123999898874453035,
+ 0.096058391865189467849,
+ 0.082550788110701737654,
+ 0.070592469906866999352,
+ 0.060059642358636300319,
+ 0.05083075757257047107,
+ 0.042787652157725676034,
+ 0.035816505604196436523,
+ 0.029808628117310126969,
+ 0.024661087314753282511,
+ 0.020277183817500123926,
+ 0.016566786254247575375,
+ 0.013446536605285730674,
+ 0.010839937168255907211,
+ 0.0086773307495391815854,
+ 0.0068957859690660035329,
+ 0.0054388997976239984331,
+ 0.0042565295990178580165,
+ 0.0033044669940348302363,
+ 0.0025440657675291729678,
+ 0.0019418357759843675814,
+ 0.0014690143599429791058,
+ 0.0011011261134519383862,
+ 0.00081754101332469493115,
+ 0.00060103987991147422573,
+ 0.00043739495615911687786,
+ 0.00031497209186021200274,
+ 0.00022435965205008549104,
+ 0.00015802788400701191949,
+ 0.00011002112846666697224,
+ 0.000075683996586201477788,
+ 0.000051421497447658802092,
+ 0.0000344921247593431977,
+ 0.000022832118109036146591,
+ 0.000014908514031870608449,
+ 9.5981941283784710776e-6,
+ 6.0899100320949039256e-6,
+ 3.8061983264644899045e-6,
+ 2.3421667208528096843e-6,
+ 1.4183067155493917523e-6,
+ 8.4473756384859863469e-7,
+ 4.9458288702754198508e-7,
+ 2.8449923659159806339e-7,
+ 1.6069394579076224911e-7,
+ 8.9071395140242387124e-8,
+ 4.8420950198072369669e-8,
+ 2.579956822953589238e-8,
+ 1.3464645522302038796e-8,
+ 6.8784610955899001111e-9,
+ 3.4371856744650090511e-9,
+ 1.6788897682161906807e-9,
+ 8.0099784479729665356e-10,
+ 3.7299501843052790038e-10,
+ 1.6939457789411646876e-10,
+ 7.4967397573818224522e-11,
+ 3.230446433325236576e-11,
+ 1.3542512912336274432e-11,
+ 5.5182369468174885821e-12,
+ 2.1835922099233609052e-12
+}; // end weights
+
+} // end namespace Dumux
+
+#endif // include guard
diff --git a/dumux/common/doubleexpintegrator.hh b/dumux/common/doubleexpintegrator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..693abc3b045146803e9ac755b834115732779f5c
--- /dev/null
+++ b/dumux/common/doubleexpintegrator.hh
@@ -0,0 +1,206 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ * *
+ * This version is modified after the original version by John D. Cook *
+ * see https://www.codeproject.com/ *
+ * Articles/31550/Fast-Numerical-Integration *
+ * which is licensed under BSD-2-clause, which reads as follows: *
+ * Copyright John D. Cook *
+ * *
+ * Redistribution and use in source and binary forms, with or without *
+ * modification, are permitted provided that the following *
+ * conditions are met: *
+ * 1. Redistributions of source code must retain the above *
+ * copyright notice, this list of conditions *
+ and the following disclaimer. *
+ * 2. Redistributions in binary form must reproduce the above *
+ * copyright notice, this list of conditions *
+ * and the following disclaimer *
+ * in the documentation and/or other materials *
+ * provided with the distribution. *
+ * *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT *
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS *
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE *
+ * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, *
+ * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL *
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE *
+ * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS *
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, *
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE *
+ * OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, *
+ * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
+ * *
+ *****************************************************************************/
+
+#ifndef DUMUX_COMMON_DOUBLEEXP_INTEGRATOR_HH
+#define DUMUX_COMMON_DOUBLEEXP_INTEGRATOR_HH
+
+#include <cmath>
+#include <limits>
+
+#include <dune/common/std/type_traits.hh>
+#include <dumux/common/doubleexpintegrationconstants.hh>
+
+namespace Dumux {
+
+/*!
+ * \brief Numerical integration in one dimension using the double exponential method of M. Mori.
+ */
+template<class Scalar>
+class DoubleExponentialIntegrator
+{
+public:
+ /*!
+ * \brief Integrate an analytic function over a finite interval
+ * \param f the integrand (invocable with a single scalar)
+ * \param a lower limit of integration
+ * \param b upper limit of integration
+ * \param targetAbsoluteError desired bound on error
+ * \param numFunctionEvaluations number of function evaluations used
+ * \param f errorEstimate estimate for error in integration
+ * \return The value of the integral
+ */
+ template<class Function,
+ typename std::enable_if_t<Dune::Std::is_invocable<Function, Scalar>::value>...>
+ static Scalar integrate(const Function& f,
+ const Scalar a,
+ const Scalar b,
+ const Scalar targetAbsoluteError,
+ int& numFunctionEvaluations,
+ Scalar& errorEstimate)
+ {
+ // Apply the linear change of variables x = ct + d
+ // $$\int_a^b f(x) dx = c \int_{-1}^1 f( ct + d ) dt$$
+ // c = (b-a)/2, d = (a+b)/2
+
+ const Scalar c = 0.5*(b - a);
+ const Scalar d = 0.5*(a + b);
+ return integrateCore_(f, c, d, targetAbsoluteError, numFunctionEvaluations, errorEstimate,
+ doubleExponentialIntegrationAbcissas, doubleExponentialIntegrationWeights);
+ }
+
+ /*!
+ * \brief Integrate an analytic function over a finite interval.
+ * \note This version overloaded to not require arguments passed in for
+ * function evaluation counts or error estimates.
+ * \param f the integrand (invocable with a single scalar)
+ * \param a lower integral bound
+ * \param b upper integral bound
+ * \param targetAbsoluteError desired absolute error in the result
+ * \return The value of the integral.
+ */
+ template<class Function,
+ typename std::enable_if_t<Dune::Std::is_invocable<Function, Scalar>::value>...>
+ static Scalar integrate(const Function& f,
+ const Scalar a,
+ const Scalar b,
+ const Scalar targetAbsoluteError)
+ {
+ int numFunctionEvaluations;
+ Scalar errorEstimate;
+ return integrate(f, a, b, targetAbsoluteError, numFunctionEvaluations, errorEstimate);
+ }
+
+private:
+ // Integrate f(cx + d) with the given integration constants
+ template<class Function>
+ static Scalar integrateCore_(const Function& f,
+ const Scalar c, // slope of change of variables
+ const Scalar d, // intercept of change of variables
+ Scalar targetAbsoluteError,
+ int& numFunctionEvaluations,
+ Scalar& errorEstimate,
+ const double* abcissas,
+ const double* weights)
+ {
+ targetAbsoluteError /= c;
+
+ // Offsets to where each level's integration constants start.
+ // The last element is not a beginning but an end.
+ static const int offsets[] = {1, 4, 7, 13, 25, 49, 97, 193};
+ static const int numLevels = sizeof(offsets)/sizeof(int) - 1;
+
+ Scalar newContribution = 0.0;
+ Scalar integral = 0.0;
+ Scalar h = 1.0;
+ errorEstimate = std::numeric_limits<Scalar>::max();
+ Scalar previousDelta, currentDelta = std::numeric_limits<Scalar>::max();
+
+ integral = f(c*abcissas[0] + d) * weights[0];
+ int i;
+ for (i = offsets[0]; i != offsets[1]; ++i)
+ integral += weights[i]*(f(c*abcissas[i] + d) + f(-c*abcissas[i] + d));
+
+ for (int level = 1; level != numLevels; ++level)
+ {
+ h *= 0.5;
+ newContribution = 0.0;
+ for (i = offsets[level]; i != offsets[level+1]; ++i)
+ newContribution += weights[i]*(f(c*abcissas[i] + d) + f(-c*abcissas[i] + d));
+ newContribution *= h;
+
+ // difference in consecutive integral estimates
+ previousDelta = currentDelta;
+ using std::abs;
+ currentDelta = abs(0.5*integral - newContribution);
+ integral = 0.5*integral + newContribution;
+
+ // Once convergence kicks in, error is approximately squared at each step.
+ // Determine whether we're in the convergent region by looking at the trend in the error.
+ if (level == 1)
+ continue; // previousDelta meaningless, so cannot check convergence.
+
+ // Exact comparison with zero is harmless here. Could possibly be replaced with
+ // a small positive upper limit on the size of currentDelta, but determining
+ // that upper limit would be difficult. At worse, the loop is executed more
+ // times than necessary. But no infinite loop can result since there is
+ // an upper bound on the loop variable.
+ if (currentDelta == 0.0)
+ break;
+
+ using std::log;
+ const Scalar rate = log( currentDelta )/log( previousDelta ); // previousDelta != 0 or would have been kicked out previously
+
+ if (rate > 1.9 && rate < 2.1)
+ {
+ // If convergence theory applied perfectly, r would be 2 in the convergence region.
+ // r close to 2 is good enough. We expect the difference between this integral estimate
+ // and the next one to be roughly delta^2.
+ errorEstimate = currentDelta*currentDelta;
+ }
+ else
+ {
+ // Not in the convergence region. Assume only that error is decreasing.
+ errorEstimate = currentDelta;
+ }
+
+ if (errorEstimate < 0.1*targetAbsoluteError)
+ break;
+ }
+
+ numFunctionEvaluations = 2*i - 1;
+ errorEstimate *= c;
+ return c*integral;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/geometry/CMakeLists.txt b/dumux/common/geometry/CMakeLists.txt
index fd5972a41baebb36223d8a7cb7fbdd7756423c19..77968db349d5c6834f0dfa6a7e1054fe8bebd64d 100644
--- a/dumux/common/geometry/CMakeLists.txt
+++ b/dumux/common/geometry/CMakeLists.txt
@@ -9,5 +9,6 @@ intersectingentities.hh
intersectspointgeometry.hh
intersectspointsimplex.hh
makegeometry.hh
+refinementquadraturerule.hh
triangulation.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/common/geometry)
diff --git a/dumux/common/geometry/refinementquadraturerule.hh b/dumux/common/geometry/refinementquadraturerule.hh
new file mode 100644
index 0000000000000000000000000000000000000000..114ec0960ad3bc70f8c813bf5df5f4db44bac39e
--- /dev/null
+++ b/dumux/common/geometry/refinementquadraturerule.hh
@@ -0,0 +1,94 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Geometry
+ * \author Timo Koch
+ * \brief A quadrature based on refinement
+ */
+
+#ifndef DUMUX_COMMON_REFINEMENT_QUADRATURERULE_HH
+#define DUMUX_COMMON_REFINEMENT_QUADRATURERULE_HH
+
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/virtualrefinement.hh>
+#include <dumux/common/math.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Geometry
+ * \brief A "quadrature" based on virtual refinement
+ */
+template<typename ct, int mydim>
+class RefinementQuadratureRule : public Dune::QuadratureRule<ct, mydim>
+{
+public:
+ //! The space dimension
+ static constexpr int dim = mydim;
+
+ //! The highest quadrature order available
+ static constexpr int highest_order = 1;
+
+ ~RefinementQuadratureRule() final {}
+
+ RefinementQuadratureRule(Dune::GeometryType type, int levels)
+ : Dune::QuadratureRule<ct, dim>(type, highest_order)
+ {
+ auto& refinement = Dune::buildRefinement<dim, ct>(type, type);
+ const auto tag = Dune::refinementLevels(levels);
+
+ // get the vertices
+ const auto numVertices = refinement.nVertices(tag);
+ std::vector<Dune::FieldVector<ct, dim>> points; points.reserve(numVertices);
+ auto vIt = refinement.vBegin(tag);
+ const auto vItEnd = refinement.vEnd(tag);
+ for (; vIt != vItEnd; ++vIt)
+ points.emplace_back(vIt.coords());
+
+ // go over all elements and get center and volume
+ const auto numElements = refinement.nElements(tag);
+ this->reserve(numElements);
+ auto eIt = refinement.eBegin(tag);
+ const auto eItEnd = refinement.eEnd(tag);
+ for (; eIt != eItEnd; ++eIt)
+ {
+ // quadrature point in the centroid of the sub-element and weight is its volume
+ const auto weight = computeVolume_(type, points, eIt.vertexIndices());
+ this->emplace_back(eIt.coords(), weight);
+ }
+ }
+
+private:
+ template<class VertexIndices>
+ ct computeVolume_(Dune::GeometryType t, const std::vector<Dune::FieldVector<ct, dim>>& points, const VertexIndices& indices) const
+ {
+ if (t != Dune::GeometryTypes::simplex(2))
+ DUNE_THROW(Dune::NotImplemented, "Only implemented for 2d simplices");
+
+ const auto ab = points[indices[1]] - points[indices[0]];
+ const auto ac = points[indices[2]] - points[indices[0]];
+ using std::abs;
+ return abs(crossProduct(ab, ac))*0.5;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/integrate.hh b/dumux/common/integrate.hh
index f458f763655027d37d1d3e093c01d8ad33f482e4..d2ece0d8b6b3b32672790bd87c128da80a161da4 100644
--- a/dumux/common/integrate.hh
+++ b/dumux/common/integrate.hh
@@ -29,11 +29,13 @@
#include <dune/geometry/quadraturerules.hh>
#include <dune/common/concept.hh>
+#include <dune/common/std/type_traits.hh>
#if HAVE_DUNE_FUNCTIONS
#include <dune/functions/gridfunctions/gridfunction.hh>
#endif
+#include <dumux/common/doubleexpintegrator.hh>
#include <dumux/discretization/evalsolution.hh>
#include <dumux/discretization/elementsolution.hh>
#include <dumux/common/typetraits/typetraits.hh>
@@ -251,6 +253,24 @@ auto integrateL2Error(const GridView& gv,
}
#endif
+/*!
+ * \brief Integrate a scalar function
+ * \param f the integrand (invocable with a single scalar)
+ * \param a lower integral bound
+ * \param b upper integral bound
+ * \param targetAbsoluteError desired absolute error in the result
+ * \return The value of the integral
+ */
+template<class Scalar, class Function,
+ typename std::enable_if_t<Dune::Std::is_invocable<Function, Scalar>::value>...>
+Scalar integrateScalarFunction(const Function& f,
+ const Scalar lowerBound,
+ const Scalar upperBound,
+ const Scalar targetAbsoluteError = 1e-13)
+{
+ return DoubleExponentialIntegrator<Scalar>::integrate(f, lowerBound, upperBound, targetAbsoluteError);
+}
+
} // end namespace Dumux
#endif
diff --git a/dumux/nonlinear/CMakeLists.txt b/dumux/nonlinear/CMakeLists.txt
index 92640a44c6e3e6da21a55c9a677c213957358e24..2c5c2050a30e68e34df4e6f3a9f6f50fcf28b3db 100644
--- a/dumux/nonlinear/CMakeLists.txt
+++ b/dumux/nonlinear/CMakeLists.txt
@@ -1,4 +1,5 @@
install(FILES
+findscalarroot.hh
newtonconvergencewriter.hh
newtonsolver.hh
staggerednewtonconvergencewriter.hh
diff --git a/dumux/nonlinear/findscalarroot.hh b/dumux/nonlinear/findscalarroot.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7c230d47fe9f5f9f57a10b3dc6349ea31235eca9
--- /dev/null
+++ b/dumux/nonlinear/findscalarroot.hh
@@ -0,0 +1,206 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Common
+ * \brief Root finding algorithms for scalar functions
+ */
+#ifndef DUMUX_COMMON_SCALAR_ROOT_FINDING_HH
+#define DUMUX_COMMON_SCALAR_ROOT_FINDING_HH
+
+#include <cmath>
+#include <limits>
+
+#include <dune/common/std/type_traits.hh>
+
+#include <dumux/common/exceptions.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/numericdifferentiation.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Common
+ * \brief Newton's root finding algorithm for scalar functions (secant method)
+ * \param xOld initial guess
+ * \param residual Residual function
+ * \param derivative Derivative of the residual
+ * \param tol Relative shift tolerance
+ * \param maxIter Maximum number of iterations
+ * \todo replace std::result_of by std::invoke_result for C++17
+ */
+template<class Scalar, class ResFunc, class DerivFunc,
+ typename std::enable_if_t<Dune::Std::is_invocable<ResFunc, Scalar>::value &&
+ Dune::Std::is_invocable<DerivFunc, Scalar>::value>...>
+Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual, const DerivFunc& derivative,
+ const Scalar tol = 1e-13, const int maxIter = 200)
+{
+ Scalar xNew = xOld;
+ Scalar r = residual(xNew);
+
+ int n = maxIter;
+ Scalar relativeShift = std::numeric_limits<Scalar>::max();
+ while (relativeShift > tol && n > 0)
+ {
+ xNew = xOld - r/derivative(xOld);
+ r = residual(xNew);
+
+ using std::abs; using std::max;
+ relativeShift = abs(xOld-xNew)/max(abs(xOld), abs(xNew));
+ xOld = xNew;
+ n--;
+ }
+
+ using std::isfinite;
+ if (!isfinite(r))
+ DUNE_THROW(NumericalProblem, "Residual is not finite: " << r << " after " << maxIter - n << " iterations!");
+
+ if (relativeShift > tol)
+ DUNE_THROW(NumericalProblem, "Scalar newton solver did not converge after " << maxIter << " iterations!");
+
+ return xNew;
+}
+
+/*!
+ * \ingroup Common
+ * \brief Newton's root finding algorithm for scalar functions (secant method)
+ * \note The derivative is numerically computed. If the derivative is know use signature with derivative function.
+ * \param xOld initial guess
+ * \param residual Residual function
+ * \param tol Relative shift tolerance
+ * \param maxIter Maximum number of iterations
+ */
+template<class Scalar, class ResFunc,
+ typename std::enable_if_t<Dune::Std::is_invocable<ResFunc, Scalar>::value>...>
+Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual,
+ const Scalar tol = 1e-13, const int maxIter = 200)
+{
+ const auto eps = NumericDifferentiation::epsilon(xOld);
+ auto derivative = [&](const auto x){ return (residual(x + eps)-residual(x))/eps; };
+ return findScalarRootNewton(xOld, residual, derivative, tol, maxIter);
+}
+
+/*!
+ * \ingroup Common
+ * \brief Brent's root finding algorithm for scalar functions
+ * \note Modified from pseudo-code on wikipedia: https://en.wikipedia.org/wiki/Brent%27s_method
+ * \note See also R.P. Brent "An algorithm with guaranteed convergence for finding a zero of a function", The Computer Journal (1971).
+ * \note This is usually more robust than Newton's method
+ * \param a Lower bound
+ * \param b Upper bound
+ * \param residual Residual function
+ * \param tol Relative shift tolerance
+ * \param maxIter Maximum number of iterations
+ */
+template<class Scalar, class ResFunc,
+ typename std::enable_if_t<Dune::Std::is_invocable<ResFunc, Scalar>::value>...>
+Scalar findScalarRootBrent(Scalar a, Scalar b, const ResFunc& residual,
+ const Scalar tol = 1e-13, const int maxIter = 200)
+{
+ // precompute the residuals (minimize function evaluations)
+ Scalar fa = residual(a);
+ Scalar fb = residual(b);
+ Scalar fs = 0.0;
+
+ // check if the root is inside the given interval
+ using std::signbit;
+ if (!signbit(fa*fb))
+ DUNE_THROW(NumericalProblem, "Brent's algorithm failed: [a,b] does not contain any, or no uniquely defined root!");
+
+ // sort bounds
+ using std::abs; using std::swap;
+ if (abs(fa) < abs(fb))
+ {
+ swap(a, b);
+ swap(fa, fb);
+ }
+
+ Scalar c = a;
+ Scalar fc = fa;
+ Scalar d = 0.0;
+ Scalar s = 0.0;
+ bool flag = true;
+
+ for (int i = 0; i < maxIter; ++i)
+ {
+ // stopping criterion
+ using std::max;
+ if (abs(b-a) < tol*max(abs(a), abs(b)))
+ return b;
+
+ // inverse quadratic interpolation
+ if (fa != fc && fb != fc)
+ {
+ const auto fab = fa-fb;
+ const auto fac = fa-fc;
+ const auto fbc = fb-fc;
+ s = a*fb*fc/(fab*fac) - b*fa*fc/(fab*fbc) + c*fa*fb/(fac*fbc);
+ }
+
+ // secant method
+ else
+ {
+ s = b - fb*(b-a)/(fb-fa);
+ }
+
+ // bisection method
+ if ( (s < (3*a + b)*0.25 || s > b)
+ || (flag && abs(s-b) >= abs(b-c)*0.5)
+ || (!flag && abs(s-b) >= abs(c-d)*0.5)
+ || (flag && abs(b-c) < tol*max(abs(b), abs(c)))
+ || (!flag && abs(c-d) < tol*max(abs(c), abs(d))) )
+ {
+ s = (a+b)*0.5;
+ flag = true;
+ }
+ else
+ flag = false;
+
+ // compute residual at new guess
+ fs = residual(s);
+ d = c;
+ c = b;
+ fc = fb;
+
+ // check on which side of the root s falls
+ if (fa*fs < 0.0)
+ {
+ b = s;
+ fb = fs;
+ }
+ else
+ {
+ a = s;
+ fa = fs;
+ }
+
+ // sort if necessary
+ if (abs(fa) < abs(fb))
+ {
+ swap(a, b);
+ swap(fa, fb);
+ }
+ }
+
+ DUNE_THROW(NumericalProblem, "Brent's algorithm didn't converge after " << maxIter << " iterations!");
+}
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 416332236d6dd28c43906415eda98911378d9bb3..b5c7428a4b7e7b3edb40d5f60c4fd3972c74014d 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -4,5 +4,6 @@ add_subdirectory(freeflow)
add_subdirectory(io)
add_subdirectory(material)
add_subdirectory(multidomain)
+add_subdirectory(nonlinear)
add_subdirectory(porousmediumflow)
add_subdirectory(discretization)
diff --git a/test/common/CMakeLists.txt b/test/common/CMakeLists.txt
index 96a770ad978cb059e815d44e68eb29eb45bab6ab..80146db8c5a314ffa848258085b72f327ca779b1 100644
--- a/test/common/CMakeLists.txt
+++ b/test/common/CMakeLists.txt
@@ -1,6 +1,7 @@
add_subdirectory(boundingboxtree)
add_subdirectory(functions)
add_subdirectory(geometry)
+add_subdirectory(integrate)
add_subdirectory(math)
add_subdirectory(parameters)
add_subdirectory(propertysystem)
@@ -8,4 +9,4 @@ add_subdirectory(spline)
add_subdirectory(timeloop)
add_subdirectory(typetraits)
-dumux_add_test(SOURCES test_partial.cc)
+dumux_add_test(SOURCES test_partial.cc LABELS unit)
diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt
index f1b5b270e8401ebd8378ea580afbdcd0209fe0df..c2e2c5c6c713490e9ef35f6d98a3e371408db367 100644
--- a/test/common/geometry/CMakeLists.txt
+++ b/test/common/geometry/CMakeLists.txt
@@ -6,3 +6,4 @@ dumux_add_test(SOURCES test_2d2d_intersection.cc LABELS unit)
dumux_add_test(SOURCES test_2d3d_intersection.cc LABELS unit)
dumux_add_test(SOURCES test_graham_convex_hull.cc LABELS unit)
dumux_add_test(SOURCES test_makegeometry.cc LABELS unit)
+dumux_add_test(SOURCES test_refinementquadraturerule.cc LABELS unit)
diff --git a/test/common/geometry/test_refinementquadraturerule.cc b/test/common/geometry/test_refinementquadraturerule.cc
new file mode 100644
index 0000000000000000000000000000000000000000..7511e351a67b87cdce9bd712de69bb261cee1955
--- /dev/null
+++ b/test/common/geometry/test_refinementquadraturerule.cc
@@ -0,0 +1,41 @@
+#include <config.h>
+
+#include <iostream>
+#include <array>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/float_cmp.hh>
+#include <dumux/common/geometry/refinementquadraturerule.hh>
+#include <dune/geometry/affinegeometry.hh>
+
+int main (int argc, char *argv[]) try
+{
+ using namespace Dumux;
+
+ // make a simplex
+ const auto simplex = Dune::AffineGeometry<double, 2, 3>(
+ Dune::GeometryTypes::simplex(2),
+ std::array<Dune::FieldVector<double, 3>, 3>{{{0.0, 0.0, 1.0}, {2.0, 3.0, 4.0}, {3.0, -1.0, 1.0}}}
+ );
+
+ // make quadrature rule
+ const auto quad = RefinementQuadratureRule<double, 2>(Dune::GeometryTypes::simplex(2), /*refinementlevels=*/3);
+
+ // integrate volume
+ double volume = 0.0;
+ for (const auto& qp : quad)
+ volume += qp.weight()*simplex.integrationElement(qp.position());
+
+ if (Dune::FloatCmp::ne(volume, simplex.volume()))
+ DUNE_THROW(Dune::Exception, "Integration yields wrong simplex volume: " << volume << ", should be " << simplex.volume());
+
+ return 0;
+}
+// //////////////////////////////////
+// Error handler
+// /////////////////////////////////
+catch (const Dune::Exception& e) {
+ std::cout << e << std::endl;
+ return 1;
+}
diff --git a/test/common/integrate/CMakeLists.txt b/test/common/integrate/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..54f255b05100a9714fe6b8ecd2fbbbb253feacaa
--- /dev/null
+++ b/test/common/integrate/CMakeLists.txt
@@ -0,0 +1 @@
+dumux_add_test(SOURCES test_integratescalar.cc LABELS unit)
diff --git a/test/common/integrate/test_integratescalar.cc b/test/common/integrate/test_integratescalar.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9841040ef4b9d6b10782b1f36d8dcd5458dee16f
--- /dev/null
+++ b/test/common/integrate/test_integratescalar.cc
@@ -0,0 +1,90 @@
+#include <config.h>
+
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/float_cmp.hh>
+#include <dumux/common/integrate.hh>
+#include <dumux/nonlinear/findscalarroot.hh>
+
+namespace Dumux {
+
+class VanGenuchtenKrw
+{
+public:
+ VanGenuchtenKrw(double alpha, double n, double l = 0.5)
+ : alpha_(alpha), n_(n), m_(1.0 - 1.0/n), l_(l)
+ {}
+
+ double operator()(double pw) const
+ {
+ using std::pow;
+ const auto pc = 1.0e5-pw;
+ const auto swe = pow(pow(alpha_*pc, n_) + 1.0, -m_);
+ const auto r = 1.0 - pow(1.0 - pow(swe, 1.0/m_), m_);
+ return pow(swe, l_)*r*r;
+ }
+
+private:
+ double alpha_, n_, m_, l_;
+};
+
+} // end namespace Dumux
+
+int main(int argc, char* argv[]) try
+{
+ using namespace Dumux;
+
+ {
+ // test some integral where we know the analytical expression
+ const auto func = [](const double x){ return std::exp(2.0*x); };
+ const auto exactIntegral = 0.5*(std::exp(10.0)-std::exp(-2.0));
+ const auto numIntegral = integrateScalarFunction(func, -1.0, 5.0, 1e-8);
+ if (Dune::FloatCmp::ne<double, Dune::FloatCmp::CmpStyle::absolute>(exactIntegral, numIntegral, 1e-8))
+ DUNE_THROW(Dune::Exception, "Wrong integral: " << numIntegral << ", should be " << exactIntegral);
+ }
+
+ {
+ // test some integral where we know the analytical expression
+ const auto func = [](const double x){ return pow(1.0 - x, 5.0)*pow(x, -1.0/3.0); };
+ const auto exactIntegral = 2187.0/5236.0;
+ const auto numIntegral = integrateScalarFunction(func, 0.0, 1.0, 1e-8);
+ if (Dune::FloatCmp::ne<double, Dune::FloatCmp::CmpStyle::absolute>(exactIntegral, numIntegral, 1e-8))
+ DUNE_THROW(Dune::Exception, "Wrong integral: " << numIntegral << ", should be " << exactIntegral);
+ }
+
+ {
+ // relative permeability function of water pressure
+ const auto krFunc = VanGenuchtenKrw(/*alpha=*/1e-4, /*n=*/3.0);
+
+ // introduce a transformed variable u(p) = int_{-1.5e6}^{p} kr(q)dq
+ const auto kirchhoffTransform = [&](const double p){
+ return integrateScalarFunction(krFunc, -1.5e6, p, 1e-16);
+ };
+
+ // compute the inverse of the transformation numerically
+ const auto inverseTransform = [&](const auto u){
+ const auto residual = [&, u](const auto& p){ return u - kirchhoffTransform(p); };
+ return findScalarRootBrent(-1.6e6, 1.0e5, residual, 1e-10, 200000);
+ };
+
+ // compute transformation and inverse for several pressures
+ for (const auto p : {0.9e5, 0.1e5, 0.0, -5.0e5, -1.5e6})
+ {
+ const auto u = kirchhoffTransform(p);
+ const auto pNum = inverseTransform(u);
+ if (Dune::FloatCmp::ne<double, Dune::FloatCmp::CmpStyle::absolute>(p, pNum, 0.001))
+ DUNE_THROW(Dune::Exception, "Inverse did not recover original pressure value: "
+ << std::setprecision(15) << pNum << ", should be " << p);
+ }
+ }
+
+ return 0;
+}
+catch (const Dune::Exception& e)
+{
+ std::cout << e << std::endl;
+ return 1;
+}
diff --git a/test/nonlinear/CMakeLists.txt b/test/nonlinear/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3640235bdcfbdd1b4fa0aacce2527f8cc8114bf9
--- /dev/null
+++ b/test/nonlinear/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(findscalarroot)
diff --git a/test/nonlinear/findscalarroot/CMakeLists.txt b/test/nonlinear/findscalarroot/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..004545b738f41d14f80c12ced54202de6cafaf80
--- /dev/null
+++ b/test/nonlinear/findscalarroot/CMakeLists.txt
@@ -0,0 +1,2 @@
+dumux_add_test(SOURCES test_findscalarroot.cc
+ LABELS unit nonlinear)
diff --git a/test/nonlinear/findscalarroot/test_findscalarroot.cc b/test/nonlinear/findscalarroot/test_findscalarroot.cc
new file mode 100644
index 0000000000000000000000000000000000000000..88dcbea99255a45c4abeb97747ddd4b0c75db6bb
--- /dev/null
+++ b/test/nonlinear/findscalarroot/test_findscalarroot.cc
@@ -0,0 +1,50 @@
+#include <config.h>
+
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/float_cmp.hh>
+#include <dumux/nonlinear/findscalarroot.hh>
+
+int main(int argc, char* argv[]) try
+{
+ using namespace Dumux;
+
+ const auto func = [](const double x){ return x*x - 5.0; };
+ const auto deriv = [](const double x){ return 2*x; };
+ const auto absTol = 1e-15;
+ const auto relTol = 1e-15;
+
+ const auto rootNewton = findScalarRootNewton(5.0, func, deriv, absTol);
+ const auto rootNewtonNumDiff = findScalarRootNewton(5.0, func, absTol);
+ const auto rootBrent = findScalarRootBrent(0.0, 5.0, func, relTol);
+ const auto rootExact = std::sqrt(5.0);
+
+ if (Dune::FloatCmp::ne(rootExact, rootNewton, absTol))
+ DUNE_THROW(Dune::Exception, "Wrong root for Newton algorithm (exact derivative)! Relative error: " << std::abs(rootExact-rootNewton)/rootExact);
+ if (Dune::FloatCmp::ne(rootExact, rootNewtonNumDiff, absTol))
+ DUNE_THROW(Dune::Exception, "Wrong root for Newton algorithm (numeric differentiation)! Relative error: " << std::abs(rootExact-rootNewtonNumDiff)/rootExact);
+ if (Dune::FloatCmp::ne(rootExact, rootBrent, absTol))
+ DUNE_THROW(Dune::Exception, "Wrong root for Brent algorithm! Relative error: " << std::abs(rootExact-rootBrent)/rootExact);
+
+ std::cout << "Positive root of f(x) = (x*x - 5) is\n" << std::setprecision(15)
+ << "-- Exact: " << rootExact << "\n"
+ << "-- Newton (exact derivative): " << rootNewton << "\n"
+ << "-- Newton (numeric derivative): " << rootNewtonNumDiff << "\n"
+ << "-- Brent: " << rootBrent << std::endl;
+
+ return 0;
+
+}
+catch (const Dune::Exception& e)
+{
+ std::cout << e << std::endl;
+ return 1;
+}
+catch (...)
+{
+ std::cout << "Unknown exception thrown!" << std::endl;
+ return 1;
+}