Skip to content

Commit 34f8176

Browse files
jbreue16jbreue16
committed
additional slope limiters: Superbee, vanLeer, vanAlbada, relaxed vanAlbada:
1 parent e939294 commit 34f8176

File tree

1 file changed

+197
-53
lines changed

1 file changed

+197
-53
lines changed

‎src/libcadet/DGSubcellLimiterFV.hpp

Lines changed: 197 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -30,95 +30,222 @@ namespace cadet
3030
{
3131
/**
3232
* @brief Implements the subcell FV slope limited reconstruction
33+
* @detail Defines the 2nd order (linear) subcell FV reconstruction by implementing its limiter.
3334
*/
3435
// todo active types required in reconstruction?
3536
class SlopeLimiterFV {
3637
public:
38+
enum class SlopeLimiterID : int
39+
{
40+
NoLimiter1stOrder = 0, //!< 1st order (constant) reconstruction, i.e. no reconstruction
41+
Minmod = 1, //!< Minmod slope limiter
42+
Superbee = 2, //!< Superbee slope limiter
43+
vanLeer = 3, //!< van Leer slope limiter
44+
vanAlbada = 4, //!< van Alaba slope limiter
45+
RelaxedVanAlbada = 6, //!< relaxed van Alaba slope limiter
46+
NoLimiterForward = -1, //!< unlimited forward FD reconstruction
47+
NoLimiterBackward = -2, //!< unlimited backward FD reconstruction
48+
};
49+
3750
virtual ~SlopeLimiterFV() {}
38-
virtual double call(double slope0, double slope1) = 0;
39-
virtual active call(active slope0, active slope1) = 0;
51+
virtual double call(double fwdSlope, double bwdSlope) = 0;
52+
virtual active call(active fwdSlope, active bwdSlope) = 0;
53+
virtual SlopeLimiterID getID() = 0;
4054
};
4155

4256
/**
43-
* @brief Implements unlimited forward slope reconstruction
57+
* @brief 1st order (constant) reconstruction
4458
*/
45-
class NoLimiterForward : public SlopeLimiterFV {
59+
class NoLimiter1stOrder : public SlopeLimiterFV {
4660

4761
public:
48-
NoLimiterForward() {}
49-
double call(double slope0, double slope1) override {
50-
return slope1;
62+
NoLimiter1stOrder() {}
63+
double call(double fwdSlope, double bwdSlope) override {
64+
return 0.0;
5165
}
52-
active call(active slope0, active slope1) override {
53-
return slope1;
66+
active call(active fwdSlope, active bwdSlope) override {
67+
return 0.0;
5468
}
69+
SlopeLimiterID getID() override { return SlopeLimiterID::NoLimiter1stOrder; }
5570
};
5671

5772
/**
58-
* @brief Implements unlimited backward slope reconstruction
73+
* @brief MINMOD reconstruction slope limiter
5974
*/
60-
class NoLimiterBackward : public SlopeLimiterFV {
75+
class Minmod : public SlopeLimiterFV {
6176

6277
public:
63-
NoLimiterBackward() {}
64-
double call(double slope0, double slope1) override {
65-
return slope0;
78+
Minmod() {}
79+
double call(double fwdSlope, double bwdSlope) override {
80+
if (std::signbit(fwdSlope) == std::signbit(bwdSlope))
81+
return std::copysign(std::min(std::abs(fwdSlope), std::abs(bwdSlope)), fwdSlope);
82+
else
83+
return 0.0;
6684
}
67-
active call(active slope0, active slope1) override {
68-
return slope0;
85+
active call(active fwdSlope, active bwdSlope) override {
86+
87+
if (std::signbit(static_cast<double>(fwdSlope)) == std::signbit(static_cast<double>(bwdSlope)))
88+
return abs(static_cast<double>(fwdSlope)) < abs(static_cast<double>(bwdSlope)) ? fwdSlope : bwdSlope;
89+
else
90+
return active(0.0);
6991
}
92+
SlopeLimiterID getID() override { return SlopeLimiterID::Minmod; }
7093
};
71-
94+
7295
/**
73-
* @brief Disables reconstruction
96+
* @brief Superbee slope limiter for reconstruction
7497
*/
75-
class NoLimiter1stOrder : public SlopeLimiterFV {
98+
class Superbee : public SlopeLimiterFV {
7699

77100
public:
78-
NoLimiter1stOrder() {}
79-
double call(double slope0, double slope1) override {
80-
return 0.0;
101+
Superbee() {}
102+
103+
double maxmod(double fwdSlope, double bwdSlope) {
104+
if (std::signbit(fwdSlope) == std::signbit(bwdSlope))
105+
return std::copysign(std::max(std::abs(fwdSlope), std::abs(bwdSlope)), fwdSlope);
106+
else
107+
return 0.0;
81108
}
82-
active call(active slope0, active slope1) override {
83-
return 0.0;
109+
active maxmod(active fwdSlope, active bwdSlope) {
110+
111+
if (std::signbit(static_cast<double>(fwdSlope)) == std::signbit(static_cast<double>(bwdSlope)))
112+
return abs(static_cast<double>(fwdSlope)) < abs(static_cast<double>(bwdSlope)) ? bwdSlope : fwdSlope;
113+
else
114+
return active(0.0);
115+
}
116+
117+
double call(double fwdSlope, double bwdSlope) override {
118+
return maxmod(minmod.call(2 * fwdSlope, bwdSlope), minmod.call(fwdSlope, 2 * bwdSlope));
119+
}
120+
active call(active fwdSlope, active bwdSlope) override {
121+
return maxmod(minmod.call(2 * fwdSlope, bwdSlope), minmod.call(fwdSlope, 2 * bwdSlope));
84122
}
123+
SlopeLimiterID getID() override { return SlopeLimiterID::Superbee; }
124+
125+
private:
126+
Minmod minmod;
85127
};
86128

87-
///**
88-
// * @brief Implements unlimited central slope reconstruction
89-
// */
90-
//class NoLimiterCentral : public SlopeLimiterFV {
91-
92-
//public:
93-
// NoLimiterCentral() {}
94-
// double call(double slope0, double slope1) override {
95-
// return slope0;
96-
// }
97-
// active call(active slope0, active slope1) override {
98-
// return slope0;
99-
// }
100-
//};
129+
/**
130+
* @brief van Leer slope limiter for reconstruction
131+
*/
132+
class vanLeer : public SlopeLimiterFV {
133+
134+
public:
135+
vanLeer() {}
101136

102-
class Minmod : public SlopeLimiterFV {
137+
double call(double fwdSlope, double bwdSlope) override {
138+
if (std::signbit(fwdSlope) == std::signbit(bwdSlope))
139+
{
140+
if (std::abs(fwdSlope) < eps && std::abs(bwdSlope) < eps) // catch zero-division problem
141+
return 0.0;
142+
else
143+
return (2.0 * fwdSlope * bwdSlope) / (fwdSlope + bwdSlope);
144+
}
145+
else
146+
return 0.0;
147+
}
148+
active call(active fwdSlope, active bwdSlope) override {
149+
if (std::signbit(static_cast<double>(fwdSlope)) == std::signbit(static_cast<double>(bwdSlope)))
150+
{
151+
if (abs(fwdSlope) < eps && abs(bwdSlope) < eps) // catch zero-division problem
152+
return 0.0;
153+
else
154+
return (2.0 * fwdSlope * bwdSlope) / (fwdSlope + bwdSlope);
155+
}
156+
else
157+
return 0.0;
158+
}
159+
SlopeLimiterID getID() override { return SlopeLimiterID::vanLeer; }
160+
161+
private:
162+
const double eps = 1e-8;
163+
};
164+
165+
/**
166+
* @brief van Albada slope limiter for reconstruction
167+
*/
168+
class vanAlbada : public SlopeLimiterFV {
103169

104170
public:
105-
Minmod() {}
106-
double call(double slope0, double slope1) override {
107-
if (std::signbit(slope0) == std::signbit(slope1))
108-
return std::copysign(std::min(std::abs(slope0), std::abs(slope1)), slope0);
171+
vanAlbada() {}
172+
173+
double call(double fwdSlope, double bwdSlope) override {
174+
if (std::signbit(fwdSlope) == std::signbit(bwdSlope))
175+
return ((fwdSlope * fwdSlope + eps) * bwdSlope + ((bwdSlope * bwdSlope + eps) * fwdSlope)) / (fwdSlope * fwdSlope + bwdSlope * bwdSlope + 2.0 * eps);
109176
else
110177
return 0.0;
111178
}
112-
active call(active slope0, active slope1) override {
179+
active call(active fwdSlope, active bwdSlope) override {
180+
if (std::signbit(static_cast<double>(fwdSlope)) == std::signbit(static_cast<double>(bwdSlope)))
181+
return ((fwdSlope * fwdSlope + eps) * bwdSlope + ((bwdSlope * bwdSlope + eps) * fwdSlope)) / (fwdSlope * fwdSlope + bwdSlope * bwdSlope + 2.0 * eps);
182+
else
183+
return 0.0;
184+
}
185+
SlopeLimiterID getID() override { return SlopeLimiterID::vanAlbada; }
113186

114-
if (std::signbit(static_cast<double>(slope0)) == std::signbit(static_cast<double>(slope1)))
115-
return abs(static_cast<double>(slope0)) < abs(static_cast<double>(slope1)) ? slope0 : slope1;
187+
private:
188+
const double eps = 1e-6; // constant chosen according to https://doi.org/10.2514/3.9465
189+
};
190+
191+
/**
192+
* @brief relaxed van Albada slope limiter for reconstruction
193+
* @detail relaxed w.r.t. required arithmetic operations, according to https://doi.org/10.1007/978-981-15-9011-5_5
194+
*/
195+
class RelaxedVanAlbada : public SlopeLimiterFV {
196+
197+
public:
198+
RelaxedVanAlbada() {}
199+
200+
double call(double fwdSlope, double bwdSlope) override {
201+
if (std::signbit(fwdSlope) == std::signbit(bwdSlope))
202+
return fwdSlope * (2.0 * fwdSlope * bwdSlope + eps) / (fwdSlope * fwdSlope + bwdSlope * bwdSlope + eps);
116203
else
117-
return active(0.0);
204+
return 0.0;
118205
}
206+
active call(active fwdSlope, active bwdSlope) override {
207+
if (std::signbit(static_cast<double>(fwdSlope)) == std::signbit(static_cast<double>(bwdSlope)))
208+
return fwdSlope * (2.0 * fwdSlope * bwdSlope + eps) / (fwdSlope * fwdSlope + bwdSlope * bwdSlope + eps);
209+
else
210+
return 0.0;
211+
}
212+
SlopeLimiterID getID() override { return SlopeLimiterID::RelaxedVanAlbada; }
213+
214+
private:
215+
const double eps = 1e-6; // constant chosen according to https://doi.org/10.2514/3.9465
119216
};
120217

121-
// todo, note: if more slope limiters are added, note that they have to be symmetrical in the current implementation (i.e. in reconstruction function reconstructedUpwindValue)
218+
/**
219+
* @brief Implements unlimited forward slope for reconstruction
220+
*/
221+
class NoLimiterForward : public SlopeLimiterFV {
222+
223+
public:
224+
NoLimiterForward() {}
225+
double call(double fwdSlope, double bwdSlope) override {
226+
return bwdSlope;
227+
}
228+
active call(active fwdSlope, active bwdSlope) override {
229+
return bwdSlope;
230+
}
231+
SlopeLimiterID getID() override { return SlopeLimiterID::NoLimiterForward; }
232+
};
233+
234+
/**
235+
* @brief Implements unlimited backward slope for reconstruction
236+
*/
237+
class NoLimiterBackward : public SlopeLimiterFV {
238+
239+
public:
240+
NoLimiterBackward() {}
241+
double call(double fwdSlope, double bwdSlope) override {
242+
return fwdSlope;
243+
}
244+
active call(active fwdSlope, active bwdSlope) override {
245+
return fwdSlope;
246+
}
247+
SlopeLimiterID getID() override { return SlopeLimiterID::NoLimiterBackward; }
248+
};
122249

123250
/**
124251
* @brief Implements the subcell FV limiting scheme for convection
@@ -169,16 +296,22 @@ namespace cadet
169296
_subcellGrid[subcell] = _subcellGrid[subcell - 1] + _LGLweights[subcell - 1];
170297

171298
_FVorder = FVorder;
299+
_slope_limiter = std::make_unique<NoLimiter1stOrder>();
300+
172301
if (_FVorder == 2) {
173302
if (limiter == "MINMOD")
174303
_slope_limiter = std::make_unique<Minmod>();
175304
else if (limiter == "UNLIMITED_FORWARD_RECONSTRUCTION")
176305
_slope_limiter = std::make_unique<NoLimiterForward>();
177-
else if (limiter == "UNLIMITED_BACKWARD_RECONSTRUCTION")
178-
_slope_limiter = std::make_unique<NoLimiterBackward>();
179-
else if (limiter == "NONE")
180-
_slope_limiter = std::make_unique<NoLimiter1stOrder>();
181-
else
306+
else if (limiter == "SUPERBEE")
307+
_slope_limiter = std::make_unique<Superbee>();
308+
else if (limiter == "VANLEER")
309+
_slope_limiter = std::make_unique<vanLeer>();
310+
else if (limiter == "VANALBADA")
311+
_slope_limiter = std::make_unique<vanAlbada>();
312+
else if (limiter == "RELAXED_VANALBADA")
313+
_slope_limiter = std::make_unique<RelaxedVanAlbada>();
314+
else if (limiter != "NONE")
182315
throw InvalidParameterException("Subcell FV slope limiter " + limiter + " unknown.");
183316

184317
switch (boundaryTreatment)
@@ -200,6 +333,17 @@ namespace cadet
200333
throw InvalidParameterException("Subcell FV order must be 1 or 2, but was specified as " + std::to_string(_FVorder));
201334
}
202335

336+
void notifyDiscontinuousSectionTransition(int direction) {
337+
338+
// Forward FD slope reconstruction depends on flow direction // todo: Koren limiter which is also non-symmetrical
339+
if (_slope_limiter->getID() == SlopeLimiterFV::SlopeLimiterID::NoLimiterForward && direction == -1)
340+
_slope_limiter = std::make_unique<NoLimiterBackward>();
341+
342+
else if (_slope_limiter->getID() == SlopeLimiterFV::SlopeLimiterID::NoLimiterBackward && direction == 1)
343+
_slope_limiter = std::make_unique<NoLimiterForward>();
344+
345+
}
346+
203347
/**
204348
* @brief Implements FV (limited) reconstruction
205349
* @details TODO

0 commit comments

Comments
 (0)