diff --git a/src/libcadet/AutoDiff.hpp b/src/libcadet/AutoDiff.hpp index 14a15369f..78816349b 100644 --- a/src/libcadet/AutoDiff.hpp +++ b/src/libcadet/AutoDiff.hpp @@ -25,7 +25,7 @@ #if defined(ACTIVE_SFAD) || defined(ACTIVE_SETFAD) - #define SFAD_DEFAULT_DIR 80 + #define SFAD_DEFAULT_DIR 100 #if defined(ACTIVE_SFAD) #include "sfad.hpp" diff --git a/src/libcadet/DGSubcellLimiterFV.hpp b/src/libcadet/DGSubcellLimiterFV.hpp new file mode 100644 index 000000000..5b6796eab --- /dev/null +++ b/src/libcadet/DGSubcellLimiterFV.hpp @@ -0,0 +1,501 @@ +// ============================================================================= +// CADET +// +// Copyright � 2008-2022: The CADET Authors +// Please see the AUTHORS and CONTRIBUTORS file. +// +// All rights reserved. This program and the accompanying materials +// are made available under the terms of the GNU Public License v3.0 (or, at +// your option, any later version) which accompanies this distribution, and +// is available at http://www.gnu.org/licenses/gpl.html +// ============================================================================= + +/** + * @file + * Implements the subcell FV limiting method for collocation DG + */ + +#ifndef LIBCADET_DGSubcellLimiterFV_HPP_ +#define LIBCADET_DGSubcellLimiterFV_HPP_ + +#include "AutoDiff.hpp" +#include "MathUtil.hpp" +#include "Memory.hpp" +#include "common/CompilerSpecific.hpp" +#include "cadet/Exceptions.hpp" + +#include + +namespace cadet +{ + /** + * @brief Implements the subcell FV slope limited reconstruction + * @detail Defines the 2nd order (linear) subcell FV reconstruction by implementing its limiter. + */ + // todo active types required in reconstruction? + class SlopeLimiterFV { + public: + enum class SlopeLimiterID : int + { + NoLimiter1stOrder = 0, //!< 1st order (constant) reconstruction, i.e. no reconstruction + Minmod = 1, //!< Minmod slope limiter + Superbee = 2, //!< Superbee slope limiter + vanLeer = 3, //!< van Leer slope limiter + vanAlbada = 4, //!< van Alaba slope limiter + //Koren = 5, //!< Koren slope limiter // todo: can we use a 3rd order approximation, i.e. can we accomodate the MUSCL scheme/k-interpolation(vanLeer) into our setting? + RelaxedVanAlbada = 6, //!< relaxed van Alaba slope limiter + //MonotonizedCentral = 7, //!< Monotonized central slope limiter // interface change required to compute central slope + NoLimiterForward = -1, //!< unlimited forward FD reconstruction + NoLimiterBackward = -2, //!< unlimited backward FD reconstruction + }; + + virtual ~SlopeLimiterFV() {} + virtual double call(double fwdSlope, double bwdSlope) = 0; + virtual active call(active fwdSlope, active bwdSlope) = 0; + virtual SlopeLimiterID getID() = 0; + }; + + /** + * @brief 1st order (constant) reconstruction + */ + class NoLimiter1stOrder : public SlopeLimiterFV { + + public: + NoLimiter1stOrder() {} + double call(double fwdSlope, double bwdSlope) override { + return 0.0; + } + active call(active fwdSlope, active bwdSlope) override { + return 0.0; + } + SlopeLimiterID getID() override { return SlopeLimiterID::NoLimiter1stOrder; } + }; + + /** + * @brief MINMOD reconstruction slope limiter + */ + class Minmod : public SlopeLimiterFV { + + public: + Minmod() {} + double call(double fwdSlope, double bwdSlope) override { + if (std::signbit(fwdSlope) == std::signbit(bwdSlope)) + return std::copysign(std::min(std::abs(fwdSlope), std::abs(bwdSlope)), fwdSlope); + else + return 0.0; + } + active call(active fwdSlope, active bwdSlope) override { + + if (std::signbit(static_cast(fwdSlope)) == std::signbit(static_cast(bwdSlope))) + return abs(static_cast(fwdSlope)) < abs(static_cast(bwdSlope)) ? fwdSlope : bwdSlope; + else + return active(0.0); + } + SlopeLimiterID getID() override { return SlopeLimiterID::Minmod; } + }; + + /** + * @brief Superbee slope limiter for reconstruction + */ + class Superbee : public SlopeLimiterFV { + + public: + Superbee() {} + + double maxmod(double fwdSlope, double bwdSlope) { + if (std::signbit(fwdSlope) == std::signbit(bwdSlope)) + return std::copysign(std::max(std::abs(fwdSlope), std::abs(bwdSlope)), fwdSlope); + else + return 0.0; + } + active maxmod(active fwdSlope, active bwdSlope) { + + if (std::signbit(static_cast(fwdSlope)) == std::signbit(static_cast(bwdSlope))) + return abs(static_cast(fwdSlope)) < abs(static_cast(bwdSlope)) ? bwdSlope : fwdSlope; + else + return active(0.0); + } + + double call(double fwdSlope, double bwdSlope) override { + return maxmod(minmod.call(2 * fwdSlope, bwdSlope), minmod.call(fwdSlope, 2 * bwdSlope)); + } + active call(active fwdSlope, active bwdSlope) override { + return maxmod(minmod.call(2 * fwdSlope, bwdSlope), minmod.call(fwdSlope, 2 * bwdSlope)); + } + SlopeLimiterID getID() override { return SlopeLimiterID::Superbee; } + + private: + Minmod minmod; + }; + + /** + * @brief van Leer slope limiter for reconstruction + */ + class vanLeer : public SlopeLimiterFV { + + public: + vanLeer() {} + + double call(double fwdSlope, double bwdSlope) override { + if (std::signbit(fwdSlope) == std::signbit(bwdSlope)) + { + if (std::abs(fwdSlope) < eps && std::abs(bwdSlope) < eps) // catch zero-division problem + return 0.0; + else + return (2.0 * fwdSlope * bwdSlope) / (fwdSlope + bwdSlope); + } + else + return 0.0; + } + active call(active fwdSlope, active bwdSlope) override { + if (std::signbit(static_cast(fwdSlope)) == std::signbit(static_cast(bwdSlope))) + { + if (abs(fwdSlope) < eps && abs(bwdSlope) < eps) // catch zero-division problem + return 0.0; + else + return (2.0 * fwdSlope * bwdSlope) / (fwdSlope + bwdSlope); + } + else + return 0.0; + } + SlopeLimiterID getID() override { return SlopeLimiterID::vanLeer; } + + private: + const double eps = 1e-8; + }; + + /** + * @brief van Albada slope limiter for reconstruction + */ + class vanAlbada : public SlopeLimiterFV { + + public: + vanAlbada() {} + + double call(double fwdSlope, double bwdSlope) override { + if (std::signbit(fwdSlope) == std::signbit(bwdSlope)) + return ((fwdSlope * fwdSlope + eps) * bwdSlope + ((bwdSlope * bwdSlope + eps) * fwdSlope)) / (fwdSlope * fwdSlope + bwdSlope * bwdSlope + 2.0 * eps); + else + return 0.0; + } + active call(active fwdSlope, active bwdSlope) override { + if (std::signbit(static_cast(fwdSlope)) == std::signbit(static_cast(bwdSlope))) + return ((fwdSlope * fwdSlope + eps) * bwdSlope + ((bwdSlope * bwdSlope + eps) * fwdSlope)) / (fwdSlope * fwdSlope + bwdSlope * bwdSlope + 2.0 * eps); + else + return 0.0; + } + SlopeLimiterID getID() override { return SlopeLimiterID::vanAlbada; } + + private: + const double eps = 1e-6; // constant chosen according to https://doi.org/10.2514/3.9465 + }; + + /** + * @brief relaxed van Albada slope limiter for reconstruction + * @detail relaxed w.r.t. required arithmetic operations, according to https://doi.org/10.1007/978-981-15-9011-5_5 + */ + class RelaxedVanAlbada : public SlopeLimiterFV { + + public: + RelaxedVanAlbada() {} + + double call(double fwdSlope, double bwdSlope) override { + if (std::signbit(fwdSlope) == std::signbit(bwdSlope)) + return fwdSlope * (2.0 * fwdSlope * bwdSlope + eps) / (fwdSlope * fwdSlope + bwdSlope * bwdSlope + eps); + else + return 0.0; + } + active call(active fwdSlope, active bwdSlope) override { + if (std::signbit(static_cast(fwdSlope)) == std::signbit(static_cast(bwdSlope))) + return fwdSlope * (2.0 * fwdSlope * bwdSlope + eps) / (fwdSlope * fwdSlope + bwdSlope * bwdSlope + eps); + else + return 0.0; + } + SlopeLimiterID getID() override { return SlopeLimiterID::RelaxedVanAlbada; } + + private: + const double eps = 1e-6; // constant chosen according to https://doi.org/10.2514/3.9465 + }; + + ///** + //* @brief Monotonized central slope limiter for reconstruction + //* @todo change interface to compute central slope? + //*/ + //class MonotonizedCentral : public SlopeLimiterFV { + + //public: + // MonotonizedCentral() {} + + // double minmod3(double fwdSlope, double cSlope, double bwdSlope) { + // if (std::signbit(fwdSlope) == std::signbit(cSlope) && std::signbit(fwdSlope) == std::signbit(bwdSlope)) + // return std::copysign(std::max(std::abs(fwdSlope), std::abs(cSlope), std::abs(bwdSlope)), fwdSlope); + // else + // return 0.0; + // } + // active minmod3(active fwdSlope, active cSlope, active bwdSlope) { + + // if (std::signbit(static_cast(fwdSlope)) == std::signbit(static_cast(cSlope)) && std::signbit(static_cast(fwdSlope)) == std::signbit(static_cast(bwdSlope))) + // { + // if (abs(static_cast(fwdSlope)) < abs(static_cast(bwdSlope))) + // { + // if (abs(static_cast(fwdSlope)) < abs(static_cast(cSlope))) + // { + // return fwdSlope; + // } + // else return cSlope; + // } + // else + // { + // if (abs(static_cast(bwdSlope)) < abs(static_cast(cSlope))) + // { + // return bwdSlope; + // } + // else return cSlope; + // } + // } + // else + // return active(0.0); + // } + + // double call(double fwdSlope, double bwdSlope) override { + // return minmod3(2.0 * fwdSlope, cSlope, 2.0 * bwdSlope); + // } + // active call(active fwdSlope, active bwdSlope) override { + // return minmod3(2.0 * fwdSlope, cSlope, 2.0 * bwdSlope); + // } + + // SlopeLimiterID getID() override { return SlopeLimiterID::MonotonizedCentral; } + //}; + + /** + * @brief Implements unlimited forward slope for reconstruction + */ + class NoLimiterForward : public SlopeLimiterFV { + + public: + NoLimiterForward() {} + double call(double fwdSlope, double bwdSlope) override { + return bwdSlope; + } + active call(active fwdSlope, active bwdSlope) override { + return bwdSlope; + } + SlopeLimiterID getID() override { return SlopeLimiterID::NoLimiterForward; } + }; + + /** + * @brief Implements unlimited backward slope for reconstruction + */ + class NoLimiterBackward : public SlopeLimiterFV { + + public: + NoLimiterBackward() {} + double call(double fwdSlope, double bwdSlope) override { + return fwdSlope; + } + active call(active fwdSlope, active bwdSlope) override { + return fwdSlope; + } + SlopeLimiterID getID() override { return SlopeLimiterID::NoLimiterBackward; } + }; + + ///** + // * @brief Implements unlimited central slope reconstruction + // */ + //class NoLimiterCentral : public SlopeLimiterFV { + + //public: + // NoLimiterCentral() {} + // double call(double fwdSlope, double bwdSlope) override { + // return fwdSlope; + // } + // active call(active fwdSlope, active bwdSlope) override { + // return fwdSlope; + // } + // SlopeLimiterID getID() override { return SlopeLimiterID::NoLimiterCentral; } + //}; + + // 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) + + /** + * @brief Implements the subcell FV limiting scheme for convection + * @details //@TODO. + */ + class DGSubcellLimiterFV + { + public: + + /** + * @brief Boundary treatment method determines how the reconstruction handles DG element boundaries. + */ + enum class BoundaryTreatment : int + { + LimiterSlope = 0, //!< Slope limiter reconstruction using neighbour interface value. + CentralSlope = 1, //!< Central slope reconstruction with interior information. + Constant = 2 //!< Constant reconstruction of boundary subcells. + }; + + /** + * @brief Creates the subcell Finite Volume scheme + */ + DGSubcellLimiterFV() : _LGLweights(nullptr), _LGLnodes(nullptr), _subcellGrid(nullptr) { } + + ~DGSubcellLimiterFV() CADET_NOEXCEPT + { + delete[] _LGLweights; + delete[] _LGLnodes; + delete[] _subcellGrid; + } + + void init(std::string limiter, const int FVorder, const int boundaryTreatment, const unsigned int nNodes, double* LGLnodes, double* invWeights, const unsigned int nCells, const unsigned int nComp) { + + _nNodes = nNodes; + _polyDeg = nNodes - 1; + _nComp = nComp; + + _LGLweights = new double[nNodes]; + for (int node = 0; node < nNodes; node++) + _LGLweights[node] = 1.0 / invWeights[node]; + + _LGLnodes = new double[nNodes]; + std::copy(LGLnodes, LGLnodes + nNodes, _LGLnodes); + + _subcellGrid = new double[nNodes + 1]; + _subcellGrid[0] = -1.0; + for (int subcell = 1; subcell < nNodes + 1; subcell++) + _subcellGrid[subcell] = _subcellGrid[subcell - 1] + _LGLweights[subcell - 1]; + + _FVorder = FVorder; + _slope_limiter = std::make_unique(); + + if (_FVorder == 2) { + if (limiter == "MINMOD") + _slope_limiter = std::make_unique(); + else if (limiter == "UNLIMITED_FORWARD_RECONSTRUCTION") + _slope_limiter = std::make_unique(); + else if (limiter == "SUPERBEE") + _slope_limiter = std::make_unique(); + else if (limiter == "VANLEER") + _slope_limiter = std::make_unique(); + else if (limiter == "VANALBADA") + _slope_limiter = std::make_unique(); + else if (limiter == "RELAXED_VANALBADA") + _slope_limiter = std::make_unique(); + else if (limiter != "NONE") + throw InvalidParameterException("Subcell FV slope limiter " + limiter + " unknown."); + + switch (boundaryTreatment) + { + case static_cast::type>(BoundaryTreatment::LimiterSlope): + _FVboundaryTreatment = BoundaryTreatment::LimiterSlope; + break; + case static_cast::type>(BoundaryTreatment::CentralSlope): + _FVboundaryTreatment = BoundaryTreatment::CentralSlope; + break; + case static_cast::type>(BoundaryTreatment::Constant): + _FVboundaryTreatment = BoundaryTreatment::Constant; + break; + default: + throw InvalidParameterException("Unknown subcell FV boundary treatment."); + } + } + else if (_FVorder != 1) + throw InvalidParameterException("Subcell FV order must be 1 or 2, but was specified as " + std::to_string(_FVorder)); + } + + void notifyDiscontinuousSectionTransition(int direction) { + + if (_slope_limiter) { // only when subcell limiter is initialized, i.e. used + // Forward FD slope reconstruction depends on flow direction // todo: Koren limiter which is also non-symmetrical + if (_slope_limiter->getID() == SlopeLimiterFV::SlopeLimiterID::NoLimiterForward && direction == -1) + _slope_limiter = std::make_unique(); + + else if (_slope_limiter->getID() == SlopeLimiterFV::SlopeLimiterID::NoLimiterBackward && direction == 1) + _slope_limiter = std::make_unique(); + } + } + + /** + * @brief Implements FV (limited) reconstruction + * @details TODO + * @param [in] leftState left neighbour state + * @param [in] centerState current state + * @param [in] rightState right neighbour state + * @param [in] subcellIdx current (center state) subcell index + * @param [in] rightInterface specifies whther right or left interface reconstruction value is returned + * @return @c reconstructed FV subcell value at left or right interface + */ + // forward backward flow unterscheidung. Check slope computation. + template + StateType reconstructedInterfaceValue(const StateType leftState, const StateType centerState, const StateType rightState, const int subcellIdx, bool rightInterface) { + + if (_FVorder == 1) // No reconstruction + return centerState; + else // Order = 2 -> reconstruction + { + // TODO what happens here with non-symmetrical slope limiters? + StateType slope; + + // Boundary cell reconstruction + if (subcellIdx == 0) + { + switch (_FVboundaryTreatment) + { + default: + case BoundaryTreatment::LimiterSlope: + slope = _slope_limiter->call((rightState - centerState) / (_LGLnodes[1] - _LGLnodes[0]), (rightState - leftState) / (_LGLnodes[1] - _LGLnodes[0])); + break; + + case BoundaryTreatment::CentralSlope: + slope = (rightState - centerState) / (_LGLnodes[1] - _LGLnodes[0]); + break; + + case BoundaryTreatment::Constant: + return centerState; + } + } + else if (subcellIdx == _nNodes - 1) + { + switch (_FVboundaryTreatment) + { + default: + case BoundaryTreatment::LimiterSlope: + slope = _slope_limiter->call((centerState - leftState) / (_LGLnodes[_nNodes - 1] - _LGLnodes[_nNodes - 2]), (rightState - leftState) / (_LGLnodes[_nNodes - 1] - _LGLnodes[_nNodes - 2])); + break; + + case BoundaryTreatment::CentralSlope: + slope = (centerState - leftState) / (_LGLnodes[_nNodes - 1] - _LGLnodes[_nNodes - 2]); + break; + + case BoundaryTreatment::Constant: + return centerState; + } + } + else // Inner subcells + slope = _slope_limiter->call((rightState - centerState) / (_LGLnodes[subcellIdx + 1] - _LGLnodes[subcellIdx]), (centerState - leftState) / (_LGLnodes[subcellIdx] - _LGLnodes[subcellIdx - 1])); + + if (rightInterface == 1) + return centerState + slope * (_subcellGrid[subcellIdx + 1] - _LGLnodes[subcellIdx]); // Return state at right interface, i.e. forward flow upwind state + else + return centerState + slope * (_subcellGrid[subcellIdx] - _LGLnodes[subcellIdx]); // Return state at left interface, i.e. backward flow upwind state + } + } + + inline double LGLweights(int idx) { return _LGLweights[idx]; } + + private: + + std::unique_ptr _slope_limiter; //!< Slope limiter for second order FV reconstruction + unsigned int _FVorder; //!< FV subcell method order + BoundaryTreatment _FVboundaryTreatment; //!< boundary treatment of FV subcell method + double* _LGLweights; //!< DG method LGL weights, i.e. subcell sizes + double* _LGLnodes; //!< DG method LGL weights, i.e. subcell sizes + int _nComp; //!< Number of liquid phase components + int _nNodes; //!< Number of DG nodes per element, i.e. FV subcells + int _polyDeg; //!< Polynomial degree of DG Ansatz + double* _subcellGrid; //!< FV subcell grid, i.e. subcell borders + }; + +} // namespace cadet + +#endif // LIBCADET_DGSubcellLimiterFV_HPP_ diff --git a/src/libcadet/SmoothnessIndicator.hpp b/src/libcadet/SmoothnessIndicator.hpp new file mode 100644 index 000000000..6ff1edf29 --- /dev/null +++ b/src/libcadet/SmoothnessIndicator.hpp @@ -0,0 +1,366 @@ +// ============================================================================= +// CADET +// +// Copyright � 2008-2022: The CADET Authors +// Please see the AUTHORS and CONTRIBUTORS file. +// +// All rights reserved. This program and the accompanying materials +// are made available under the terms of the GNU Public License v3.0 (or, at +// your option, any later version) which accompanies this distribution, and +// is available at http://www.gnu.org/licenses/gpl.html +// ============================================================================= + +/** + * @file + * Implements smooothness indicators used in oscillation suppression for DG + */ + + +#ifndef LIBCADET_SMOOTHNESSINDICATOR_HPP_ +#define LIBCADET_SMOOTHNESSINDICATOR_HPP_ + +#include "AutoDiff.hpp" +#include +#include "ParamReaderHelper.hpp" + +namespace cadet +{ + /** + * @brief indicator to detect and assess DG elements troubled by oscillations and strong gradients + */ + class SmoothnessIndicator { + public: + virtual ~SmoothnessIndicator() {} + // todo: no templates for virtual functions. Is there another way to overload? + + /** + * @brief smoothness detector for DG element + * @param [in] localC pointer to first node of current DG element + * @param [in] strideNode node stride in state vector + * @param [in] strideLiquid liquid phase stride in state vector + * @param [in] cellIdx index of DG element + */ + virtual double calcSmoothness(const double* const localC, const int strideNode, const int strideLiquid, const int cellIdx) = 0; + /** + * @brief smoothness detector for DG element + * @param [in] localC pointer to first node of current DG element + * @param [in] strideNode node stride in state vector + * @param [in] strideLiquid liquid phase stride in state vector + * @param [in] cellIdx index of DG element + */ + virtual double calcSmoothness(const active* const localC, const int strideNode, const int strideLiquid, const int cellIdx) = 0; + /** + * @brief tells if and how much time relaxation is applied, i.e. change in smoothness indication is limited over time + */ + inline virtual double timeRelaxation() = 0; + + virtual int configure(IParameterProvider& parameterProvider) = 0; + }; + + /** + * @brief dummy indicator that identifies every element as troubled with max weight 1.0 + */ + class AllElementsIndicator : public SmoothnessIndicator { + + public: + + AllElementsIndicator() {} + + double calcSmoothness(const double* const localC, const int strideNode, const int strideLiquid, const int cellIdx) override { return 1.0; } + double calcSmoothness(const active* const localC, const int strideNode, const int strideLiquid, const int cellIdx) override { return 1.0; } + inline double timeRelaxation() override { return false; } + int configure(IParameterProvider& parameterProvider) override { return 1; }; + + }; + + /** + * @brief indicates troubled elements based on the relative modal energy of the two highes modes + */ + class ModalEnergyIndicator : public SmoothnessIndicator { + + public: + + ModalEnergyIndicator() {} + ModalEnergyIndicator(int polyDeg, const Eigen::MatrixXd& modalVanInv, int numModes, double nodalCoefThreshold, double nodalCoefShift, bool timeRelax) + : _polyDeg(polyDeg), _numModes(numModes), _modalVanInv(modalVanInv), _nodalCoefThreshold(nodalCoefThreshold), _nodalCoefShift(nodalCoefShift), _timeRelaxation(timeRelax) + { + _nNodes = _polyDeg + 1; + _modalCoeff.resize(_nNodes); + _modalCoeff.setZero(); + } + + inline double timeRelaxation() override { return _timeRelaxation; } + int configure(IParameterProvider& parameterProvider) override { return 1; }; + + double calcSmoothness(const double* const localC, const int strideNode, const int strideLiquid, const int cellIdx) override + { + // TODO: default constants, store them. + const double _s = 9.2102; + + Eigen::Map> _C(localC, _nNodes, Eigen::InnerStride(strideNode)); + + if (_C.cwiseAbs().maxCoeff() < _nodalCoefThreshold) + return 0.0; + else + _modalCoeff = _modalVanInv * (_C + _nodalCoefShift * Eigen::VectorXd::Ones(_nNodes)); + + double modalEnergySquareSum = _modalCoeff.cwiseAbs2().sum(); + double modeEnergy = (modalEnergySquareSum == 0.0) ? 0.0 : std::pow(_modalCoeff[_polyDeg], 2.0) / modalEnergySquareSum; + + if (_numModes == 2) { + modalEnergySquareSum = _modalCoeff.segment(0, _polyDeg).cwiseAbs2().sum(); + modeEnergy = std::max(modeEnergy, (modalEnergySquareSum == 0.0) ? 0.0 : std::pow(_modalCoeff[_polyDeg - 1], 2.0) / modalEnergySquareSum); + } + + const double _T = 0.5 * std::pow(10.0, -1.8 * std::pow(static_cast(_nNodes), 0.25)); // todo store this constant + return 1.0 / (1.0 + exp(-_s / _T * (modeEnergy - _T))); + } + double calcSmoothness(const active* const localC, const int strideNode, const int strideLiquid, const int cellIdx) override + { + // TODO: default constants, store them. + const double _s = 9.2102; + + Eigen::Map, 0, Eigen::InnerStride> _C(localC, _nNodes, Eigen::InnerStride(strideNode)); + + if (_C.template cast().cwiseAbs().maxCoeff() < _nodalCoefThreshold) + return 0.0; + else + _modalCoeff = _modalVanInv * (_C.template cast() + _nodalCoefShift * Eigen::VectorXd::Ones(_nNodes)); + + double modalEnergySquareSum = _modalCoeff.cwiseAbs2().sum(); + double modeEnergy = (modalEnergySquareSum == 0.0) ? 0.0 : std::pow(_modalCoeff[_polyDeg], 2.0) / modalEnergySquareSum; + + if (_numModes == 2) { + modalEnergySquareSum = _modalCoeff.segment(0, _polyDeg).cwiseAbs2().sum(); + modeEnergy = std::max(modeEnergy, (modalEnergySquareSum == 0.0) ? 0.0 : std::pow(_modalCoeff[_polyDeg - 1], 2.0) / modalEnergySquareSum); + } + + const double _T = 0.5 * std::pow(10.0, -1.8 * std::pow(static_cast(_nNodes), 0.25)); // 10-1.8(N + 1)0.2; // todo store this constant + return 1.0 / (1.0 + exp(-_s / _T * (modeEnergy - _T))); + } + + private: + + int _nNodes; //!< Number of polynomial interpolation nodes of DG ansatz + int _polyDeg; //!< Polynomial degree of DG ansatz + + double _nodalCoefThreshold; //!< Threshold for nodal polynomial coefficients + double _nodalCoefShift; //!< Shift for nodal polynomial coefficients (To get c(z) >= eps > 0 at constant 0 concentration values) + Eigen::VectorXd _modalCoeff; //!< Modal polynomial coefficient memory buffer + Eigen::MatrixXd _modalVanInv; //!< Inverse Vandermonde matrix of modal polynomial basis + int _numModes; //!< Number of modes to be considered in modal energy comparison + double _timeRelaxation; //!< Time relaxation of smoothness indicator, i.e. limit sudden changes by _timeRelaxation * old_value + }; + + /** + * @brief indicates troubled elements based on the relative modal energy of the two highes modes but gives a (binary) hard switch instead of a weight + */ + class ModalEnergyIndicatorHardSwitch : public SmoothnessIndicator { + + public: + + ModalEnergyIndicatorHardSwitch(int polyDeg, const Eigen::MatrixXd& modalVanInv, int numModes, double nodalCoefThreshold, double nodalCoefShift, bool timeRelax, double FVswitchAt) + : _FVswitchAt(FVswitchAt) + { + _ModalEnergyIndicator = std::make_unique(polyDeg, modalVanInv, numModes, nodalCoefThreshold, nodalCoefShift, timeRelax); + } + + inline double timeRelaxation() { return _ModalEnergyIndicator->timeRelaxation(); } + int configure(IParameterProvider& parameterProvider) override { return 1; }; + + double calcSmoothness(const double* const localC, const int strideNode, const int strideLiquid, const int cellIdx) override + { + return _ModalEnergyIndicator->calcSmoothness(localC, strideNode, strideLiquid, cellIdx) >= _FVswitchAt ? 1.0 : 0.0; + } + double calcSmoothness(const active* const localC, const int strideNode, const int strideLiquid, const int cellIdx) override + { + return _ModalEnergyIndicator->calcSmoothness(localC, strideNode, strideLiquid, cellIdx) >= _FVswitchAt ? 1.0 : 0.0; + } + + private: + + std::unique_ptr _ModalEnergyIndicator; + double _FVswitchAt; + }; + + /** + * @brief indicates troubled elements using slope limiters based on element-averaged and interface slopes + */ + class SlopeIndicator : public SmoothnessIndicator { + + /** + * @brief slope limiters for three (consecutive) slopes + */ + class IndicatorLimiter { + public: + virtual ~IndicatorLimiter() {} + virtual double call(double u0, double u1, double u2) = 0; + virtual active call(active u0, active u1, active u2) = 0; + }; + + /** + * @brief MINMOD slope limiter + */ + class MinmodIndicator : public IndicatorLimiter { + + public: + MinmodIndicator() {} + double call(double u0, double u1, double u2) override { + bool sgn = std::signbit(u0); + if (sgn == std::signbit(u1) && sgn == std::signbit(u2)) + return std::copysign(std::min(std::min(std::abs(u0), std::abs(u1)), std::abs(u1)), u0); + else + return 0.0; + } + active call(active u0, active u1, active u2) override { + bool sgn = std::signbit(static_cast(u0)); + if (sgn == std::signbit(static_cast(u1)) && sgn == std::signbit(static_cast(u2))) + { + if (abs(static_cast(u0)) < abs(static_cast(u1))) + { + if (abs(static_cast(u0)) < abs(static_cast(u2))) + return (sgn) ? -u0 : u0; + else + return (sgn) ? -u2 : u2; + } + else if (abs(static_cast(u1)) < abs(static_cast(u2))) + return (sgn) ? -u1 : u1; + else + return (sgn) ? -u2 : u2; + } + else + return active(0.0); + } + }; + + class TVBMinmodIndicator : public IndicatorLimiter { + + public: + TVBMinmodIndicator() {} + double call(double u0, double u1, double u2) override { + if (std::abs(u0) <= M * h * h) + return u0; + else + return minmod.call(u0, u1, u2); + } + active call(active u0, active u1, active u2) override { + if (std::abs(static_cast(u0)) <= M * h * h) + return u0; + else + return minmod.call(u0, u1, u2); + } + private: + MinmodIndicator minmod; + active h; + active M; + }; + + public: + + SlopeIndicator(int polyDeg, int nCells, double deltaZ, Eigen::VectorXd LGLweights, Eigen::MatrixXd polyDerM, std::string limiter, bool timeRelax) + : _polyDeg(polyDeg), _nCells(nCells), _deltaZ(deltaZ), _LGLweights(LGLweights), _polyDerM(polyDerM), _timeRelaxation(timeRelax) + { + _nNodes = _polyDeg + 1; + } + + inline double timeRelaxation() override { return _timeRelaxation; } // todo time relaxation required in this indicator which is either 1 or 0? + int configure(IParameterProvider& paramProvider) override + { + _deltaZ = paramProvider.getDouble("COL_LENGTH") / _nCells; + return 1; + }; + + double calcSmoothness(const double* const localC, const int strideNode, const int strideLiquid, const int cellIdx) override + { + Eigen::Map, 0, Eigen::InnerStride> _C(localC, _nNodes, Eigen::InnerStride(strideNode)); + + if (cellIdx > 0 && cellIdx < _nCells - 1) // todo boundary treatment + { + // todo store mass matrix and LGL weights additionally to respective inverse + // todo overwrite values instead of recalculation // todo deltaZ ParamType + _pAvg0 = 1.0 / static_cast(_deltaZ) * (_LGLweights.array() * _C.segment((cellIdx - 1) * _nNodes, _nNodes).array()).sum(); + _pAvg1 = 1.0 / static_cast(_deltaZ) * (_LGLweights.array() * _C.segment(cellIdx * _nNodes, _nNodes).array()).sum(); + _pAvg2 = 1.0 / static_cast(_deltaZ) * (_LGLweights.array() * _C.segment((cellIdx + 1) * _nNodes, _nNodes).array()).sum(); + + _uTilde = _C[cellIdx * _nNodes + _nNodes - 1] - _pAvg1; // average minus inner interface value on right face + _u2Tilde = _pAvg1 - _C[cellIdx * _nNodes]; // average minus inner interface value on left face + + double trigger1 = indicatorLimiter->call(_uTilde, _pAvg2 - _pAvg1, _pAvg1 - _pAvg0); + double trigger2 = indicatorLimiter->call(_u2Tilde, _pAvg2 - _pAvg1, _pAvg1 - _pAvg0); + + double M2 = (_polyDerM * _polyDerM * _C.segment(cellIdx * _nNodes, _nNodes)).maxCoeff(); + //double u_x = (_polyDerM * _C.segment(cell * _nNodes, _nNodes)).maxCoeff(); + double M_ = 2.0 / 3.0 * M2; + //double hmpf = 2.0 / 9.0 * (3.0 + 10.0 * M2) * M2 * _deltaZ / (_deltaZ + 2.0 * u_x * _deltaZ); + + // reconstruct if cell is troubled, i.e. potential oscillations + if (abs(trigger1 - _uTilde) > 1e-8 || abs(trigger2 - _u2Tilde) > 1e-8) + { + if (abs(_uTilde) > M_ * _deltaZ * _deltaZ || abs(_u2Tilde) > M_ * _deltaZ * _deltaZ) + return 1.0; + } + } + return 0.0; + } + double calcSmoothness(const active* const localC, const int strideNode, const int strideLiquid, const int cellIdx) override + { + Eigen::Map, 0, Eigen::InnerStride> _C(localC, _nNodes, Eigen::InnerStride(strideNode)); + + if (cellIdx > 0 && cellIdx < _nCells - 1) // todo boundary treatment + { + // todo store mass matrix and LGL weights additionally to respective inverse + // todo overwrite values instead of recalculation // todo deltaZ ParamType + _pAvg0 = 1.0 / static_cast(_deltaZ) * (_LGLweights.array() * _C.segment((cellIdx - 1) * _nNodes, _nNodes).template cast().array()).sum(); + _pAvg1 = 1.0 / static_cast(_deltaZ) * (_LGLweights.array() * _C.segment(cellIdx * _nNodes, _nNodes).template cast().array()).sum(); + _pAvg2 = 1.0 / static_cast(_deltaZ) * (_LGLweights.array() * _C.segment((cellIdx + 1) * _nNodes, _nNodes).template cast().array()).sum(); + + _uTilde = static_cast(_C[cellIdx * _nNodes + _nNodes - 1]) - _pAvg1; // average minus inner interface value on right face + _u2Tilde = _pAvg1 - static_cast(_C[cellIdx * _nNodes]); // average minus inner interface value on left face + + double trigger1 = indicatorLimiter->call(_uTilde, _pAvg2 - _pAvg1, _pAvg1 - _pAvg0); + double trigger2 = indicatorLimiter->call(_u2Tilde, _pAvg2 - _pAvg1, _pAvg1 - _pAvg0); + + double M2 = (_polyDerM * _polyDerM * _C.segment(cellIdx * _nNodes, _nNodes).template cast()).maxCoeff(); + //double u_x = (_polyDerM * _C.segment(cell * _nNodes, _nNodes)).maxCoeff(); + double M_ = 2.0 / 3.0 * M2; + //double hmpf = 2.0 / 9.0 * (3.0 + 10.0 * M2) * M2 * _deltaZ / (_deltaZ + 2.0 * u_x * _deltaZ); + + // reconstruct if cell is troubled, i.e. potential oscillations + if (abs(trigger1 - _uTilde) > 1e-8 || abs(trigger2 - _u2Tilde) > 1e-8) + { + if (abs(_uTilde) > M_ * _deltaZ * _deltaZ || abs(_u2Tilde) > M_ * _deltaZ * _deltaZ) + return 1.0; + } + } + return 0.0; + } + + double getCellAverage(int cellIdx) { + switch (cellIdx) + { + case 0: return _pAvg0; case 1: return _pAvg1; case 2: return _pAvg0; default: return 0.0; + } + } + + private: + + int _nNodes; //!< Number of polynomial interpolation nodes of DG ansatz + int _polyDeg; //!< Polynomial degree of DG ansatz + int _nCells; //!< Number of DG cells/elements + double _deltaZ; + Eigen::VectorXd _LGLweights; + Eigen::MatrixXd _polyDerM; + + std::unique_ptr indicatorLimiter; + + double _pAvg0; + double _pAvg1; + double _pAvg2; + double _uTilde; + double _u2Tilde; + + double _timeRelaxation; + }; +} + +#endif // LIBCADET_SMOOTHNESSINDICATOR_HPP_ diff --git a/src/libcadet/Weno_DG.hpp b/src/libcadet/WenoDG.hpp similarity index 54% rename from src/libcadet/Weno_DG.hpp rename to src/libcadet/WenoDG.hpp index ec473922c..de44e3091 100644 --- a/src/libcadet/Weno_DG.hpp +++ b/src/libcadet/WenoDG.hpp @@ -29,7 +29,6 @@ namespace cadet { - //template // todo class WENOLimiter { public: virtual ~WENOLimiter() {} @@ -60,12 +59,8 @@ namespace cadet return 0.0; } double call(active u0, active u1, active u2) override { - if (std::signbit(u0.getValue()) == std::signbit(u1.getValue()) && std::signbit(u0.getValue()) == std::signbit(u2.getValue())) { - if (std::signbit(u0.getValue())) // i.e. negative sign - return -static_cast(u0); - else - return static_cast(u0); - } + if (std::signbit(u0.getValue()) == std::signbit(u1.getValue()) && std::signbit(u0.getValue()) == std::signbit(u2.getValue())) + return std::copysign(std::min(std::abs(u0.getValue()), std::min(std::abs(u1.getValue()), std::abs(u2.getValue()))), u0.getValue()); else return 0.0; } @@ -107,21 +102,10 @@ namespace cadet { public: - /** - * @brief Boundary treatment method determines how the reconstruction handles insufficient available elements (i.e., less elements available than stencil size) - */ - enum class BoundaryTreatment : int - { - ReduceOrder = 0, //!< Reduce the order of the WENO method such that the stencil is small enough - ZeroWeights = 1, //!< Set weights of WENO method to 0 for unavailable elements - ZeroWeightsForPnotZero = 2, //!< Set weights of WENO method to 0 for unavailable elements, except for the first cell (order is reduced to 1) - LargeGhostNodes = 3, - }; - /** * @brief Creates the WENO scheme */ - //WenoDG() : _order(maxOrder()), _boundaryTreatment(BoundaryTreatment::ReduceOrder), _intermediateValues(3 * maxOrder() * sizeof(active)) { } + WenoDG() : _weno_gamma(nullptr), _w(nullptr), _wHat(nullptr), _beta(nullptr), weno_limiter(nullptr) { } ~WenoDG() CADET_NOEXCEPT { @@ -129,32 +113,85 @@ namespace cadet delete[] _w; delete[] _wHat; delete[] _beta; - delete[] _troubled_cells; } - void init(const double eps, const double r, const double gamma, const unsigned int nCells, const unsigned int nComp, const bool trackTroubledCells = false) { + void init(std::string limiter, const double eps, const double r, const double gamma, double* invWeights, Eigen::MatrixXd polyDerM, const unsigned int nCells, const unsigned int nNodes, const unsigned int nComp) { + + _nComp = nComp; + _nCells = nCells; + _nNodes = nNodes; - _weno_gamma = new double[3]{ (1 - gamma) / 2.0, gamma, (1 - gamma) / 2.0 }; + _polyDerM = polyDerM; + _LGLweights.resize(nNodes); + for (int node = 0; node < nNodes; node++) + _LGLweights[node] = 1.0 / invWeights[node]; + + _weno_gamma = new double[3]{ (1.0 - gamma) / 2.0, gamma, (1.0 - gamma) / 2.0 }; _weno_r = r; _weno_eps = eps; - if (trackTroubledCells) - _troubled_cells = new double[nCells * nComp]; - else - _troubled_cells = nullptr; + + if (limiter == "MINMOD") + weno_limiter = std::make_unique(); + else if (limiter == "TVB_MINMOD") + weno_limiter = std::make_unique(); } - inline double* troubledCells() const CADET_NOEXCEPT { return &_troubled_cells[0]; } + void setDeltaZ(active deltaZ) { _deltaZ = deltaZ; } + + template + void calcSmoothness(const StateType* const localC, const int strideNode, const int comp, const int cell, double* troubled_cells) + { + troubled_cells[comp + cell * _nComp] = 2.0; + + Eigen::Map, 0, InnerStride> _C(localC, _nNodes, InnerStride(strideNode)); + + if (cell > 0 && cell < _nCells - 1) // todo boundary treatment + { + // todo store mass matrix and LGL weights additionally to respective inverse + // todo overwrite values instead of recalculation // todo deltaZ ParamType + _pAvg0 = 1.0 / static_cast(_deltaZ) * (_LGLweights.array() * _C.segment((cell - 1) * _nNodes, _nNodes).template cast().array()).sum(); + _pAvg1 = 1.0 / static_cast(_deltaZ) * (_LGLweights.array() * _C.segment(cell * _nNodes, _nNodes).template cast().array()).sum(); + _pAvg2 = 1.0 / static_cast(_deltaZ) * (_LGLweights.array() * _C.segment((cell + 1) * _nNodes, _nNodes).template cast().array()).sum(); + + _uTilde = _C[cell * _nNodes + _nNodes - 1] - _pAvg1; // average minus inner interface value on right face + _u2Tilde = _pAvg1 - _C[cell * _nNodes]; // average minus inner interface value on left face + + double trigger1 = weno_limiter->call(_uTilde, _pAvg2 - _pAvg1, _pAvg1 - _pAvg0); + double trigger2 = weno_limiter->call(_u2Tilde, _pAvg2 - _pAvg1, _pAvg1 - _pAvg0); + + double M2 = (_polyDerM * _polyDerM * _C.segment(cell * _nNodes, _nNodes).template cast()).maxCoeff(); + //double u_x = (_polyDerM * _C.segment(cell * _nNodes, _nNodes)).maxCoeff(); + double M_ = 2.0 / 3.0 * M2; + //double hmpf = 2.0 / 9.0 * (3.0 + 10.0 * M2) * M2 * _deltaZ / (_deltaZ + 2.0 * u_x * _deltaZ); + + // reconstruct if cell is troubled, i.e. potential oscillations + if (abs(trigger1 - _uTilde) > 1e-8 || abs(trigger2 - _u2Tilde) > 1e-8) + { + if (abs(_uTilde) > M_ * _deltaZ * _deltaZ || abs(_u2Tilde) > M_ * _deltaZ * _deltaZ) + // mark troubled cell + troubled_cells[comp + cell * _nComp] = 1.0; + } + } + } //private: // indicator for troubled cells std::unique_ptr weno_limiter; + int _nNodes; + int _nComp; + int _nCells; + // WENO parameters double* _weno_gamma; double _weno_r = 2.0; double _weno_eps = 1e-6; + // weights, smoothness indicator, integral values. + Eigen::VectorXd _LGLweights; //!< DG method LGL weights, i.e. subcell sizes + Eigen::MatrixXd _polyDerM; //!< DG method LGL weights, i.e. subcell sizes + active _deltaZ; active* _wHat = new active[3]{ 0.0, 0.0, 0.0 }; active* _w = new active[3]{ 0.0, 0.0, 0.0 }; active _wSum = 0.0; @@ -167,8 +204,6 @@ namespace cadet active _avgUdeltaMinus = 0.0; active _uTilde = 0.0; active _u2Tilde = 0.0; - // mark troubled cells (smoothness indicator) - double* _troubled_cells; }; diff --git a/src/libcadet/model/LumpedRateModelWithPoresDG-InitialConditions.cpp b/src/libcadet/model/LumpedRateModelWithPoresDG-InitialConditions.cpp index e3941ee4e..05de64559 100644 --- a/src/libcadet/model/LumpedRateModelWithPoresDG-InitialConditions.cpp +++ b/src/libcadet/model/LumpedRateModelWithPoresDG-InitialConditions.cpp @@ -343,7 +343,7 @@ namespace cadet LinearBufferAllocator tlmAlloc = threadLocalMem.get(); // Reuse memory of band matrix for dense matrix - linalg::DenseMatrixView fullJacobianMatrix(_globalJacDisc.valuePtr() + _globalJacDisc.outerIndexPtr()[idxr.offsetCp(ParticleTypeIndex{ type }) - idxr.offsetC() + pblk], nullptr, mask.len, mask.len); + linalg::DenseMatrixView fullJacobianMatrix(_globalJacDisc.valuePtr() + _globalJacDisc.outerIndexPtr()[idxr.offsetCp(ParticleTypeIndex{ type }) - idxr.offsetC() + pblk * idxr.strideParBlock(pblk)], nullptr, mask.len, mask.len); // z coordinate (column length normed to 1) of current node - needed in externally dependent adsorption kinetic const double z = _convDispOp.relativeCoordinate(pblk); @@ -362,7 +362,7 @@ namespace cadet double* const fullX = static_cast(fullXBuffer); BufferedArray jacobianMemBuffer = tlmAlloc.array(probSize * probSize); - linalg::DenseMatrixView jacobianMatrix(static_cast(jacobianMemBuffer), _globalJacDisc.outerIndexPtr() + pblk * probSize, probSize, probSize); + linalg::DenseMatrixView jacobianMatrix(static_cast(jacobianMemBuffer), /*_globalJacDisc.outerIndexPtr() + pblk * probSize*/nullptr, probSize, probSize); BufferedArray conservedQuantsBuffer = tlmAlloc.array(numActiveComp); double* const conservedQuants = static_cast(conservedQuantsBuffer); diff --git a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp index 87a438aa4..d73532afe 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp @@ -17,8 +17,6 @@ #include "ParamReaderHelper.hpp" #include "AdUtils.hpp" #include "SimulationTypes.hpp" -//#include "model/parts/AxialConvectionDispersionKernelDG.hpp" // todo radial flow DG and outsource residual implementation to kernel -//#include "model/parts/RadialConvectionDispersionKernelDG.hpp" #include "model/ParameterDependence.hpp" #include "SensParamUtil.hpp" #include "ConfigurationHelper.hpp" @@ -44,8 +42,7 @@ namespace parts * @brief Creates an AxialConvectionDispersionOperatorBaseDG */ AxialConvectionDispersionOperatorBaseDG::AxialConvectionDispersionOperatorBaseDG() : - /*_stencilMemory(sizeof(active) * Weno::maxStencilSize()), _wenoDerivatives(new double[Weno::maxStencilSize()]), _weno(),*/ // todo weno - _dispersionDep(nullptr) + _dispersionDep(nullptr), _DGjacAxDispBlocks(nullptr) { } @@ -105,7 +102,7 @@ bool AxialConvectionDispersionOperatorBaseDG::configureModelDiscretization(IPara invMMatrix(); derivativeMatrix(); - if(polynomial_integration_mode == 2) // use Gauss quadrature for exact integration + if (polynomial_integration_mode == 2) // use Gauss quadrature for exact integration _invMM = gaussQuadratureMMatrix(_nodes, _nNodes).inverse(); if (paramProvider.exists("COL_DISPERSION_DEP")) @@ -122,12 +119,68 @@ bool AxialConvectionDispersionOperatorBaseDG::configureModelDiscretization(IPara paramProvider.pushScope("discretization"); _OSmode = 0; + _maxBlending = 1.0; + _blendingThreshold = 0.0; + _troubledCells = nullptr; + // Read oscillation suppression settings and apply them if (paramProvider.exists("oscillation_suppression")) { paramProvider.pushScope("oscillation_suppression"); - std::string osMethod = paramProvider.getString("METHOD"); + /* Configure smoothness/troubled cell indicator, i.e. shock detection */ + + const bool write_smoothness_indicator = paramProvider.exists("RETURN_SMOOTHNESS_INDICATOR") ? static_cast(paramProvider.getInt("RETURN_SMOOTHNESS_INDICATOR")) : false; + if (write_smoothness_indicator) + { + _troubledCells = new double[nCells * nComp]; + std::fill(_troubledCells, _troubledCells + nCells * nComp, 0.0); + } + + std::string smoothness_indicator = paramProvider.exists("SMOOTHNESS_INDICATOR") ? paramProvider.getString("SMOOTHNESS_INDICATOR") : "MODAL_ENERGY_LEGENDRE"; + + if (smoothness_indicator == "MODAL_ENERGY_LEGENDRE" || smoothness_indicator == "MODAL_ENERGY_LEGENDRE_HARD_SWITCH") + { + _blendingThreshold = paramProvider.exists("BLENDING_THRESHOLD") ? paramProvider.getDouble("BLENDING_THRESHOLD") : 0.0; // 0.01; // todo default value + if (_blendingThreshold < 0.0 || _blendingThreshold > 1.0) + throw InvalidParameterException("Blending threshold must be >= 0.0 and <= 1.0, but was specified as " + std::to_string(_blendingThreshold)); + + double nodalCoefThreshold = paramProvider.exists("NODAL_COEFFICIENT_THRESHOLD") ? paramProvider.getDouble("NODAL_COEFFICIENT_THRESHOLD") : 0.0; // 1e-12; // todo default value + double nodalCoefShift = paramProvider.exists("NODAL_COEFFICIENT_SHIFT") ? paramProvider.getDouble("NODAL_COEFFICIENT_SHIFT") : 0.0; // 0.1; // todo default value + int numModes = paramProvider.exists("NUM_MODES") ? paramProvider.getInt("NUM_MODES") : 2; // todo default value + if (numModes != 1 && numModes != 2) + throw InvalidParameterException("Number of modes to compare modal energy must be 1 or 2, but was specified as " + std::to_string(numModes)); + + double timeRelaxation = paramProvider.exists("TIME_RELAXATION") ? paramProvider.getDouble("TIME_RELAXATION") : 0.0; // todo default value (0.7 was chosen in https://doi.org/10.1016/j.jcp.2021.110580) + + if (smoothness_indicator == "MODAL_ENERGY_LEGENDRE") + _smoothnessIndicator = std::make_unique(_polyDeg, getVandermonde_LEGENDRE().inverse(), numModes, nodalCoefThreshold, nodalCoefShift, timeRelaxation); + else if (smoothness_indicator == "MODAL_ENERGY_LEGENDRE_HARD_SWITCH") + { + if (paramProvider.exists("FV_SWITCH_AT")) + { + double FVswitchAt = paramProvider.getDouble("FV_SWITCH_AT"); + _smoothnessIndicator = std::make_unique(_polyDeg, getVandermonde_LEGENDRE().inverse(), numModes, nodalCoefThreshold, nodalCoefShift, timeRelaxation, FVswitchAt); + } + else + throw InvalidParameterException("Field FV_SWITCH_AT must be specifiec for MODAL_ENERGY_LEGENDRE_HARD_SWITCH smoothness indicator"); + } + } + else if (smoothness_indicator == "SLOPE_LIMITER") + { + const std::string limiter = paramProvider.exists("INDICATOR_LIMITER") ? paramProvider.getString("INDICATOR_LIMITER") : "MINMOD"; // todo choose default + // todo time relaxation not required for this indicator (either 1 or zero) ? + //double timeRelaxation = paramProvider.exists("TIME_RELAXATION") ? paramProvider.getDouble("TIME_RELAXATION") : 0.0; // todo default value (0.7 was chosen in https://doi.org/10.1016/j.jcp.2021.110580) + double timeRelaxation = 0.0; + // note that deltaZ is set in configure() + _smoothnessIndicator = std::make_unique(_polyDeg, _nCells, 0.0, _invWeights.cwiseInverse(), _polyDerM, limiter, timeRelaxation); + } + else if (smoothness_indicator == "ALL_ELEMENTS") + _smoothnessIndicator = std::make_unique(); + else + throw InvalidParameterException("Unknown smoothness indicator " + smoothness_indicator); + + std::string osMethod = paramProvider.getString("METHOD"); if (osMethod == "WENO") { @@ -137,25 +190,49 @@ bool AxialConvectionDispersionOperatorBaseDG::configureModelDiscretization(IPara double eps = paramProvider.exists("WENO_EPS") ? paramProvider.getDouble("WENO_EPS") : 1e-6; double r = paramProvider.exists("WENO_R") ? paramProvider.getDouble("WENO_R") : 2.0; double gamma = paramProvider.exists("WENO_GAMMA") ? paramProvider.getDouble("WENO_GAMMA") : 0.998; - bool write_smoothness_indicator = paramProvider.exists("RETURN_SMOOTHNESS_INDICATOR") ? static_cast(paramProvider.getInt("RETURN_SMOOTHNESS_INDICATOR")) : false; - _weno.init(eps, r, gamma, _nCells, _nComp, write_smoothness_indicator); + std::string limiter = paramProvider.exists("WENO_LIMITER") ? paramProvider.getString("WENO_LIMITER") : "MINMOD"; + + _weno.init(limiter, eps, r, gamma, &_invWeights[0], _polyDerM, _nCells, _nNodes, _nComp); } - else if (osMethod == "SUBCELL_LIMITING") + else if (osMethod == "SUBCELL_FV_LIMITING" || osMethod == "SUBCELL_FV") { - LOG(Debug) << "Subcell limiting used as oscillation suppression mechanism."; + if (osMethod == "SUBCELL_FV_LIMITING") { + LOG(Debug) << "Subcell FV limiting used as oscillation suppression mechanism for advection."; + _OSmode = 2; + } + else if (osMethod == "SUBCELL_FV") { + LOG(Debug) << "Subcell FV for advection and dispersion used on DG grid."; + _OSmode = 4; + } + + const int _FVorder = paramProvider.exists("SUBCELL_FV_ORDER") ? paramProvider.getInt("SUBCELL_FV_ORDER") : 1; // todo choose default + if (_FVorder <= 0) + throw InvalidParameterException("Subcell FV order must be > 0, but was specified as " + std::to_string(_FVorder)); + + _maxBlending = paramProvider.exists("SUBCELL_FV_MAX_BLENDING") ? paramProvider.getDouble("SUBCELL_FV_MAX_BLENDING") : 1.0; // todo choose default + if(_maxBlending < 0.0 || _maxBlending > 1.0) + throw InvalidParameterException("Low order blending factor must be >= 0.0 and <= 1.0, but was specified as " + std::to_string(_maxBlending)); - _OSmode = 2; - int order = paramProvider.getInt("FV_ORDER"); - // @TODO Subcell limiting - throw InvalidParameterException("Oscillation suppression mechanism " + osMethod + " not implemented yet"); + const std::string limiter = paramProvider.exists("SUBCELL_FV_LIMITER") ? paramProvider.getString("SUBCELL_FV_LIMITER") : "MINMOD"; // todo choose default + const int _FVboundaryTreatment = paramProvider.exists("SUBCELL_FV_BOUNDARY_TREATMENT") ? paramProvider.getInt("SUBCELL_FV_BOUNDARY_TREATMENT") : 0; // todo choose default + + _subcellLimiter.init(limiter, _FVorder, _FVboundaryTreatment, _nNodes, &_nodes[0], &_invWeights[0], _nCells, _nComp); } else if (osMethod != "NONE") throw InvalidParameterException("Unknown oscillation suppression mechanism " + osMethod + " in oscillation_suppression_mode"); + _calc_smoothness_indicator = write_smoothness_indicator || (_OSmode != 0 && _OSmode != 4); + paramProvider.popScope(); } + else + { + _smoothnessIndicator = nullptr; + _calc_smoothness_indicator = false; + } + paramProvider.popScope(); return true; @@ -174,6 +251,10 @@ bool AxialConvectionDispersionOperatorBaseDG::configure(UnitOpIdx unitOpIdx, IPa // Read geometry parameters _colLength = paramProvider.getDouble("COL_LENGTH"); _deltaZ = _colLength / _nCells; + if (_OSmode == 1) + _weno.setDeltaZ(_deltaZ); + if (_OSmode) + _smoothnessIndicator->configure(paramProvider); /* compute dispersion jacobian blocks(without parameters except element spacing, i.e. static entries) */ // we only need unique dispersion blocks, which are given by cells 1, 2, nCol for inexact integration DG and by cells 1, 2, 3, nCol-1, nCol for eaxct integration DG @@ -319,6 +400,8 @@ bool AxialConvectionDispersionOperatorBaseDG::notifyDiscontinuousSectionTransiti _curVelocity *= -1.0; } + _subcellLimiter.notifyDiscontinuousSectionTransition(_dir); + // Remaining case: _velocity is empty and _crossSection <= 0.0 // _curVelocity is goverend by network flow rate provided in setFlowRates(). // Direction never changes (always forward, that is, _dir = 1)- @@ -419,41 +502,36 @@ int AxialConvectionDispersionOperatorBaseDG::residualImpl(const IModel& model, d { for (unsigned int comp = 0; comp < _nComp; comp++) { - // create Eigen objects + _boundary[0] = y[comp]; // copy inlet DOFs to ghost node + + // Create Eigen objects Eigen::Map, 0, InnerStride> _C(y + offsetC() + comp, _nPoints, InnerStride(_strideNode)); Eigen::Map, 0, InnerStride> _resC(res + offsetC() + comp, _nPoints, InnerStride(_strideNode)); Eigen::Map, 0, InnerStride<>> _h(reinterpret_cast(&_h[0]), _nPoints, InnerStride<>(1)); Eigen::Map, 0, InnerStride<>> _g(reinterpret_cast(&_g[0]), _nPoints, InnerStride<>(1)); - ////ResidualType _pAvg0 = reinterpret_cast(_weno._pAvg0); - ////ResidualType _pAvg1 = reinterpret_cast(_weno._pAvg1); - ////ResidualType _pAvg2 = reinterpret_cast(_weno._pAvg2); - - //// calculate smoothness indicator, if oscillation suppression is active - //if (_OSmode != 0) - //{ - // if (_OSmode == 1) // WENO - // { - // for (unsigned int cell = 0; cell < _nCells; cell++) - // { - // if (cell > 0 && cell < _nCells - 1) // todo boundary treatment - // { - // // todo store weights additionally to inverse - // // todo overwrite values instead of recalculation - // _weno._pAvg0 = 1.0 / static_cast(_deltaZ) * (_invWeights.cwiseInverse().cast().array() * _C.segment((cell - 1) * _nNodes, _nNodes).array()).sum(); - // _weno._pAvg1 = 1.0 / static_cast(_deltaZ) * (_invWeights.cwiseInverse().cast().array() * _C.segment(cell * _nNodes, _nNodes).array()).sum(); - // _weno._pAvg2 = 1.0 / static_cast(_deltaZ) * (_invWeights.cwiseInverse().cast().array() * _C.segment((cell + 1) * _nNodes, _nNodes).array()).sum(); - - - // // mark troubled cell - // if (true) - // { - // int hm = 0; - // } - // } - // } - // } - //} + // Calculate smoothness indicator if oscillation suppression is active + if (_calc_smoothness_indicator) + { + double oldIndicator = 0.0; + + for (unsigned int cell = 0; cell < _nCells; cell++) + { + oldIndicator = _troubledCells[comp + cell * _nComp]; + + _troubledCells[comp + cell * _nComp] = std::min(_maxBlending, _smoothnessIndicator->calcSmoothness(y + offsetC() + comp + cell * _strideCell, _strideNode, _nComp, cell)); + if (_troubledCells[comp + cell * _nComp] < _blendingThreshold) + _troubledCells[comp + cell * _nComp] = 0.0; + else if (1.0 - _troubledCells[comp + cell * _nComp] < _blendingThreshold) + _troubledCells[comp + cell * _nComp] = 1.0; + else if (_smoothnessIndicator->timeRelaxation()) // note: 0.0 <-> no time relaxation + _troubledCells[comp + cell * _nComp] = std::max(_smoothnessIndicator->timeRelaxation() * oldIndicator, _troubledCells[comp + cell * _nComp]); + } + //if (_OSmode == 1) // WENO + //{ + // // todo: calculate and store reconstructed polynomial in h, calculate g from h and then set h = 2.0 / deltaZ * (-u * h + d_ax * (-2.0 / deltaZ) * g); + //} + } // Add time derivative to bulk residual if (yDot) @@ -467,49 +545,73 @@ int AxialConvectionDispersionOperatorBaseDG::residualImpl(const IModel& model, d const ParamType u = static_cast(_curVelocity); const ParamType d_ax = static_cast(getSectionDependentSlice(_colDispersion, _nComp, secIdx)[comp]); + double DGblending = 1.0 - _maxBlending; + + if (_OSmode == 2 && _maxBlending > 0.0) // Element-wise lower order subcell FV blending for advection, treat dispersion with high order DG. + subcellFVconvectionIntegral(_troubledCells + comp, y + offsetC() + comp, res + offsetC() + comp); + // else if (_OSmode == 3) // todo ? subcell limiting with subelement-wise blending + + else if (_OSmode == 4) // Subcell FV for both, advection and dispersion and without blending, i.e. pure FV // todo enable blending? -> ensure mass conservation for dispersive flux across DG elements + { + subcellFVconvDispIntegral(y + offsetC() + comp, res + offsetC() + comp, d_ax); + continue; + } + // ===================================// - // reset cache // + // Reset cache // // ===================================// _h.setZero(); _g.setZero(); - _boundary[0] = y[comp]; // copy inlet DOFs to ghost node - // ======================================// - // solve auxiliary system g = d c / d x // - // ======================================// + // ========================================// + // Compute auxiliary system g = d c / d x // + // ========================================// - // DG volume integral in strong form + // DG volume integral in strong form volumeIntegral(_C, _g); - // calculate numerical flux values c* + // Calculate numerical flux values c* InterfaceFluxAuxiliary(y + offsetC() + comp, _strideNode, _strideCell); // DG surface integral in strong form - surfaceIntegral(y + offsetC() + comp, reinterpret_cast(&_g[0]), - _strideNode, _strideCell, 1u, _nNodes); + surfaceIntegral(y + offsetC() + comp, reinterpret_cast(&_g[0]), _strideNode, _strideCell, 1u, _nNodes); - // ======================================// - // solve main equation RHS d h / d x // - // ======================================// + // ========================================// + // Compute main equation RHS d h / d x // + // ========================================// - // calculate the substitute h(g(c), c) and apply inverse mapping jacobian (reference space) - _h = 2.0 / static_cast(_deltaZ) * (-u * _C + d_ax * (-2.0 / static_cast(_deltaZ)) * _g).template cast(); + // Calculate the substitute h(g(c), c) and apply inverse mapping jacobian (reference space) + if (_OSmode != 0) + { + if (DGblending < 1.0) + { + for (int cell = 0; cell < _nCells; cell++) + _h.segment(cell * _nNodes, _nNodes) = 2.0 / static_cast(_deltaZ) * (-u * _C.segment(cell * _nNodes, _nNodes) * (1.0 - _troubledCells[comp + cell * _nComp])).template cast(); + _h += (2.0 / static_cast(_deltaZ) * d_ax * (-2.0 / static_cast(_deltaZ)) * _g).template cast(); + } + else + _h = 2.0 / static_cast(_deltaZ) * (-u * _C * DGblending + d_ax * (-2.0 / static_cast(_deltaZ)) * _g).template cast(); + } + else + _h = 2.0 / static_cast(_deltaZ) * (-u * _C + d_ax * (-2.0 / static_cast(_deltaZ)) * _g).template cast(); // DG volume integral in strong form volumeIntegral(_h, _resC); - // update boundary values for auxiliary variable g (solid wall) + // Update boundary values for auxiliary variable g (solid wall) calcBoundaryValues(); - // calculate numerical flux values h* + // Calculate numerical flux values h* InterfaceFlux(y + offsetC() + comp, d_ax); + // By leaving the following out, the surface integral computes B (h^* - \hat{h}) with \hat{h} = h - \alpha u c, i.e. \hat{h} equals the substitute h minus the subcell (convection) part. + // This way we can leave out any interface flux computations at DG-element boundaries in the subcell FV scheme. + // _h += 2.0 / static_cast(_deltaZ) * (-u * _C * FVblending).template cast(); + // DG surface integral in strong form - surfaceIntegral(&_h[0], res + offsetC() + comp, - 1u, _nNodes, _strideNode, _strideCell); + surfaceIntegral(&_h[0], res + offsetC() + comp, 1u, _nNodes, _strideNode, _strideCell); } - return 0; } /** @@ -553,9 +655,9 @@ unsigned int AxialConvectionDispersionOperatorBaseDG::nConvDispEntries(bool pure void model::parts::AxialConvectionDispersionOperatorBaseDG::convDispJacPattern(std::vector& tripletList, const int bulkOffset) { if (_exactInt) - ConvDispModalPattern(tripletList, bulkOffset); + ConvDispExIntDGPattern(tripletList, bulkOffset); else - ConvDispNodalPattern(tripletList, bulkOffset); + ConvDispCollocationDGPattern(tripletList, bulkOffset); } /** * @brief Multiplies the time derivative Jacobian @f$ \frac{\partial F}{\partial \dot{y}}\left(t, y, \dot{y}\right) @f$ with a given vector diff --git a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp index 7c821ec3e..7bd845bb7 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp @@ -21,7 +21,9 @@ #include "ParamIdUtil.hpp" #include "AutoDiff.hpp" #include "Memory.hpp" -#include "Weno_DG.hpp" +#include "SmoothnessIndicator.hpp" +#include "WenoDG.hpp" +#include "DGSubcellLimiterFV.hpp" #include "SimulationTypes.hpp" #include #include "linalg/BandedEigenSparseRowIterator.hpp" @@ -124,13 +126,11 @@ namespace cadet inline unsigned int nNodes() const CADET_NOEXCEPT { return _nNodes; } inline unsigned int nPoints() const CADET_NOEXCEPT { return _nPoints; } inline bool exactInt() const CADET_NOEXCEPT { return _exactInt; } - inline bool hasSmoothnessIndicator() const CADET_NOEXCEPT { return static_cast(_OSmode); } // only zero if no oscillation suppression + inline bool hasSmoothnessIndicator() const CADET_NOEXCEPT { return _calc_smoothness_indicator; } inline double* smoothnessIndicator() const CADET_NOEXCEPT { - if (_OSmode == 1) - return _weno.troubledCells(); - //else if (_OSmode == 2) // todo subcell limiting - // return _subcell.troubledCells(); + if (hasSmoothnessIndicator()) + return _troubledCells; else return nullptr; } @@ -191,8 +191,13 @@ namespace cadet // non-linear oscillation prevention mechanism int _OSmode; //!< oscillation suppression mode; 0 : none, 1 : WENO, 2 : Subcell limiting + double _maxBlending; //!< maximal blending coefficient of oscillation suppression mechanism + double _blendingThreshold; //!< Threshold to clip-off blending coefficient (in both directions) WenoDG _weno; //!< WENO operator - //SucellLimiter _subcellLimiter; // todo + DGSubcellLimiterFV _subcellLimiter; //!< FV subcell limiting operator + std::unique_ptr _smoothnessIndicator; //!< smoothness/troubled-element indicator + double* _troubledCells; //!< troubled-element indicator values + bool _calc_smoothness_indicator; //!< Determines whether or not the smoothness indicator is evaluated // Simulation parameters active _colLength; //!< Column length \f$ L \f$ @@ -204,17 +209,12 @@ namespace cadet active _curVelocity; //!< Current interstitial velocity \f$ u \f$ in this time section int _dir; //!< Current flow direction in this time section - // needed? + // todo still needed? int _curSection; //!< current section index bool _newStaticJac; //!< determines wether static analytical jacobian needs to be computed (every section) - // todo weno - //ArrayPool _stencilMemory; //!< Provides memory for the stencil - //double* _wenoDerivatives; //!< Holds derivatives of the WENO scheme - //Weno _weno; //!< The WENO scheme implementation - //double _wenoEpsilon; //!< The @f$ \varepsilon @f$ of the WENO scheme (prevents division by zero) - bool _dispersionCompIndep; //!< Determines whether dispersion is component independent + IParameterParameterDependence* _dispersionDep; /* =================================================================================== @@ -249,6 +249,7 @@ namespace cadet /** * @brief computes the Legendre-Gauss-Lobatto nodes and (inverse) quadrature weights + * @detail inexact LGL-quadrature leads to a diagonal mass matrix (mass lumping), defined by the quadrature weights */ void lglNodesWeights() { // tolerance and max #iterations for Newton iteration @@ -463,7 +464,7 @@ namespace cadet massMatrix(i, j) = (aux.array() * LGweigths.array()).sum(); } } - + return massMatrix; } @@ -910,16 +911,165 @@ namespace cadet else { // collocated numerical integration -> diagonal mass matrix for (unsigned int Cell = 0; Cell < _nCells; Cell++) { // strong surface integral -> M^-1 B [state - state*] - stateDer[Cell * strideCell_stateDer] // first cell, node + stateDer[Cell * strideCell_stateDer] // first node -= static_cast(_invWeights[0] * (state[Cell * strideCell_state] - _surfaceFlux(Cell))); - stateDer[Cell * strideCell_stateDer + _polyDeg * strideNode_stateDer] // last cell, node + stateDer[Cell * strideCell_stateDer + _polyDeg * strideNode_stateDer] // last node += static_cast(_invWeights[_polyDeg] * (state[Cell * strideCell_state + _polyDeg * strideNode_state] - _surfaceFlux(Cell + 1))); } } } + /** + * @brief calculates the subcell FV fluxes for convection and dispersion + * @param [in] blending pointer to DG element-wise blending coeffcients + * @param [in] _C pointer to state vector at component of interest + * @param [in] stateDer pointer to residual vector at component of interest + * @param [in] d_ax axial dispersion coefficient of current component + */ + template + void subcellFVconvDispIntegral(const StateType* _C, ResidualType* stateDer, ParamType d_ax) + { + + double blending = 1.0; // no blending with DG possible for this implementation using FD approximation of dispersion + const StateType* localC = _C; + StateType upwindState = 0.0; + ResidualType* localStateDer = stateDer; + ResidualType flux = 0.0; + ResidualType derivative = 0.0; + + bool forward = (_dir == 1); + // Subcell inner faces + for (unsigned int cell = 0; cell < _nCells; cell++, localC += _strideNode, localStateDer += _strideNode) { + + for (int interface = 1; interface < _nNodes; interface++, localC += _strideNode, localStateDer += _strideNode) { + + if (cell == 0 && interface == 1) + upwindState = _subcellLimiter.reconstructedInterfaceValue(static_cast(_boundary[0]), localC[0], localC[_strideNode], interface - 1, forward); + else if (cell == _nCells - 1 && interface == _nNodes - 1) + upwindState = _subcellLimiter.reconstructedInterfaceValue(*(localC - _strideNode), localC[0], localC[0], interface - 1, forward); + else + upwindState = _subcellLimiter.reconstructedInterfaceValue(*(localC - _strideNode), localC[0], localC[_strideNode], interface - 1, forward); + + // Approximate derivative by central finite difference. Only central (2nd order) for N_d = 1, otherwise some weighted central. We have distance between the two points as h = 0.5 * summed cell sizes * reference element mapping + derivative = (localC[_strideNode] - localC[0]) / (0.5 * (_subcellLimiter.LGLweights(interface - 1) + _subcellLimiter.LGLweights(interface)) * 0.5 * static_cast(_deltaZ)); + + flux = blending * (2.0 / static_cast(_deltaZ)) * // blending coefficient and mapping + (static_cast(_curVelocity) * upwindState // upwind flux for convection + - d_ax * derivative); // central FD approximation of c_x + + // add flux divided by respective cell volume. Note that mapping was already applied + localStateDer[0] += flux * _invWeights[interface - 1]; + localStateDer[_strideNode] -= flux * _invWeights[interface]; + } + } + + // DG element faces + // Inlet boundary BC + flux = blending * (2.0 / static_cast(_deltaZ)) * static_cast(_curVelocity * _boundary[0]); + if (_dir == 1) + stateDer[0] -= flux * _invWeights[0]; + else + stateDer[(_nPoints - 1) * _strideNode] -= flux * _invWeights[_polyDeg]; + + // Inner DG element faces + localC = _C + _strideNode * (_nNodes - 1); + localStateDer = stateDer + _strideNode * (_nNodes - 1); + + for (unsigned int elemFace = 1; elemFace < _nCells; elemFace++, localC += _nNodes * _strideNode, localStateDer += _nNodes * _strideNode) { + + upwindState = _subcellLimiter.reconstructedInterfaceValue(*(localC - _strideNode), localC[0], localC[_strideNode], _polyDeg, forward); + + derivative = (localC[_strideNode] - localC[0]) / (0.5 * (_subcellLimiter.LGLweights(0) + _subcellLimiter.LGLweights(_polyDeg)) * 0.5 * static_cast(_deltaZ)); + + flux = blending * (2.0 / static_cast(_deltaZ)) * + (static_cast(_curVelocity) * upwindState // upwind flux for convection + - d_ax * derivative); // central FD approximation of c_x + + localStateDer[0] += flux * _invWeights[_polyDeg]; + localStateDer[_strideNode] -= flux * _invWeights[0]; + } + // outlet boundary BC + flux = blending * (2.0 / static_cast(_deltaZ)) * static_cast(_curVelocity) * (_dir == 1 ? _C[_nCells * (_strideNode * _nNodes) - _strideNode] : _C[0]); // upwind flux for convection + + if (_dir == 1) + stateDer[_nCells * (_strideNode * _nNodes) - _strideNode] += flux * _invWeights[_polyDeg]; + else + stateDer[0] += flux * _invWeights[0]; + + } + /** + * @brief calculates the subcell FV fluxes for convection + * @detail Only inner subcell FV fluxes are computed, DG-element fluxes need to be computed elsewhere. + * @param [in] blending pointer to DG element-wise blending coeffcients + * @param [in] _C pointer to state vector at component of interest + * @param [in] stateDer pointer to residual vector at component of interest + * @param [in] d_ax axial dispersion coefficient of current component + */ + template + void subcellFVconvectionIntegral(double* blending, const StateType* _C, ResidualType* stateDer) + { + const StateType* localC = _C; + StateType upwindState = 0.0; + ResidualType* localStateDer = stateDer; + ResidualType flux = 0.0; + + bool forward = (_dir == 1); // todo test backward flow + // inner subcell interfaces + for (unsigned int cell = 0; cell < _nCells; cell++, localC += _strideNode, localStateDer += _strideNode) { + + for (int interface = 1; interface < _nNodes; interface++, localC += _strideNode, localStateDer += _strideNode) { + + // upwind state for convection + if (cell == 0 && interface == 1) + upwindState = _subcellLimiter.reconstructedInterfaceValue(static_cast(_boundary[0]), localC[0], localC[_strideNode], interface - 1, forward); + else if (cell == _nCells - 1 && interface == _nNodes - 1) + upwindState = _subcellLimiter.reconstructedInterfaceValue(*(localC - _strideNode), localC[0], localC[0], interface - 1, forward); + else + upwindState = _subcellLimiter.reconstructedInterfaceValue(*(localC - _strideNode), localC[0], localC[_strideNode], interface - 1, forward); + + flux = blending[cell * _nComp] * (2.0 / static_cast(_deltaZ)) * static_cast(_curVelocity) * upwindState; + + // add flux divided by respective cell volume. Note that mapping was already applied + localStateDer[0] += flux * _invWeights[interface - 1]; + localStateDer[_strideNode] -= flux * _invWeights[interface]; + } + } + + /* DG element interfaces, i.e. subcell FV outer/element boundary faces */ + // Not needed when computed in DG surface integral. + + //// Inlet boundary BC + //flux = troubledCells[0] * (2.0 / static_cast(_deltaZ)) * static_cast(_curVelocity * _boundary[0]); + //if (_dir == 1) + // stateDer[0] -= flux * _invWeights[0]; + //else + // stateDer[(_nPoints - 1) * _strideNode] -= flux * _invWeights[_polyDeg]; + + //// Inner DG element faces + //localC = _C + _strideNode * (_nNodes - 1); + //localStateDer = stateDer + _strideNode * (_nNodes - 1); + + //for (unsigned int elemFace = 1; elemFace < _nCells; elemFace++, localC += _nNodes * _strideNode, localStateDer += _nNodes * _strideNode) { + + // upwindState = _subcellLimiter.reconstructedInterfaceValue(*(localC - _strideNode), localC[0], localC[_strideNode], _polyDeg, forward); + + // flux = troubledCells[elemFace-1] * (2.0 / static_cast(_deltaZ)) * static_cast(_curVelocity) * upwindState; + + // localStateDer[0] += flux * _invWeights[_polyDeg]; + // localStateDer[_strideNode] -= flux * _invWeights[0]; + //} + //// outlet boundary BC + //flux = troubledCells[_nCells-1] * (2.0 / static_cast(_deltaZ)) * static_cast(_curVelocity) * (_dir == 1 ? _C[_nCells * (_strideNode * _nNodes) - _strideNode] : _C[0]); + + //if (_dir == 1) + // stateDer[_nCells * (_strideNode * _nNodes) - _strideNode] += flux * _invWeights[_polyDeg]; + //else + // stateDer[0] += flux * _invWeights[0]; + + } + /** * @brief computes ghost nodes to implement boundary conditions * @detail to implement Danckwert boundary conditions, we only need to set the solid wall BC values for auxiliary variable @@ -939,7 +1089,7 @@ namespace cadet /** * @brief sets the sparsity pattern of the convection dispersion Jacobian for the nodal DG scheme */ - int ConvDispNodalPattern(std::vector& tripletList, const int offC = 0) { + int ConvDispCollocationDGPattern(std::vector& tripletList, const int offC = 0) { /*======================================================*/ /* Define Convection Jacobian Block */ @@ -1067,7 +1217,7 @@ namespace cadet /** * @brief sets the sparsity pattern of the convection dispersion Jacobian for the exact integration (her: modal) DG scheme */ - int ConvDispModalPattern(std::vector& tripletList, const int offC = 0) { + int ConvDispExIntDGPattern(std::vector& tripletList, const int offC = 0) { /*======================================================*/ /* Define Convection Jacobian Block */ @@ -1655,4 +1805,4 @@ namespace cadet } // namespace model } // namespace cadet -#endif // LIBCADET_CONVECTIONDISPERSIONOPERATORDG_HPP_ \ No newline at end of file +#endif // LIBCADET_CONVECTIONDISPERSIONOPERATORDG_HPP_