3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
staggeredupwindmethods.hh
Go to the documentation of this file.
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 2 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_UPWINDING_METHODS_HH
26#define DUMUX_UPWINDING_METHODS_HH
27
28#include <cmath>
29#include <functional>
30#include <iostream>
31#include <string>
32
35
36namespace Dumux {
37
42enum class TvdApproach
43{
44 none, uniform, li, hou
45};
46
52{
54};
55
60template<class Scalar, int upwindSchemeOrder>
62{
63public:
64 StaggeredUpwindMethods(const std::string& paramGroup = "")
65 {
66 upwindWeight_ = getParamFromGroup<Scalar>(paramGroup, "Flux.UpwindWeight");
67
68 if (upwindSchemeOrder > 1)
69 {
70 // Read the runtime parameters
71 tvdApproach_ = tvdApproachFromString(getParamFromGroup<std::string>(paramGroup, "Flux.TvdApproach", "Uniform"));
72 differencingScheme_ = differencingSchemeFromString(getParamFromGroup<std::string>(paramGroup, "Flux.DifferencingScheme", "Minmod"));
73
74 // Assign the limiter_ depending on the differencing scheme
75 switch (differencingScheme_)
76 {
78 {
79 limiter_ = this->vanleer;
80 break;
81 }
83 {
84 if (tvdApproach_ == TvdApproach::hou)
85 DUNE_THROW(ParameterException, "\nDifferencing scheme (Van Albada) is not implemented for the Tvd approach (Hou).");
86 else
87 limiter_ = this->vanalbada;
88 break;
89 }
91 {
92 limiter_ = this->minmod;
93 break;
94 }
96 {
97 limiter_ = this->superbee;
98 break;
99 }
101 {
102 limiter_ = this->umist;
103 break;
104 }
106 {
107 limiter_ = this->mclimiter;
108 break;
109 }
111 {
112 limiter_ = this->wahyd;
113 break;
114 }
115 default:
116 {
117 DUNE_THROW(ParameterException, "\nDifferencing scheme " << static_cast<std::string>(differencingSchemeToString(differencingScheme_)) <<
118 " is not implemented.\n"); // should never be reached
119 break;
120 }
121 }
122
123 if (!hasParamInGroup(paramGroup, "Flux.TvdApproach"))
124 {
125 std::cout << "No TvdApproach specified. Defaulting to the Uniform method." << "\n";
126 std::cout << "Other available TVD approaches for uniform (and nonuniform) grids are as follows: \n"
127 << " " << tvdApproachToString(TvdApproach::uniform) << ": assumes a Uniform cell size distribution\n"
128 << " " << tvdApproachToString(TvdApproach::li) << ": Li's approach for nonuniform cell sizes\n"
129 << " " << tvdApproachToString(TvdApproach::hou) << ": Hou's approach for nonuniform cell sizes \n";
130 std::cout << "Each approach can be specified as written above in the Flux group under the title TvdApproach in your input file. \n";
131 }
132 if (!hasParamInGroup(paramGroup, "Flux.DifferencingScheme"))
133 {
134
135 std::cout << "No DifferencingScheme specified. Defaulting to the Minmod scheme." << "\n";
136 std::cout << "Other available Differencing Schemes are as follows: \n"
137 << " " << differencingSchemeToString(DifferencingScheme::vanleer) << ": The Vanleer flux limiter\n"
138 << " " << differencingSchemeToString(DifferencingScheme::vanalbada) << ": The Vanalbada flux limiter\n"
139 << " " << differencingSchemeToString(DifferencingScheme::minmod) << ": The Minmod flux limiter\n"
140 << " " << differencingSchemeToString(DifferencingScheme::superbee) << ": The Superbee flux limiter\n"
141 << " " << differencingSchemeToString(DifferencingScheme::umist) << ": The Umist flux limiter\n"
142 << " " << differencingSchemeToString(DifferencingScheme::mclimiter) << ": The Mclimiter flux limiter\n"
143 << " " << differencingSchemeToString(DifferencingScheme::wahyd) << ": The Wahyd flux limiter";
144 std::cout << "Each scheme can be specified as written above in the Flux group under the variable DifferencingScheme in your input file. \n";
145 }
146
147 std::cout << "Using the tvdApproach \"" << tvdApproachToString(tvdApproach_)
148 << "\" and the differencing Scheme \" " << differencingSchemeToString(differencingScheme_) << "\" \n";
149 }
150 else
151 {
152 // If the runtime parameters are not specified we will use upwind
153 tvdApproach_ = TvdApproach::none;
154 differencingScheme_ = DifferencingScheme::none;
155 }
156 }
157
163 {
164 if (tvd == "Uniform") return TvdApproach::uniform;
165 if (tvd == "Li") return TvdApproach::li;
166 if (tvd == "Hou") return TvdApproach::hou;
167 DUNE_THROW(ParameterException, "\nThis tvd approach : \"" << tvd << "\" is not implemented.\n"
168 << "The available TVD approaches for uniform (and nonuniform) grids are as follows: \n"
169 << tvdApproachToString(TvdApproach::uniform) << ": assumes a uniform cell size distribution\n"
170 << tvdApproachToString(TvdApproach::li) << ": li's approach for nonuniform cell sizes\n"
171 << tvdApproachToString(TvdApproach::hou) << ": hou's approach for nonuniform cell sizes");
172 }
173
178 {
179 switch (tvd)
180 {
181 case TvdApproach::uniform: return "Uniform";
182 case TvdApproach::li: return "Li";
183 case TvdApproach::hou: return "Hou";
184 default: return "Invalid"; // should never be reached
185 }
186 }
187
193 {
194 if (differencingScheme == "Vanleer") return DifferencingScheme::vanleer;
195 if (differencingScheme == "Vanalbada") return DifferencingScheme::vanalbada;
196 if (differencingScheme == "Minmod") return DifferencingScheme::minmod;
197 if (differencingScheme == "Superbee") return DifferencingScheme::superbee;
198 if (differencingScheme == "Umist") return DifferencingScheme::umist;
199 if (differencingScheme == "Mclimiter") return DifferencingScheme::mclimiter;
200 if (differencingScheme == "Wahyd") return DifferencingScheme::wahyd;
201 DUNE_THROW(ParameterException, "\nThis differencing scheme: \"" << differencingScheme << "\" is not implemented.\n"
202 << "The available differencing schemes are as follows: \n"
203 << differencingSchemeToString(DifferencingScheme::vanleer) << ": The Vanleer flux limiter\n"
204 << differencingSchemeToString(DifferencingScheme::vanalbada) << ": The VanAlbada flux limiter\n"
205 << differencingSchemeToString(DifferencingScheme::minmod) << ": The Min-Mod flux limiter\n"
206 << differencingSchemeToString(DifferencingScheme::superbee) << ": The SuperBEE flux limiter\n"
207 << differencingSchemeToString(DifferencingScheme::umist) << ": The UMist flux limiter\n"
208 << differencingSchemeToString(DifferencingScheme::mclimiter) << ": The McLimiter flux limiter\n"
209 << differencingSchemeToString(DifferencingScheme::wahyd) << ": The Wahyd flux limiter");
210 }
211
216 {
217 switch (differencingScheme)
218 {
219 case DifferencingScheme::vanleer: return "Vanleer";
220 case DifferencingScheme::vanalbada: return "Vanalbada";
221 case DifferencingScheme::minmod: return "Minmod";
222 case DifferencingScheme::superbee: return "Superbee";
223 case DifferencingScheme::umist: return "Umist";
224 case DifferencingScheme::mclimiter: return "Mclimiter";
225 case DifferencingScheme::wahyd: return "Wahyd";
226 default: return "Invalid"; // should never be reached
227 }
228 }
229
233 Scalar upwind(const Scalar downstreamMomentum,
234 const Scalar upstreamMomentum) const
235 {
236 return (upwindWeight_ * upstreamMomentum + (1.0 - upwindWeight_) * downstreamMomentum);
237 }
238
243 Scalar tvd(const std::array<Scalar,3>& momenta,
244 const std::array<Scalar,3>& distances,
245 const bool selfIsUpstream,
246 const TvdApproach tvdApproach) const
247 {
248 Scalar momentum = 0.0;
249 switch(tvdApproach)
250 {
252 {
253 momentum += tvdUniform(momenta, distances, selfIsUpstream);
254 break;
255 }
256 case TvdApproach::li :
257 {
258 momentum += tvdLi(momenta, distances, selfIsUpstream);
259 break;
260 }
261 case TvdApproach::hou :
262 {
263 momentum += tvdHou(momenta, distances, selfIsUpstream);
264 break;
265 }
266 default:
267 {
268 DUNE_THROW(ParameterException, "\nThis Tvd Approach is not implemented.\n");
269 break;
270 }
271 }
272 return momentum;
273 }
274
280 Scalar tvdUniform(const std::array<Scalar,3>& momenta,
281 const std::array<Scalar,3>& distances,
282 const bool selfIsUpstream) const
283 {
284 using std::isfinite;
285 const Scalar downstreamMomentum = momenta[0];
286 const Scalar upstreamMomentum = momenta[1];
287 const Scalar upUpstreamMomentum = momenta[2];
288 const Scalar ratio = (upstreamMomentum - upUpstreamMomentum) / (downstreamMomentum - upstreamMomentum);
289
290 // If the velocity field is uniform (like at the first newton step) we get a NaN
291 if(ratio > 0.0 && isfinite(ratio))
292 {
293 const Scalar secondOrderTerm = 0.5 * limiter_(ratio, 2.0) * (downstreamMomentum - upstreamMomentum);
294 return (upstreamMomentum + secondOrderTerm);
295 }
296 else
297 return upstreamMomentum;
298 }
299
307 Scalar tvdLi(const std::array<Scalar,3>& momenta,
308 const std::array<Scalar,3>& distances,
309 const bool selfIsUpstream) const
310 {
311 using std::isfinite;
312 const Scalar downstreamMomentum = momenta[0];
313 const Scalar upstreamMomentum = momenta[1];
314 const Scalar upUpstreamMomentum = momenta[2];
315 const Scalar upstreamToDownstreamDistance = distances[0];
316 const Scalar upUpstreamToUpstreamDistance = distances[1];
317
318 // selfIsUpstream is required to get the correct sign because upUpstreamToUpstreamDistance is always positive
319 const Scalar upUpstreamGradient = (upstreamMomentum - upUpstreamMomentum) / upUpstreamToUpstreamDistance * selfIsUpstream;
320
321 // Distance between the upUpstream node and the position where it should be if the grid were uniform.
322 const Scalar correctionDistance = upUpstreamToUpstreamDistance - upstreamToDownstreamDistance;
323 const Scalar reconstrutedUpUpstreamVelocity = upUpstreamMomentum + upUpstreamGradient * correctionDistance;
324 const Scalar ratio = (upstreamMomentum - reconstrutedUpUpstreamVelocity) / (downstreamMomentum - upstreamMomentum);
325
326 // If the velocity field is uniform (like at the first newton step) we get a NaN
327 if(ratio > 0.0 && isfinite(ratio))
328 {
329 const Scalar secondOrderTerm = 0.5 * limiter_(ratio, 2.0) * (downstreamMomentum - upstreamMomentum);
330 return (upstreamMomentum + secondOrderTerm);
331 }
332 else
333 return upstreamMomentum;
334 }
335
342 Scalar tvdHou(const std::array<Scalar,3>& momenta,
343 const std::array<Scalar,3>& distances,
344 const bool selfIsUpstream) const
345 {
346 using std::isfinite;
347 const Scalar downstreamMomentum = momenta[0];
348 const Scalar upstreamMomentum = momenta[1];
349 const Scalar upUpstreamMomentum = momenta[2];
350 const Scalar upstreamToDownstreamDistance = distances[0];
351 const Scalar upUpstreamToUpstreamDistance = distances[1];
352 const Scalar downstreamStaggeredCellSize = distances[2];
353 const Scalar ratio = (upstreamMomentum - upUpstreamMomentum) / (downstreamMomentum - upstreamMomentum)
354 * upstreamToDownstreamDistance / upUpstreamToUpstreamDistance;
355
356 // If the velocity field is uniform (like at the first newton step) we get a NaN
357 if(ratio > 0.0 && isfinite(ratio))
358 {
359 const Scalar upstreamStaggeredCellSize = 0.5 * (upstreamToDownstreamDistance + upUpstreamToUpstreamDistance);
360 const Scalar R = (upstreamStaggeredCellSize + downstreamStaggeredCellSize) / upstreamStaggeredCellSize;
361 const Scalar secondOrderTerm = limiter_(ratio, R) / R * (downstreamMomentum - upstreamMomentum);
362 return (upstreamMomentum + secondOrderTerm);
363 }
364 else
365 return upstreamMomentum;
366 }
367
373 static Scalar vanleer(const Scalar r, const Scalar R)
374 {
375 return R * r / (R - 1.0 + r);
376 }
377
381 static Scalar vanalbada(const Scalar r, const Scalar R)
382 {
383 return r * (r + 1.0) / (1.0 + r * r);
384 }
385
389 static Scalar minmod(const Scalar r, const Scalar R)
390 {
391 using std::min;
392 return min(r, 1.0);
393 }
394
400 static Scalar superbee(const Scalar r, const Scalar R)
401 {
402 using std::min;
403 using std::max;
404 return max(min(R * r, 1.0), min(r, R));
405 }
406
410 static Scalar umist(const Scalar r, const Scalar R)
411 {
412 using std::min;
413 return min({R * r, (r * (5.0 - R) + R - 1.0) / 4.0, (r * (R - 1.0) + 5.0 - R) / 4.0, R});
414 }
415
416 /*
417 * \brief Monotonized-Central limiter [Van Leer 1977]
418 */
419 static Scalar mclimiter(const Scalar r, const Scalar R)
420 {
421 using std::min;
422 return min({R * r, (r + 1.0) / 2.0, R});
423 }
424
428 static Scalar wahyd(const Scalar r, const Scalar R)
429 {
430 using std::min;
431 return r > 1 ? min((r + R * r * r) / (R + r * r), R)
432 : vanleer(r, R);
433 }
434
437 {
438 return tvdApproach_;
439 }
440
443 {
444 return differencingScheme_;
445 }
446
447private:
448 TvdApproach tvdApproach_;
449 DifferencingScheme differencingScheme_;
450 Scalar upwindWeight_;
451
452 std::function<Scalar(const Scalar, const Scalar)> limiter_;
453};
454
455} // end namespace Dumux
456
457#endif
Some exceptions thrown in DuMux
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
TvdApproach
Available Tvd approaches.
Definition: staggeredupwindmethods.hh:43
DifferencingScheme
Available differencing schemes.
Definition: staggeredupwindmethods.hh:52
bool hasParamInGroup(const std::string &paramGroup, const std::string &param)
Check whether a key exists in the parameter tree with a model group prefix.
Definition: parameters.hh:177
Definition: adapt.hh:29
Exception thrown if a run-time parameter is not specified correctly.
Definition: exceptions.hh:60
This file contains different higher order methods for approximating the velocity.
Definition: staggeredupwindmethods.hh:62
TvdApproach tvdApproachFromString(const std::string &tvd)
Convenience function to convert user input given as std::string to the corresponding enum class used ...
Definition: staggeredupwindmethods.hh:162
static Scalar vanalbada(const Scalar r, const Scalar R)
Van Albada flux limiter function [Van Albada et al. 1982].
Definition: staggeredupwindmethods.hh:381
StaggeredUpwindMethods(const std::string &paramGroup="")
Definition: staggeredupwindmethods.hh:64
Scalar upwind(const Scalar downstreamMomentum, const Scalar upstreamMomentum) const
Upwind Method.
Definition: staggeredupwindmethods.hh:233
const DifferencingScheme & differencingScheme() const
Returns the differencing scheme.
Definition: staggeredupwindmethods.hh:442
std::string tvdApproachToString(TvdApproach tvd)
return the name of the TVD approach
Definition: staggeredupwindmethods.hh:177
static Scalar umist(const Scalar r, const Scalar R)
UMIST flux limiter function [Lien and Leschziner 1993].
Definition: staggeredupwindmethods.hh:410
Scalar tvdUniform(const std::array< Scalar, 3 > &momenta, const std::array< Scalar, 3 > &distances, const bool selfIsUpstream) const
Tvd Scheme: Total Variation Diminishing.
Definition: staggeredupwindmethods.hh:280
Scalar tvdLi(const std::array< Scalar, 3 > &momenta, const std::array< Scalar, 3 > &distances, const bool selfIsUpstream) const
Tvd Scheme: Total Variation Diminishing.
Definition: staggeredupwindmethods.hh:307
static Scalar mclimiter(const Scalar r, const Scalar R)
Definition: staggeredupwindmethods.hh:419
DifferencingScheme differencingSchemeFromString(const std::string &differencingScheme)
Convenience function to convert user input given as std::string to the corresponding enum class used ...
Definition: staggeredupwindmethods.hh:192
Scalar tvd(const std::array< Scalar, 3 > &momenta, const std::array< Scalar, 3 > &distances, const bool selfIsUpstream, const TvdApproach tvdApproach) const
Tvd Scheme: Total Variation Diminishing.
Definition: staggeredupwindmethods.hh:243
Scalar tvdHou(const std::array< Scalar, 3 > &momenta, const std::array< Scalar, 3 > &distances, const bool selfIsUpstream) const
Tvd Scheme: Total Variation Diminishing.
Definition: staggeredupwindmethods.hh:342
static Scalar vanleer(const Scalar r, const Scalar R)
Van Leer flux limiter function [Van Leer 1974].
Definition: staggeredupwindmethods.hh:373
const TvdApproach & tvdApproach() const
Returns the Tvd approach.
Definition: staggeredupwindmethods.hh:436
std::string differencingSchemeToString(DifferencingScheme differencingScheme)
return the name of the Discretization Method
Definition: staggeredupwindmethods.hh:215
static Scalar superbee(const Scalar r, const Scalar R)
SUPERBEE flux limiter function [Roe 1985].
Definition: staggeredupwindmethods.hh:400
static Scalar minmod(const Scalar r, const Scalar R)
MinMod flux limiter function [Roe 1985].
Definition: staggeredupwindmethods.hh:389
static Scalar wahyd(const Scalar r, const Scalar R)
WAHYD Scheme [Hou, Simons, Hinkelmann 2007];.
Definition: staggeredupwindmethods.hh:428