Skip to content

Commit

Permalink
Add Gauss quadrature for DG
Browse files Browse the repository at this point in the history
  • Loading branch information
@jbreue16
jbreue16 committed Jun 18, 2024
1 parent 2925af5 commit b54811e
Show file tree
Hide file tree
Showing 5 changed files with 187 additions and 15 deletions.
5 changes: 3 additions & 2 deletions src/libcadet/model/GeneralRateModelDG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<bool>(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix

const std::vector<int> nParCell = paramProvider.getIntArray("NPARCELL");
const std::vector<int> parPolyDeg = paramProvider.getIntArray("PARPOLYDEG");
Expand Down Expand Up @@ -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;

Expand Down
5 changes: 3 additions & 2 deletions src/libcadet/model/LumpedRateModelWithPoresDG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<bool>(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix

const std::vector<int> nBound = paramProvider.getIntArray("NBOUND");
if (nBound.size() < _disc.nComp)
Expand Down Expand Up @@ -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;

Expand Down
5 changes: 3 additions & 2 deletions src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<bool>(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!");
Expand Down Expand Up @@ -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;

Expand Down
12 changes: 8 additions & 4 deletions src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<bool>(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix
_nCells = nCells;
_polyDeg = polyDeg;
_nNodes = _polyDeg + 1u;
Expand Down Expand Up @@ -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");
Expand Down Expand Up @@ -718,9 +722,9 @@ int ConvectionDispersionOperatorDG<Operator>::requiredADdirs() const CADET_NOEXC
* @return @c true if configuration went fine, @c false otherwise
*/
template <typename Operator>
bool ConvectionDispersionOperatorDG<Operator>::configureModelDiscretization(IParameterProvider& paramProvider, const IConfigHelper& helper, unsigned int nComp, bool exact_integration, unsigned int nCells, unsigned int polyDeg, unsigned int strideNode)
bool ConvectionDispersionOperatorDG<Operator>::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
Expand Down
175 changes: 170 additions & 5 deletions src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<ParameterId, active*>& parameters);
bool notifyDiscontinuousSectionTransition(double t, unsigned int secIdx, Eigen::MatrixXd& jacInlet);

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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<double>::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!
Expand Down Expand Up @@ -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()));
Expand Down Expand Up @@ -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<ParameterId, active*>& parameters);
bool notifyDiscontinuousSectionTransition(double t, unsigned int secIdx, const AdJacobianParams& adJac);

Expand Down

0 comments on commit b54811e

Please sign in to comment.