diff --git a/src/libcadet/model/GeneralRateModelDG.cpp b/src/libcadet/model/GeneralRateModelDG.cpp index 8287216dd..3a1368994 100644 --- a/src/libcadet/model/GeneralRateModelDG.cpp +++ b/src/libcadet/model/GeneralRateModelDG.cpp @@ -144,7 +144,8 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP LOG(Warning) << "Polynomial degree > 2 in bulk discretization (cf. POLYDEG) is always recommended for performance reasons."; _disc.nNodes = _disc.polyDeg + 1u; _disc.nPoints = _disc.nNodes * _disc.nCol; - _disc.exactInt = paramProvider.getBool("EXACT_INTEGRATION"); + const int polynomial_integration_mode = paramProvider.getInt("EXACT_INTEGRATION"); + _disc.exactInt = static_cast(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix const std::vector nParCell = paramProvider.getIntArray("NPARCELL"); const std::vector parPolyDeg = paramProvider.getIntArray("PARPOLYDEG"); @@ -446,7 +447,7 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP } unsigned int strideColNode = _disc.nComp; - const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, helper, _disc.nComp, _disc.exactInt, _disc.nCol, _disc.polyDeg, strideColNode); + const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, helper, _disc.nComp, polynomial_integration_mode, _disc.nCol, _disc.polyDeg, strideColNode); _disc.curSection = -1; diff --git a/src/libcadet/model/LumpedRateModelWithPoresDG.cpp b/src/libcadet/model/LumpedRateModelWithPoresDG.cpp index f76564116..9ed195197 100644 --- a/src/libcadet/model/LumpedRateModelWithPoresDG.cpp +++ b/src/libcadet/model/LumpedRateModelWithPoresDG.cpp @@ -114,7 +114,8 @@ bool LumpedRateModelWithPoresDG::configureModelDiscretization(IParameterProvider _disc.nNodes = _disc.polyDeg + 1; _disc.nPoints = _disc.nNodes * _disc.nCol; - _disc.exactInt = paramProvider.getBool("EXACT_INTEGRATION"); + const int polynomial_integration_mode = paramProvider.getInt("EXACT_INTEGRATION"); + _disc.exactInt = static_cast(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix const std::vector nBound = paramProvider.getIntArray("NBOUND"); if (nBound.size() < _disc.nComp) @@ -222,7 +223,7 @@ bool LumpedRateModelWithPoresDG::configureModelDiscretization(IParameterProvider paramProvider.popScope(); const unsigned int strideNode = _disc.nComp; - const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, helper, _disc.nComp, _disc.exactInt, _disc.nCol, _disc.polyDeg, strideNode); + const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, helper, _disc.nComp, polynomial_integration_mode, _disc.nCol, _disc.polyDeg, strideNode); _disc.curSection = -1; diff --git a/src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp b/src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp index 62bb8e358..1587c2b49 100644 --- a/src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp +++ b/src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp @@ -96,7 +96,8 @@ namespace cadet paramProvider.pushScope("discretization"); - _disc.exactInt = paramProvider.getBool("EXACT_INTEGRATION"); + const int polynomial_integration_mode = paramProvider.getInt("EXACT_INTEGRATION"); + _disc.exactInt = static_cast(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix if (paramProvider.getInt("POLYDEG") < 1) throw InvalidParameterException("Polynomial degree must be at least 1!"); @@ -143,7 +144,7 @@ namespace cadet paramProvider.popScope(); const unsigned int strideNode = _disc.nComp + _disc.strideBound; - const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, helper, _disc.nComp, _disc.exactInt, _disc.nCol, _disc.polyDeg, strideNode); + const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, helper, _disc.nComp, polynomial_integration_mode, _disc.nCol, _disc.polyDeg, strideNode); _disc.curSection = -1; diff --git a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp index ce7deabe0..ce0a6f629 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp @@ -70,10 +70,11 @@ AxialConvectionDispersionOperatorBaseDG::~AxialConvectionDispersionOperatorBaseD * @param [in] strideNode node stride in state vector * @return @c true if configuration went fine, @c false otherwise */ -bool AxialConvectionDispersionOperatorBaseDG::configureModelDiscretization(IParameterProvider& paramProvider, const IConfigHelper& helper, unsigned int nComp, bool exact_integration, unsigned int nCells, unsigned int polyDeg, unsigned int strideNode) +bool AxialConvectionDispersionOperatorBaseDG::configureModelDiscretization(IParameterProvider& paramProvider, const IConfigHelper& helper, unsigned int nComp, int polynomial_integration_mode, unsigned int nCells, unsigned int polyDeg, unsigned int strideNode) { + _nComp = nComp; - _exactInt = exact_integration; + _exactInt = static_cast(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix _nCells = nCells; _polyDeg = polyDeg; _nNodes = _polyDeg + 1u; @@ -105,6 +106,9 @@ bool AxialConvectionDispersionOperatorBaseDG::configureModelDiscretization(IPara invMMatrix(); derivativeMatrix(); + if(polynomial_integration_mode == 2) // use Gauss quadrature for exact integration + _invMM = gaussQuadratureMMatrix(_nodes, _nNodes).inverse(); + if (paramProvider.exists("COL_DISPERSION_DEP")) { const std::string paramDepName = paramProvider.getString("COL_DISPERSION_DEP"); @@ -718,9 +722,9 @@ int ConvectionDispersionOperatorDG::requiredADdirs() const CADET_NOEXC * @return @c true if configuration went fine, @c false otherwise */ template -bool ConvectionDispersionOperatorDG::configureModelDiscretization(IParameterProvider& paramProvider, const IConfigHelper& helper, unsigned int nComp, bool exact_integration, unsigned int nCells, unsigned int polyDeg, unsigned int strideNode) +bool ConvectionDispersionOperatorDG::configureModelDiscretization(IParameterProvider& paramProvider, const IConfigHelper& helper, unsigned int nComp, int polynomial_integration_mode, unsigned int nCells, unsigned int polyDeg, unsigned int strideNode) { - const bool retVal = _baseOp.configureModelDiscretization(paramProvider, helper, nComp, exact_integration, nCells, polyDeg, strideNode); + const bool retVal = _baseOp.configureModelDiscretization(paramProvider, helper, nComp, polynomial_integration_mode, nCells, polyDeg, strideNode); // todo: manage jacobians in convDispOp instead of unitOp ? //// Allocate memory diff --git a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp index 6dfc75204..f5a4a1b8a 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp @@ -83,7 +83,7 @@ namespace cadet void setFlowRates(const active& in, const active& out, const active& colPorosity) CADET_NOEXCEPT; - bool configureModelDiscretization(IParameterProvider& paramProvider, const IConfigHelper& helper, unsigned int nComp, bool exact_integration, unsigned int nCells, unsigned int polyDeg, unsigned int strideNode); + bool configureModelDiscretization(IParameterProvider& paramProvider, const IConfigHelper& helper, unsigned int nComp, int polynomial_integration_mode, unsigned int nCells, unsigned int polyDeg, unsigned int strideNode); bool configure(UnitOpIdx unitOpIdx, IParameterProvider& paramProvider, std::unordered_map& parameters); bool notifyDiscontinuousSectionTransition(double t, unsigned int secIdx, Eigen::MatrixXd& jacInlet); @@ -228,7 +228,6 @@ 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 @@ -281,6 +280,172 @@ namespace cadet _invWeights = _invWeights.cwiseInverse(); } + /** + * @brief computes the Legendre polynomial and its derivative + * @param [in, out] leg Legendre polynomial + * @param [in, out] leg Legendre polynomial derivative + * @param [in] N polynomial degree + * @param [in] x evaluation point + */ + void legendrePolynomialAndDerivative(double& leg, double& legDer, const int N, const double x) { + + switch (N) { + case 0: + leg = 1.0; + legDer = 0.0; + break; + case 1: + leg = x; + legDer = 1.0; + break; + default: + double leg_2 = 1.0; + double leg_1 = x; + double legDer_2 = 0.0; + double legDer_1 = 1.0; + + for (int k = 2; k <= N; k++) { + leg = (2.0 * k - 1.0) / k * x * leg_1 - (k - 1.0) / k * leg_2; + legDer = legDer_2 + (2.0 * k - 1.0) * leg_1; + leg_2 = leg_1; + leg_1 = leg; + legDer_2 = legDer_1; + legDer_1 = legDer; + } + } + } + + /** + * @brief computes the Legendre-Gauss nodes and quadrature weights + * @param [in, out] nodes Legendre Gauss nodes + * @param [in, out] weights Legendre Gauss quadrature weights + * @param [in] N polynomial degree + */ + void lgNodesWeights(Eigen::VectorXd& nodes, Eigen::VectorXd& weights, const int N) { + // tolerance and max #iterations for Newton iteration + int nIterations = 10; + double tolerance = 1e-15; + + switch (N) { + case 0: + nodes[0] = 0.0; + weights[0] = 2.0; + break; + case 1: + nodes[0] = -std::sqrt(1.0 / 3.0); + weights[0] = 1; + nodes[1] = -nodes[0]; + weights[1] = weights[0]; + break; + default: + + double leg = 0.0; + double legDer = 0.0; + double delta = 0.0; + + for (int j = 0; j <= std::floor((N + 1) / 2) - 1; j++) + { + nodes[j] = -std::cos((2.0 * j + 1.0) / (2.0 * N + 2.0) * M_PI); + for (int k = 0; k <= nIterations; k++) + { + legendrePolynomialAndDerivative(leg, legDer, N + 1, nodes[j]); + delta = -leg / legDer; + nodes[j] = nodes[j] + delta; + if (std::abs(delta) <= tolerance * std::abs(nodes[j])) + break; + } + legendrePolynomialAndDerivative(leg, legDer, N + 1, nodes[j]); + nodes[N - j] = -nodes[j]; + weights[j] = 2.0 / ((1.0 - std::pow(nodes[j], 2.0)) * std::pow(legDer, 2.0)); + weights[N - j] = weights[j]; + } + + if (N % 2 == 0) + { + legendrePolynomialAndDerivative(leg, legDer, N + 1, 0.0); + nodes[N / 2] = 0.0; + weights[N / 2] = 2.0 / std::pow(legDer, 2.0); + } + } + } + + /** + * @brief evaluates a Lagrange polynomial built on input nodes at a set of points + * @detail can be used to establish quadrature rules + * @param [in] j index of Lagrange basis function + * @param [in] intNodes interpolation nodes the Lagrange basis is constructed with + * @param [in] evalNodes nodes the Lagrange basis is evaluated at + */ + Eigen::VectorXd evalLagrangeBasis(const int j, const Eigen::VectorXd intNodes, const Eigen::VectorXd evalNodes) { + + const int nIntNodes = intNodes.size(); + const int nEvalNodes = evalNodes.size(); + VectorXd evalEll = VectorXd::Zero(nEvalNodes); + + double nominator = 1.0; + double denominator = 1.0; + + for (int i = 0; i < nIntNodes; i++) + if (i != j) + denominator *= (intNodes[j] - intNodes[i]); + + for (int k = 0; k < nEvalNodes; k++) + { + for (int i = 0; i < nIntNodes; i++) + { + if (i != j) + { + if (std::abs(evalNodes[k] - intNodes[i]) < std::numeric_limits::epsilon()) + { + nominator = denominator; + break; + } + else + nominator *= (evalNodes[k] - intNodes[i]); + } + } + evalEll[k] = nominator / denominator; + nominator = 1.0; + } + + return evalEll; + } + + /** + * @brief calculates the Gauss quadrature mass matrix for LGL interpolation polynomial on LG points + * @detail exact integration of polynomials up to order 2N - 1 + * @param [in] LGLnodes Lagrange Gauss Lobatto nodes + * @param [in] nGaussNodes number of Gauss quadrature nodes + */ + Eigen::MatrixXd gaussQuadratureMMatrix(const Eigen::VectorXd LGLnodes, const int nLGNodes) { + + const int Ldegree = nLGNodes - 1; // Legendre polynomial degree + const int nLGLnodes = LGLnodes.size(); + + MatrixXd evalEll = MatrixXd::Zero(nLGNodes, nLGNodes); + MatrixXd massMatrix = MatrixXd::Zero(nLGNodes, nLGNodes); + + VectorXd LGnodes = VectorXd::Zero(nLGNodes); + VectorXd LGweigths = VectorXd::Zero(nLGNodes); + lgNodesWeights(LGnodes, LGweigths, Ldegree); + + for (int i = 0; i < nLGLnodes; i++) + evalEll.row(i) = evalLagrangeBasis(i, LGLnodes, LGnodes); + + VectorXd aux = VectorXd::Zero(nLGNodes); + + for (int i = 0; i < nLGLnodes; i++) + { + for (int j = 0; j < nLGLnodes; j++) + { + aux = evalEll.row(i).array() * evalEll.row(j).array(); + massMatrix(i, j) = (aux.array() * LGweigths.array()).sum(); + } + } + + return massMatrix; + } + /** * @brief computation of barycentric weights for fast polynomial evaluation * @param [in] baryWeights vector to store barycentric weights. Must already be initialized with ones! @@ -362,8 +527,8 @@ namespace cadet return V; } /** - * @brief calculates mass matrix for exact polynomial integration - * @detail exact polynomial integration leads to a full mass matrix + * @brief calculates mass matrix via transformation to orthonormal Legendre (modal) basis + * @detail exact integration for integrals of the form \int_E \ell_i \ell_j d\xi */ void invMMatrix() { _invMM = (getVandermonde_LEGENDRE() * (getVandermonde_LEGENDRE().transpose())); @@ -1392,7 +1557,7 @@ namespace cadet void setFlowRates(const active& in, const active& out, const active& colPorosity) CADET_NOEXCEPT; - bool configureModelDiscretization(IParameterProvider& paramProvider, const IConfigHelper& helper, unsigned int nComp, bool exact_integration, unsigned int nCells, unsigned int polyDeg, unsigned int strideNode); + bool configureModelDiscretization(IParameterProvider& paramProvider, const IConfigHelper& helper, unsigned int nComp, int polynomial_integration_mode, unsigned int nCells, unsigned int polyDeg, unsigned int strideNode); bool configure(UnitOpIdx unitOpIdx, IParameterProvider& paramProvider, std::unordered_map& parameters); bool notifyDiscontinuousSectionTransition(double t, unsigned int secIdx, const AdJacobianParams& adJac);