| 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 RichardsTests | ||
| 10 | * \brief Helper functions for analytical solutions | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_RICHARDS_BENCHMARKS_ANALYTICAL_HH | ||
| 14 | #define DUMUX_RICHARDS_BENCHMARKS_ANALYTICAL_HH | ||
| 15 | |||
| 16 | #include <cmath> | ||
| 17 | #include <algorithm> | ||
| 18 | |||
| 19 | #include <dumux/io/format.hh> | ||
| 20 | #include <dumux/io/gnuplotinterface.hh> | ||
| 21 | |||
| 22 | #include <dumux/common/math.hh> | ||
| 23 | #include <dumux/common/parameters.hh> | ||
| 24 | #include <dumux/common/parameters.hh> | ||
| 25 | #include <dumux/common/integrate.hh> | ||
| 26 | |||
| 27 | #include <dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh> | ||
| 28 | #include <dumux/material/components/simpleh2o.hh> | ||
| 29 | |||
| 30 | namespace Dumux { | ||
| 31 | |||
| 32 | /*! | ||
| 33 | * \file | ||
| 34 | * \ingroup RichardsTests | ||
| 35 | * \brief Analytical solution for infiltration scenario | ||
| 36 | * \note Root-soil benchmark paper Schnepf et al. (case M2.1, Eq. 4) https://doi.org/10.3389/fpls.2020.00316 | ||
| 37 | * \note based on Vanderborght 2005 (see Fig. 4abc and Eq. 56-60) https://doi.org/10.2113/4.1.206 | ||
| 38 | * | ||
| 39 | * Assumption is that the infiltration front is sufficiently far | ||
| 40 | * away from the boundary. The boundary conditions for the analytical | ||
| 41 | * solution are prescribed at z = ±∞. | ||
| 42 | */ | ||
| 43 | class AnalyticalSolutionM21 | ||
| 44 | { | ||
| 45 | using PcKrSwCurve = FluidMatrix::VanGenuchtenNoReg<double>; | ||
| 46 | public: | ||
| 47 | 3 | AnalyticalSolutionM21() | |
| 48 |
7/18✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 3 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 3 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
|
6 | : pcKrSwCurve_("SpatialParams") |
| 49 | { | ||
| 50 | 3 | gravity_ = 9.81; | |
| 51 | 3 | density_ = Components::SimpleH2O<double>::liquidDensity(0,0); | |
| 52 | 3 | viscosity_ = Components::SimpleH2O<double>::liquidViscosity(0,0); | |
| 53 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | porosity_ = getParam<double>("SpatialParams.Porosity"); |
| 54 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | permeability_ = getParam<double>("SpatialParams.Permeability"); |
| 55 | |||
| 56 |
1/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
3 | const auto initialHead = getParam<double>("Problem.InitialHeadInCm", -400.0)*0.01; // cm to m |
| 57 | 3 | const auto pwInitial = initialHead*gravity_*density_ + 1.0e5; | |
| 58 | 3 | const auto pcInitial = 1.0e5 - pwInitial; | |
| 59 | 3 | const auto swIntial = pcKrSwCurve_.sw(pcInitial); | |
| 60 | |||
| 61 | // compute θ_i, K_i, θ_sur, K_sur, θ_a | ||
| 62 | 3 | waterContentInitial_ = swIntial*porosity_; | |
| 63 | 3 | conductivityInitial_ = conductivity(waterContentInitial_); | |
| 64 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | waterContentSurface_ = getParam<double>("Analytical.WaterContentSurface"); |
| 65 | 3 | conductivitySurface_ = conductivity(waterContentSurface_); | |
| 66 | 3 | waterContentRef_ = 0.5*(waterContentSurface_ + waterContentInitial_); | |
| 67 | |||
| 68 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
3 | std::cout << std::endl << "Computed analytical solution (Vanderborght 2005 4abc, Schnepf 2020 M2.1) with:\n"; |
| 69 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | const auto soilType = getParam<std::string>("SpatialParams.Name"); |
| 70 |
2/4✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
|
6 | std::cout << "Material: " << soilType << "\n"; |
| 71 |
3/8✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
3 | std::cout << Fmt::format("-> θ_i: {}\n-> θ_sur: {}\n", waterContentInitial_, waterContentSurface_); |
| 72 |
3/8✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
3 | std::cout << Fmt::format("-> θ_a: {}\n", waterContentRef_); |
| 73 |
3/10✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
3 | std::cout << Fmt::format("-> K_i: {}\n-> K_sur: {}\n", conductivityInitial_, conductivitySurface_); |
| 74 | |||
| 75 | // create a table for fast evaluation of the inverse | ||
| 76 | // for the inverse we need to go from high water content to low water content so | ||
| 77 | // that deltaEta is monotonically increasing which is a requirement for the value key | ||
| 78 | // for the table interpolation | ||
| 79 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
|
3 | waterContents_ = linspace(waterContentSurface_ - 1e-8, waterContentInitial_ + 1e-8, 10000); |
| 80 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
6 | deltaEtas_.resize(waterContents_.size()); |
| 81 | 3 | std::transform(waterContents_.begin(), waterContents_.end(), deltaEtas_.begin(), | |
| 82 | 12 | [this](auto theta){ return deltaEtaExact(theta); }); | |
| 83 | |||
| 84 | // plot delta eta over water content | ||
| 85 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
|
6 | GnuplotInterface<double> gnuplot(false); |
| 86 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | gnuplot.setOpenPlotWindow(false); |
| 87 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | gnuplot.resetPlot(); |
| 88 |
5/12✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
6 | gnuplot.setXlabel("water content"); |
| 89 |
5/12✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
6 | gnuplot.setYlabel("delta eta"); |
| 90 |
9/24✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 3 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 3 times.
✓ Branch 19 taken 3 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 3 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
|
12 | gnuplot.addDataSetToPlot(waterContents_, deltaEtas_, "theta_vs_deltaeta_" + soilType + ".dat", "w l t 'analytical'"); |
| 91 |
5/12✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
9 | gnuplot.plot("infiltration_theta_vs_deltaeta"); |
| 92 | |||
| 93 | // TODO how to choose these? Is there an exact value? | ||
| 94 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | waterContentRefDepth_ = getParam<double>("Analytical.WaterContentRefDepth"); |
| 95 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | waterContentRefTime_ = getParam<double>("Analytical.WaterContentRefTime"); |
| 96 |
3/8✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
6 | std::cout << Fmt::format("-> η_a: {}\n", eta(waterContentRefDepth_, waterContentRefTime_)); |
| 97 | 3 | } | |
| 98 | |||
| 99 | // analytical solution for given x and t | ||
| 100 | 2700 | double waterContent(double x /*in cm*/, double t /*in days*/) const | |
| 101 | { | ||
| 102 | 8100 | const auto deltaEta = eta(x, t) - eta(waterContentRefDepth_, waterContentRefTime_); | |
| 103 | 5400 | return interpolate<InterpolationPolicy::LinearTable>(deltaEta, deltaEtas_, waterContents_); | |
| 104 | } | ||
| 105 | |||
| 106 | double deltaEta(double x /*in cm*/, double t /*in days*/) const | ||
| 107 | 8100 | { return eta(x, t) - eta(waterContentRefDepth_, waterContentRefTime_); } | |
| 108 | |||
| 109 | ✗ | double initialWaterContent() const | |
| 110 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | { return waterContentInitial_; } |
| 111 | |||
| 112 | double surfaceWaterContent() const | ||
| 113 | { return waterContentSurface_; } | ||
| 114 | |||
| 115 | ✗ | double refWaterContent() const | |
| 116 | 9 | { return waterContentRef_; } | |
| 117 | |||
| 118 | auto waterContentAndDeltaEtaPlot() const | ||
| 119 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
6 | { return std::tie(waterContents_, deltaEtas_); } |
| 120 | |||
| 121 | private: | ||
| 122 | // hydraulic conductivity Kf (including relative permeability) in cm/day | ||
| 123 | 21123750 | double conductivity(double waterContent) const | |
| 124 | { | ||
| 125 | 21123750 | const auto sw = waterContent/porosity_; | |
| 126 | 21123750 | const auto krw = pcKrSwCurve_.krw(sw); | |
| 127 | 21123750 | return krw*permeability_*gravity_*density_/viscosity_*100*86400; | |
| 128 | } | ||
| 129 | |||
| 130 | // D = Kf*dh/dtheta Van Genuchten 1980 Eq. 10 in cm^2/day | ||
| 131 | 10561872 | double diffusivity(double waterContent) const | |
| 132 | { | ||
| 133 | // pc = -h*gravity_*density_ | ||
| 134 | // h = -pc/(gravity_*density_) | ||
| 135 | 10561872 | const auto sw = waterContent/porosity_; | |
| 136 | 10561872 | const auto dh_dpc = -1.0/(0.01*gravity_*density_); | |
| 137 | 10561872 | const auto dpc_dsw = pcKrSwCurve_.dpc_dsw(sw); | |
| 138 | 10561872 | const auto dsw_dtheta = 1.0/porosity_; | |
| 139 | 10561872 | const auto dh_dtheta = dh_dpc*dpc_dsw*dsw_dtheta; | |
| 140 | 10561872 | return conductivity(waterContent)*dh_dtheta; | |
| 141 | } | ||
| 142 | |||
| 143 | // integral travelling wave (Eq. 4 Schnepf 2020) | ||
| 144 | double integral(double waterContent) const | ||
| 145 | { | ||
| 146 | 31715616 | const auto integrand = [this](auto theta){ | |
| 147 | 10561872 | return diffusivity(theta)/( | |
| 148 | 10561872 | (conductivitySurface_ - conductivityInitial_)*(theta - waterContentInitial_) | |
| 149 | 10561872 | - (conductivity(theta) - conductivityInitial_)*(waterContentSurface_ - waterContentInitial_) | |
| 150 | 10561872 | ); | |
| 151 | 30000 | }; | |
| 152 | |||
| 153 | 60000 | return integrateScalarFunction(integrand, waterContent, waterContentRef_); | |
| 154 | } | ||
| 155 | |||
| 156 | // relative transformed coordinate position (Eq. 4 Schnepf 2020) | ||
| 157 | double deltaEtaExact(double waterContent) const | ||
| 158 | { | ||
| 159 | 30000 | return (waterContentSurface_ - waterContentInitial_)*integral(waterContent); | |
| 160 | } | ||
| 161 | |||
| 162 | // transformed coordinate | ||
| 163 | double eta(double x /*in cm*/, double t /*in days*/) const | ||
| 164 | { | ||
| 165 |
1/2✓ Branch 1 taken 2700 times.
✗ Branch 2 not taken.
|
5403 | return std::fabs(x) - (conductivitySurface_ - conductivityInitial_)/(waterContentSurface_ - waterContentInitial_)*t; |
| 166 | } | ||
| 167 | |||
| 168 | PcKrSwCurve pcKrSwCurve_; | ||
| 169 | double waterContentSurface_, waterContentInitial_, waterContentRef_; | ||
| 170 | double waterContentRefDepth_, waterContentRefTime_; | ||
| 171 | double conductivitySurface_, conductivityInitial_; | ||
| 172 | double porosity_, permeability_, viscosity_, density_, gravity_; | ||
| 173 | |||
| 174 | std::vector<double> waterContents_; | ||
| 175 | std::vector<double> deltaEtas_; | ||
| 176 | }; | ||
| 177 | |||
| 178 | /*! | ||
| 179 | * \file | ||
| 180 | * \ingroup RichardsTests | ||
| 181 | * \brief Analytical solution for evaporation scenario | ||
| 182 | * \note Root-soil benchmark paper Schnepf et al. (case M2.2) https://doi.org/10.3389/fpls.2020.00316 | ||
| 183 | * \note based on Vanderborght 2005 (see Fig. 5abcd and Eq. 39-47) https://doi.org/10.2113/4.1.206 | ||
| 184 | * \note Derivation in Lockington 1994 https://doi.org/10.1029/93WR03411 | ||
| 185 | * | ||
| 186 | * The analytical solution is an approximative solution for an infinite domain | ||
| 187 | * and is neglects graviational effects. Lockington 1994 estimates less than 1% error | ||
| 188 | * against a given exact solution (without gravity) for a simplified pc-sw relationship. | ||
| 189 | */ | ||
| 190 | class AnalyticalSolutionM22 | ||
| 191 | { | ||
| 192 | using PcKrSwCurve = FluidMatrix::VanGenuchtenNoReg<double>; | ||
| 193 | public: | ||
| 194 | 4 | AnalyticalSolutionM22() | |
| 195 |
4/10✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
8 | : pcKrSwCurve_("SpatialParams") |
| 196 | { | ||
| 197 | 4 | gravity_ = 9.81; | |
| 198 | 4 | density_ = Components::SimpleH2O<double>::liquidDensity(0,0); | |
| 199 | 4 | viscosity_ = Components::SimpleH2O<double>::liquidViscosity(0,0); | |
| 200 | 4 | porosity_ = getParam<double>("SpatialParams.Porosity"); | |
| 201 | 4 | permeability_ = getParam<double>("SpatialParams.Permeability"); | |
| 202 | |||
| 203 | 4 | const auto initialHead = getParam<double>("Problem.InitialHeadInCm", -40.0)*0.01; // cm to m | |
| 204 | 4 | potentialEvaporationRate_ = getParam<double>("Problem.SurfaceFluxMilliMeterPerDay"); | |
| 205 | |||
| 206 | 4 | const auto pwInitial = initialHead*gravity_*density_ + 1.0e5; | |
| 207 | 4 | const auto pcInitial = 1.0e5 - pwInitial; | |
| 208 | 4 | const auto swIntial = pcKrSwCurve_.sw(pcInitial); | |
| 209 | 4 | waterContentInitial_ = swIntial*porosity_; | |
| 210 | |||
| 211 | 4 | const auto pwSurface = -10000*0.01*gravity_*density_ + 1.0e5; | |
| 212 | 4 | const auto pcSurface = 1.0e5 - pwSurface; | |
| 213 | 4 | const auto swSurface = pcKrSwCurve_.sw(pcSurface); | |
| 214 | 4 | waterContentSurface_ = swSurface*porosity_; | |
| 215 | |||
| 216 | 4 | std::cout << std::endl << "Computed analytical solution (Vanderborght 2005 5abcd, Schnepf 2020 M2.2) with:\n"; | |
| 217 | 4 | const auto soilType = getParam<std::string>("SpatialParams.Name"); | |
| 218 |
2/4✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
|
8 | std::cout << "Material: " << soilType << "\n"; |
| 219 |
3/8✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
4 | std::cout << Fmt::format("-> θ_i: {}\n-> θ_sur: {}\n", waterContentInitial_, waterContentSurface_); |
| 220 | |||
| 221 | // water content from "effective" water content, Vanderborght 2005 Eq. 40 | ||
| 222 | 4 | const auto theta = [this](auto thetaEff){ | |
| 223 | 4152 | return thetaEff*(waterContentInitial_ - waterContentSurface_) + waterContentSurface_; | |
| 224 | 4 | }; | |
| 225 | |||
| 226 | 4 | const auto dIntegral = integrateScalarFunction( | |
| 227 | 676 | [&,this](auto thetaEff){ return diffusivity(theta(thetaEff)); }, 0.0, 1.0 | |
| 228 | 4 | ); | |
| 229 | |||
| 230 |
3/8✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
4 | std::cout << Fmt::format("D(0), D(1): {}, {}\n", diffusivity(waterContentSurface_), diffusivity(waterContentInitial_)); |
| 231 |
3/8✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
4 | std::cout << Fmt::format("Kf(0), Kf(1): {}, {}\n", conductivity(waterContentSurface_), conductivity(waterContentInitial_)); |
| 232 |
3/8✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
4 | std::cout << Fmt::format("-> Integral D(Θ): {}\n", dIntegral); |
| 233 | |||
| 234 | 4 | const auto dThetaIntegral = integrateScalarFunction( | |
| 235 | 772 | [&,this](auto thetaEff){ return thetaEff*diffusivity(theta(thetaEff)); }, 0.0, 1.0 | |
| 236 | 4 | ); | |
| 237 | |||
| 238 |
1/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
4 | std::cout << Fmt::format("0D(0), 0.5D(0.5) 1D(1): {}, {}, {}\n", |
| 239 |
2/4✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
|
8 | 0.0*diffusivity(theta(0.0)), 0.5*diffusivity(theta(0.5)), 1.0*diffusivity(theta(1.0))); |
| 240 |
3/8✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
4 | std::cout << Fmt::format("-> Integral Θ D(Θ): {}\n", dThetaIntegral); |
| 241 | |||
| 242 | // Eq. 43 Vanderborght 2005 | ||
| 243 | 4 | const auto beta = dThetaIntegral*dThetaIntegral/(dIntegral*dIntegral); | |
| 244 | |||
| 245 | 4 | const auto dThetaFactorIntegral = integrateScalarFunction( | |
| 246 | 2504 | [&,this](auto thetaEff){ | |
| 247 | 1252 | const auto weight = (1.0-beta*thetaEff); | |
| 248 | 1252 | return weight*weight*diffusivity(theta(thetaEff)); | |
| 249 | }, 0.0, 1.0 | ||
| 250 | 4 | ); | |
| 251 | |||
| 252 |
3/8✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
4 | std::cout << Fmt::format("-> Integral (1.0 - βΘ)² D(θ): {}\n", dThetaFactorIntegral); |
| 253 | |||
| 254 | // Eq. 42 Vanderborght 2005 | ||
| 255 | 4 | const auto alpha = dThetaFactorIntegral/dIntegral; | |
| 256 | |||
| 257 | // Eq. 41 Vanderborght 2005 | ||
| 258 | 4 | const auto bracketFactor = 1.0 - alpha/((1.0 - beta)*(1.0 - beta)); | |
| 259 | 4 | const auto mu = 3*beta*(1.0 + std::sqrt(1.0 - 14.0/9.0*bracketFactor)) / ( 2*(beta - 1.0)*bracketFactor ); | |
| 260 | |||
| 261 |
3/8✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
4 | std::cout << Fmt::format("-> µ: {}, α: {}, β: {}\n", mu, alpha, beta); |
| 262 | |||
| 263 | // Eq. 39 Vanderborght 2005 | ||
| 264 | 4 | desorptivity_ = (waterContentInitial_ - waterContentSurface_)*std::sqrt(4.0/mu*dIntegral); | |
| 265 | |||
| 266 |
3/8✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
4 | std::cout << Fmt::format("-> Desorptivity, Sw(θ_i, θ_sur): {}\n", desorptivity_); |
| 267 | |||
| 268 | // Eq. 44 & 45 Vanderborght 2005 | ||
| 269 | 4 | tPotential_ = desorptivity_*desorptivity_/(2.0*potentialEvaporationRate_*0.1*potentialEvaporationRate_*0.1); | |
| 270 | 4 | tPrime_ = 0.5*tPotential_; | |
| 271 | |||
| 272 |
3/10✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
4 | std::cout << Fmt::format("-> Critical time (t_pot): {} days\n", tPotential_); |
| 273 | |||
| 274 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
|
4 | std::cout << std::endl; |
| 275 | 4 | } | |
| 276 | |||
| 277 | // Evaporationrate in mm/day (Eq. 46 & 47 Vanderborght 2005) | ||
| 278 | 931 | double evaporationRate(double t /*in days*/) const | |
| 279 | { | ||
| 280 |
2/2✓ Branch 0 taken 99 times.
✓ Branch 1 taken 832 times.
|
931 | if (t < tPotential_) |
| 281 | 99 | return potentialEvaporationRate_; | |
| 282 | else | ||
| 283 | 832 | return 0.5*desorptivity_/std::sqrt(tPrime_ + t - tPotential_)*10; | |
| 284 | } | ||
| 285 | |||
| 286 | private: | ||
| 287 | // hydraulic conductivity Kf (including relative permeability) in cm/day | ||
| 288 | 4168 | double conductivity(double waterContent) const | |
| 289 | { | ||
| 290 | 4168 | const auto sw = waterContent/porosity_; | |
| 291 | 4168 | const auto krw = pcKrSwCurve_.krw(sw); | |
| 292 | 4168 | return krw*permeability_*gravity_*density_/viscosity_*100*86400; | |
| 293 | } | ||
| 294 | |||
| 295 | // D = Kf*dh/dtheta Van Genuchten 1980 Eq. 10 in cm^2/day | ||
| 296 | 4160 | double diffusivity(double waterContent) const | |
| 297 | { | ||
| 298 | // pc = -h*gravity_*density_ | ||
| 299 | // h = -pc/(gravity_*density_) | ||
| 300 | 4160 | const auto sw = waterContent/porosity_; | |
| 301 | 4160 | const auto dh_dpc = -1.0/(0.01*gravity_*density_); | |
| 302 | 4160 | const auto dpc_dsw = pcKrSwCurve_.dpc_dsw(sw); | |
| 303 | 4160 | const auto dsw_dtheta = 1.0/porosity_; | |
| 304 | 4160 | const auto dh_dtheta = dh_dpc*dpc_dsw*dsw_dtheta; | |
| 305 | 4160 | return conductivity(waterContent)*dh_dtheta; | |
| 306 | } | ||
| 307 | |||
| 308 | PcKrSwCurve pcKrSwCurve_; | ||
| 309 | double waterContentSurface_, waterContentInitial_; | ||
| 310 | double desorptivity_; | ||
| 311 | double tPrime_, tPotential_; | ||
| 312 | double potentialEvaporationRate_; | ||
| 313 | double porosity_, permeability_, viscosity_, density_, gravity_; | ||
| 314 | }; | ||
| 315 | |||
| 316 | } // end namespace Dumux | ||
| 317 | |||
| 318 | #endif | ||
| 319 |