|
| 1 | +// ============================================================================= |
| 2 | +// CADET |
| 3 | +// |
| 4 | +// Copyright © 2008-2022: The CADET Authors |
| 5 | +// Please see the AUTHORS and CONTRIBUTORS file. |
| 6 | +// |
| 7 | +// All rights reserved. This program and the accompanying materials |
| 8 | +// are made available under the terms of the GNU Public License v3.0 (or, at |
| 9 | +// your option, any later version) which accompanies this distribution, and |
| 10 | +// is available at http://www.gnu.org/licenses/gpl.html |
| 11 | +// ============================================================================= |
| 12 | + |
| 13 | +/** |
| 14 | + * @file |
| 15 | + * Implements the WENO method |
| 16 | + */ |
| 17 | + |
| 18 | +#ifndef LIBCADET_WENODG_HPP_ |
| 19 | +#define LIBCADET_WENODG_HPP_ |
| 20 | + |
| 21 | +#include "AutoDiff.hpp" |
| 22 | +#include "MathUtil.hpp" |
| 23 | +#include "Memory.hpp" |
| 24 | +#include "common/CompilerSpecific.hpp" |
| 25 | +#include "cadet/Exceptions.hpp" |
| 26 | + |
| 27 | +#include <algorithm> |
| 28 | + |
| 29 | +namespace cadet |
| 30 | +{ |
| 31 | + |
| 32 | + //template<typename StateType> // todo |
| 33 | + class WENOLimiter { |
| 34 | + public: |
| 35 | + virtual ~WENOLimiter() {} |
| 36 | + virtual double call(double u0, double u1, double u2) = 0; |
| 37 | + virtual double call(active u0, active u1, active u2) = 0; |
| 38 | + }; |
| 39 | + |
| 40 | + class FullLimiter : public WENOLimiter { |
| 41 | + |
| 42 | + public: |
| 43 | + FullLimiter() {} |
| 44 | + double call(double u0, double u1, double u2) override { |
| 45 | + return 0.0; |
| 46 | + } |
| 47 | + double call(active u0, active u1, active u2) override { |
| 48 | + return 0.0; |
| 49 | + } |
| 50 | + }; |
| 51 | + |
| 52 | + class MinmodWENO : public WENOLimiter { |
| 53 | + |
| 54 | + public: |
| 55 | + MinmodWENO() {} |
| 56 | + double call(double u0, double u1, double u2) override { |
| 57 | + if (std::signbit(u0) == std::signbit(u1) && std::signbit(u0) == std::signbit(u2)) |
| 58 | + return std::copysign(std::min(std::abs(u0), std::min(std::abs(u1), std::abs(u2))), u0); |
| 59 | + else |
| 60 | + return 0.0; |
| 61 | + } |
| 62 | + double call(active u0, active u1, active u2) override { |
| 63 | + if (std::signbit(u0.getValue()) == std::signbit(u1.getValue()) && std::signbit(u0.getValue()) == std::signbit(u2.getValue())) { |
| 64 | + if (std::signbit(u0.getValue())) // i.e. negative sign |
| 65 | + return -static_cast<double>(u0); |
| 66 | + else |
| 67 | + return static_cast<double>(u0); |
| 68 | + } |
| 69 | + else |
| 70 | + return 0.0; |
| 71 | + } |
| 72 | + }; |
| 73 | + |
| 74 | + class TVBMinmodWENO : public WENOLimiter { |
| 75 | + |
| 76 | + public: |
| 77 | + TVBMinmodWENO() {} |
| 78 | + double call(double u0, double u1, double u2) override { |
| 79 | + if (std::abs(u0) <= M * h * h) |
| 80 | + return u0; |
| 81 | + else |
| 82 | + return minmod.call(u0, u1, u2); |
| 83 | + } |
| 84 | + double call(active u0, active u1, active u2) override { |
| 85 | + if (std::abs(u0.getValue()) <= M * h * h) |
| 86 | + return static_cast<double>(u0); |
| 87 | + else |
| 88 | + return minmod.call(u0, u1, u2); |
| 89 | + } |
| 90 | + private: |
| 91 | + MinmodWENO minmod; |
| 92 | + active h; |
| 93 | + active M; |
| 94 | + }; |
| 95 | + |
| 96 | + /** |
| 97 | + * @brief Implements the WENO scheme for convection |
| 98 | + * @details This scheme is based on upwind stencils and provides WENO methods 1-1, 2-3, and 3-5. |
| 99 | + * In general, WENO achieves order \f$ r \f$ using a stencil with \f$ (2r-1) \f$ points |
| 100 | + * that is subdivided into \f$ r \f$ smaller stencils having \f$ r \f$ points each. |
| 101 | + * WENO combines all substencils with an estimate of the smoothness of the solution (also obtained from the |
| 102 | + * substencils) in order to achieve a non-oscillatory high order reconstruction of the face values given |
| 103 | + * volume averages (cell boundary fluxes in finite volume schemes). |
| 104 | + * For details see \cite Liu1994 and \cite Jiang1996. |
| 105 | + */ |
| 106 | + class WenoDG |
| 107 | + { |
| 108 | + public: |
| 109 | + |
| 110 | + /** |
| 111 | + * @brief Boundary treatment method determines how the reconstruction handles insufficient available elements (i.e., less elements available than stencil size) |
| 112 | + */ |
| 113 | + enum class BoundaryTreatment : int |
| 114 | + { |
| 115 | + ReduceOrder = 0, //!< Reduce the order of the WENO method such that the stencil is small enough |
| 116 | + ZeroWeights = 1, //!< Set weights of WENO method to 0 for unavailable elements |
| 117 | + ZeroWeightsForPnotZero = 2, //!< Set weights of WENO method to 0 for unavailable elements, except for the first cell (order is reduced to 1) |
| 118 | + LargeGhostNodes = 3, |
| 119 | + }; |
| 120 | + |
| 121 | + /** |
| 122 | + * @brief Creates the WENO scheme |
| 123 | + */ |
| 124 | + //WenoDG() : _order(maxOrder()), _boundaryTreatment(BoundaryTreatment::ReduceOrder), _intermediateValues(3 * maxOrder() * sizeof(active)) { } |
| 125 | + |
| 126 | + ~WenoDG() CADET_NOEXCEPT |
| 127 | + { |
| 128 | + delete[] _weno_gamma; |
| 129 | + delete[] _w; |
| 130 | + delete[] _wHat; |
| 131 | + delete[] _beta; |
| 132 | + delete[] _troubled_cells; |
| 133 | + } |
| 134 | + |
| 135 | + void init(const double eps, const double r, const double gamma, const unsigned int nCells, const unsigned int nComp, const bool trackTroubledCells = false) { |
| 136 | + |
| 137 | + _weno_gamma = new double[3]{ (1 - gamma) / 2.0, gamma, (1 - gamma) / 2.0 }; |
| 138 | + _weno_r = r; |
| 139 | + _weno_eps = eps; |
| 140 | + if (trackTroubledCells) |
| 141 | + _troubled_cells = new double[nCells * nComp]; |
| 142 | + else |
| 143 | + _troubled_cells = nullptr; |
| 144 | + } |
| 145 | + |
| 146 | + inline double* troubledCells() const CADET_NOEXCEPT { return &_troubled_cells[0]; } |
| 147 | + |
| 148 | + //private: |
| 149 | + |
| 150 | + // indicator for troubled cells |
| 151 | + std::unique_ptr<WENOLimiter> weno_limiter; |
| 152 | + |
| 153 | + // WENO parameters |
| 154 | + double* _weno_gamma; |
| 155 | + double _weno_r = 2.0; |
| 156 | + double _weno_eps = 1e-6; |
| 157 | + // weights, smoothness indicator, integral values. |
| 158 | + active* _wHat = new active[3]{ 0.0, 0.0, 0.0 }; |
| 159 | + active* _w = new active[3]{ 0.0, 0.0, 0.0 }; |
| 160 | + active _wSum = 0.0; |
| 161 | + active* _beta = new active[3]{ 0.0, 0.0, 0.0 }; |
| 162 | + active _pAvg0 = 0.0; |
| 163 | + active _pAvg1 = 0.0; |
| 164 | + active _pAvg2 = 0.0; |
| 165 | + // limiter inputs |
| 166 | + active _avgUdeltaPlus = 0.0; |
| 167 | + active _avgUdeltaMinus = 0.0; |
| 168 | + active _uTilde = 0.0; |
| 169 | + active _u2Tilde = 0.0; |
| 170 | + // mark troubled cells (smoothness indicator) |
| 171 | + double* _troubled_cells; |
| 172 | + |
| 173 | + }; |
| 174 | + |
| 175 | +} // namespace cadet |
| 176 | + |
| 177 | +#endif // LIBCADET_WENODG_HPP_ |
0 commit comments