From c3ff5dbdd1fb1f9c2c2b4493b3f84471b5d4c7c8 Mon Sep 17 00:00:00 2001 From: Jan Breuer <74359367+jbreue16@users.noreply.github.com> Date: Tue, 27 Feb 2024 16:05:51 +0100 Subject: [PATCH] Add particle GSM discretization (#173) --- src/libcadet/model/GeneralRateModelDG.cpp | 289 ++++++++++++++-------- src/libcadet/model/GeneralRateModelDG.hpp | 79 ++++-- src/libcadet/model/parts/DGToolbox.cpp | 57 ++++- src/libcadet/model/parts/DGToolbox.hpp | 20 +- 4 files changed, 317 insertions(+), 128 deletions(-) diff --git a/src/libcadet/model/GeneralRateModelDG.cpp b/src/libcadet/model/GeneralRateModelDG.cpp index 496fa83dd..c90252d3e 100644 --- a/src/libcadet/model/GeneralRateModelDG.cpp +++ b/src/libcadet/model/GeneralRateModelDG.cpp @@ -191,6 +191,7 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP _disc.parPolyDeg = new unsigned int[_disc.nParType]; _disc.nParCell = new unsigned int[_disc.nParType]; _disc.parExactInt = new bool[_disc.nParType]; + _disc.parGSM = new bool[_disc.nParType]; if (paramProvider.exists("PAR_POLYDEG")) { @@ -247,6 +248,23 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP else throw InvalidParameterException("Specify field PAR_POLYDEG (or NPAR)"); + if (paramProvider.exists("PAR_GSM")) + { + std::vector parGSM = paramProvider.getBoolArray("PAR_GSM"); + if (parGSM.size() == 1) + std::fill(_disc.parGSM, _disc.parGSM + _disc.nParType, parGSM[0]); + else + std::copy_n(parGSM.begin(), _disc.nParType, _disc.parGSM); + for (int type = 0; type < _disc.nParType; type++) + if (_disc.parGSM[type] && _disc.nParCell[type] != 1) + throw InvalidParameterException("Field PAR_NELEM must equal one to use a GSM discretization in the corresponding particle type"); + } + else // Use GSM as default for particle discretization + { + for (int type = 0; type < _disc.nParType; type++) + _disc.parGSM[type] = (_disc.nParCell[type] == 1); + } + // Compute discretization operators and initialize containers _disc.initializeDG(); @@ -590,7 +608,7 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP _globalJacDisc.resize(numDofs(), numDofs()); _globalJac.resize(numDofs(), numDofs()); // pattern is set in configure(), after surface diffusion is read - //FDJac = MatrixXd::Zero(numDofs(), numDofs()); // todo delete + // FDJac = MatrixXd::Zero(numDofs(), numDofs()); // todo delete // Set whether analytic Jacobian is used useAnalyticJacobian(analyticJac); @@ -1302,7 +1320,7 @@ int GeneralRateModelDG::residualParticle(double t, unsigned int parType, unsigne LinearBufferAllocator tlmAlloc = threadLocalMem.get(); - // special case: individual treatment of time derivatives in particle mass balance at inner particle boundary node + // Special case: individual treatment of time derivatives in particle mass balance at inner particle boundary node bool specialCase = !_disc.parExactInt[parType] && (_parGeomSurfToVol[parType] != _disc.SurfVolRatioSlab && _parCoreRadius[parType] == 0.0); // Prepare parameters @@ -1310,7 +1328,7 @@ int GeneralRateModelDG::residualParticle(double t, unsigned int parType, unsigne // Ordering of particle surface diffusion: // bnd0comp0, bnd0comp1, bnd0comp2, bnd1comp0, bnd1comp1, bnd1comp2 - active const* const _parSurfDiff = getSectionDependentSlice(_parSurfDiffusion, _disc.strideBound[_disc.nParType], secIdx) + _disc.nBoundBeforeType[parType]; + active const* const parSurfDiff = getSectionDependentSlice(_parSurfDiffusion, _disc.strideBound[_disc.nParType], secIdx) + _disc.nBoundBeforeType[parType]; // z coordinate (column length normed to 1) of current node - needed in externally dependent adsorption kinetic const double z = _convDispOp.relativeCoordinate(colNode); @@ -1344,11 +1362,11 @@ int GeneralRateModelDG::residualParticle(double t, unsigned int parType, unsigne const ColumnPosition colPos{ z, 0.0, r }; // Handle time derivatives, binding, dynamic reactions. - // if special case: Dont add time derivatives to inner boundary node for DG discretized mass balance equations. + // Special case: Dont add time derivatives to inner boundary node for DG discretized mass balance equations. // This can be achieved by setting yDot pointer to null before passing to residual kernel, and adding only the time derivative for dynamic binding // TODO Check Treatment of reactions (do we need yDot then?) - if (cadet_unlikely(par == 0 && specialCase)) { - + if (cadet_unlikely(par == 0 && specialCase)) + { parts::cell::residualKernel( t, secIdx, colPos, local_y, nullptr, local_res, jac, cellResParams, tlmAlloc // TODO Check Treatment of reactions (do we need yDot then?) ); @@ -1364,12 +1382,12 @@ int GeneralRateModelDG::residualParticle(double t, unsigned int parType, unsigne if (cellResParams.qsReaction[idx]) continue; - // for kinetic bindings and surface diffusion, we have an additional DG-discretized mass balance eq. + // For kinetic bindings and surface diffusion, we have an additional DG-discretized mass balance eq. // -> add time derivate at inner bonudary node only without surface diffusion else if (_hasSurfaceDiffusion[parType]) continue; - // some bound states might still not be effected by surface diffusion - else if (_parSurfDiff[idx] != 0.0) + // Some bound states might still not be effected by surface diffusion + else if (parSurfDiff[idx] != 0.0) continue; // Add time derivative to solid phase @@ -1378,27 +1396,26 @@ int GeneralRateModelDG::residualParticle(double t, unsigned int parType, unsigne } } } - else { - + else + { parts::cell::residualKernel( t, secIdx, colPos, local_y, local_yDot, local_res, jac, cellResParams, tlmAlloc ); - } - // move rowiterator to next particle node + // Move rowiterator to next particle node jac += idxr.strideParNode(parType); } // We still need to handle transport/diffusion - // get pointers to the particle block of the current column node, particle type + // Get pointers to the particle block of the current column node, particle type const StateType* c_p = yBase + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }); ResidualType* resC_p = resBase + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }); // Mobile phase RHS - // get film diffusion flux at current node to compute boundary condition + // Get film diffusion flux at current node to compute boundary condition active const* const filmDiff = getSectionDependentSlice(_filmDiffusion, _disc.nComp * _disc.nParType, secIdx) + parType * _disc.nComp; for (unsigned int comp = 0; comp < _disc.nComp; comp++) { _disc.localFlux[comp] = filmDiff[comp] * (yBase[idxr.offsetC() + colNode * idxr.strideColNode() + comp] @@ -1413,98 +1430,159 @@ int GeneralRateModelDG::residualParticle(double t, unsigned int parType, unsigne int strideParLiquid = idxr.strideParLiquid(); int strideParNode = idxr.strideParNode(parType); - for (unsigned int comp = 0; comp < nComp; comp++) + if (_disc.parGSM[parType]) // GSM implementation { - // component dependent (through access factor) inverse Beta_P - ParamType invBetaP = (1.0 - static_cast(_parPorosity[parType])) / (static_cast(_poreAccessFactor[_disc.nComp * parType + comp]) * static_cast(_parPorosity[parType])); - - // =====================================================================================================// - // solve auxiliary systems d_p g_p + d_s beta_p sum g_s= d (d_p c_p + d_s beta_p sum c_s) / d xi // - // =====================================================================================================// - // component-wise! strides - unsigned int strideCell = nNodes; - unsigned int strideNode = 1u; - // reset cache for auxiliary variable - Eigen::Map, 0, InnerStride<>> _g_p(reinterpret_cast(&_disc.g_p[parType][0]), nPoints, InnerStride<>(1)); - Eigen::Map, 0, InnerStride<>> _g_pSum(reinterpret_cast(&_disc.g_pSum[parType][0]), nPoints, InnerStride<>(1)); - _g_p.setZero(); - _g_pSum.setZero(); - - Eigen::Map, 0, InnerStride> cp(c_p + comp, _disc.nParPoints[parType], InnerStride(idxr.strideParNode(parType))); + for (unsigned int comp = 0; comp < nComp; comp++) + { - // handle surface diffusion: Compute auxiliary variable; For kinetic bindings: add additional mass balance to residual of respective bound state - if (_hasSurfaceDiffusion[parType]) { + // ====================================================================================// + // solve GSM-discretized particle mass balance // + // ====================================================================================// - for (int bnd = 0; bnd < _disc.nBound[parType * _disc.nComp + comp]; bnd++) { + Eigen::Map, 0, InnerStride> resCp(resC_p + comp, nPoints, InnerStride(idxr.strideParNode(parType))); + Eigen::Map, 0, InnerStride> Cp(c_p + comp, nPoints, InnerStride(idxr.strideParNode(parType))); - if (_parSurfDiff[getOffsetSurfDiff(parType, comp, bnd)] != 0.0) { // some bound states might still not be effected by surface diffusion + // Use auxiliary variable to get c^p + \sum 1 / \Beta_p c^s + Eigen::Map, 0, InnerStride<>> sum_cp_cs(reinterpret_cast(&_disc.g_pSum[parType][0]), nPoints, InnerStride<>(1)); + sum_cp_cs = static_cast(parDiff[comp]) * Cp.template cast(); - // get solid phase vector - Eigen::Map, 0, InnerStride> q_p(c_p + strideParLiquid + idxr.offsetBoundComp(ParticleTypeIndex{ parType }, ComponentIndex{ comp }) + bnd, - nPoints, InnerStride(strideParNode)); - // compute g_s = d c_s / d xi - solve_auxiliary_DG(parType, q_p, strideCell, strideNode, comp); - // apply invBeta_p, d_s and add to sum -> gSum += d_s * invBeta_p * (D c - M^-1 B [c - c^*]) - _g_pSum += _g_p.template cast() * invBetaP * static_cast(_parSurfDiff[getOffsetSurfDiff(parType, comp, bnd)]); + for (int bnd = 0; bnd < _disc.nBound[parType * _disc.nComp + comp]; bnd++) + { + if (parSurfDiff[getOffsetSurfDiff(parType, comp, bnd)] != 0.0) // some bound states might still not be effected by surface diffusion + { + Eigen::Map, 0, InnerStride> c_s(c_p + _disc.nComp + idxr.offsetBoundComp(ParticleTypeIndex{ parType }, ComponentIndex{ comp }) + bnd, nPoints, InnerStride(idxr.strideParNode(parType))); + ParamType invBetaP = (1.0 - static_cast(_parPorosity[parType])) / (static_cast(_poreAccessFactor[_disc.nComp * parType + comp]) * static_cast(_parPorosity[parType])); + sum_cp_cs += invBetaP * static_cast(parSurfDiff[getOffsetSurfDiff(parType, comp, bnd)]) * c_s; /* For kinetic bindings with surface diffusion: add the additional DG-discretized particle mass balance equations to residual */ - if (!qsReaction[bnd]) { - + if (!qsReaction[bnd]) + { // Eigen access to current bound state residual Eigen::Map, 0, InnerStride> resCs(resBase + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) + idxr.strideParLiquid() + idxr.offsetBoundComp(ParticleTypeIndex{ parType }, ComponentIndex{ comp }) + bnd, nPoints, InnerStride(idxr.strideParNode(parType))); - // promote auxiliary variable storage from double to active if required - // @todo is there a more efficient or elegant solution? - if (std::is_same::value && std::is_same::value) - vectorPromoter(reinterpret_cast(&_g_p[0]), nPoints); // reinterpret_cast only required because statement is scanned also when StateType != double + // Use auxiliary variable to get \Beta_p D_s c^s + Eigen::Map, 0, InnerStride<>> c_s_modified(reinterpret_cast(&_disc.g_p[parType][0]), nPoints, InnerStride<>(1)); - // Access auxiliary variable as ResidualType - Eigen::Map, 0, InnerStride<>> _g_p_ResType(reinterpret_cast(&_disc.g_p[parType][0]), nPoints, InnerStride<>(1)); + // Apply squared inverse mapping and surface diffusion + c_s_modified = 2.0 / static_cast(_disc.deltaR[_disc.offsetMetric[parType]]) * 2.0 / static_cast(_disc.deltaR[_disc.offsetMetric[parType]]) * + static_cast(parSurfDiff[getOffsetSurfDiff(parType, comp, bnd)]) * c_s; - applyParInvMap(_g_p_ResType, parType); - _g_p_ResType *= static_cast(_parSurfDiff[getOffsetSurfDiff(parType, comp, bnd)]); + Eigen::Map, 0, InnerStride> c_s_modified_const(&c_s_modified[0], nPoints, InnerStride(1)); + parGSMVolumeIntegral(parType, c_s_modified_const, resCs); - // Eigen access to auxiliary variable of current bound state - Eigen::Map, 0, InnerStride> _g_p_ResType_const(&_g_p_ResType[0], nPoints, InnerStride(1)); + // Leave out the surface integral as we only have one element, i.e. we apply BC with zeros + } + } + } + + // Apply squared inverse mapping + sum_cp_cs *= 2.0 / static_cast(_disc.deltaR[_disc.offsetMetric[parType]]) * 2.0 / static_cast(_disc.deltaR[_disc.offsetMetric[parType]]); - // adds - D_r * gs to the residual, including metric part.->res = invMap^2* [ -D_r * (d_s c^s) ] - parVolumeIntegral(parType, false, _g_p_ResType_const, resCs); + Eigen::Map, 0, InnerStride> sum_cp_cs_const(&sum_cp_cs[0], nPoints, InnerStride(1)); + parGSMVolumeIntegral(parType, sum_cp_cs_const, resCp); - // adds M^-1 B (gs - gs^*) to the residual -> res = invMap^2 * [ - D_r * (d_s c^s) + M^-1 B (gs - gs^*) ] - parSurfaceIntegral(parType, _g_p_ResType_const, resCs, strideCell, strideNode, false, comp, true); + // Pass sum_cp_cs_const to match the DGSEM interface; nullptr might also be feasible + parSurfaceIntegral(parType, sum_cp_cs_const, resCp, nNodes, 1u, false, comp); + + } + } + else // DGSEM implementation + { + for (unsigned int comp = 0; comp < nComp; comp++) + { + // Component dependent (through access factor) inverse Beta_P + ParamType invBetaP = (1.0 - static_cast(_parPorosity[parType])) / (static_cast(_poreAccessFactor[_disc.nComp * parType + comp]) * static_cast(_parPorosity[parType])); + + // =====================================================================================================// + // Solve auxiliary systems d_p g_p + d_s beta_p sum g_s= d (d_p c_p + d_s beta_p sum c_s) / d xi // + // =====================================================================================================// + // Component-wise! strides + unsigned int strideCell = nNodes; + unsigned int strideNode = 1u; + // Reset cache for auxiliary variable + Eigen::Map, 0, InnerStride<>> _g_p(reinterpret_cast(&_disc.g_p[parType][0]), nPoints, InnerStride<>(1)); + Eigen::Map, 0, InnerStride<>> _g_pSum(reinterpret_cast(&_disc.g_pSum[parType][0]), nPoints, InnerStride<>(1)); + _g_p.setZero(); + _g_pSum.setZero(); + + Eigen::Map, 0, InnerStride> cp(c_p + comp, _disc.nParPoints[parType], InnerStride(idxr.strideParNode(parType))); + + // Handle surface diffusion: Compute auxiliary variable; For kinetic bindings: add additional mass balance to residual of respective bound state + if (_hasSurfaceDiffusion[parType]) + { + for (int bnd = 0; bnd < _disc.nBound[parType * _disc.nComp + comp]; bnd++) + { + if (parSurfDiff[getOffsetSurfDiff(parType, comp, bnd)] != 0.0) // some bound states might still not be effected by surface diffusion + { + // Get solid phase vector + Eigen::Map, 0, InnerStride> c_s(c_p + strideParLiquid + idxr.offsetBoundComp(ParticleTypeIndex{ parType }, ComponentIndex{ comp }) + bnd, + nPoints, InnerStride(strideParNode)); + // Compute g_s = d c_s / d xi + solve_auxiliary_DG(parType, c_s, strideCell, strideNode, comp); + // Apply invBeta_p, d_s and add to sum -> gSum += d_s * invBeta_p * (D c - M^-1 B [c - c^*]) + _g_pSum += _g_p.template cast() * invBetaP * static_cast(parSurfDiff[getOffsetSurfDiff(parType, comp, bnd)]); + + /* For kinetic bindings with surface diffusion: add the additional DG-discretized particle mass balance equations to residual */ + + if (!qsReaction[bnd]) + { + // Eigen access to current bound state residual + Eigen::Map, 0, InnerStride> resCs(resBase + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) + idxr.strideParLiquid() + idxr.offsetBoundComp(ParticleTypeIndex{ parType }, ComponentIndex{ comp }) + bnd, + nPoints, InnerStride(idxr.strideParNode(parType))); + + // Promote auxiliary variable storage from double to active if required + // @todo is there a more efficient or elegant solution? + if (std::is_same::value && std::is_same::value) + vectorPromoter(reinterpret_cast(&_g_p[0]), nPoints); // reinterpret_cast only required because statement is scanned also when StateType != double + + // Access auxiliary variable as ResidualType + Eigen::Map, 0, InnerStride<>> _g_p_ResType(reinterpret_cast(&_disc.g_p[parType][0]), nPoints, InnerStride<>(1)); + + applyParInvMap(_g_p_ResType, parType); + _g_p_ResType *= static_cast(parSurfDiff[getOffsetSurfDiff(parType, comp, bnd)]); + + // Eigen access to auxiliary variable of current bound state + Eigen::Map, 0, InnerStride> _g_p_ResType_const(&_g_p_ResType[0], nPoints, InnerStride(1)); + + // Add - D_r * gs to the residual, including metric part.->res = invMap^2* [ -D_r * (d_s c^s) ] + parVolumeIntegral(parType, false, _g_p_ResType_const, resCs); + + // Add M^-1 B (gs - gs^*) to the residual -> res = invMap^2 * [ - D_r * (d_s c^s) + M^-1 B (gs - gs^*) ] + parSurfaceIntegral(parType, _g_p_ResType_const, resCs, strideCell, strideNode, false, comp, true); + } } } } - } - // compute g_p = d c_p / d xi - solve_auxiliary_DG(parType, cp, strideCell, strideNode, comp); + // Compute g_p = d c_p / d xi + solve_auxiliary_DG(parType, cp, strideCell, strideNode, comp); + + // Add particle diffusion part to auxiliary variable sum -> gSum += d_p * (D c - M^-1 B [c - c^*]) + _g_pSum += _g_p * static_cast(parDiff[comp]); - // add particle diffusion part to auxiliary variable sum -> gSum += d_p * (D c - M^-1 B [c - c^*]) - _g_pSum += _g_p * static_cast(parDiff[comp]); + // apply squared inverse mapping to sum of bound state auxiliary variables -> gSum = - invMap^2 * (d_p * c^p + sum_mi d_s invBeta_p c^s) + applyParInvMap(_g_pSum, parType); - // apply squared inverse mapping to sum of bound state auxiliary variables -> gSum = - invMap^2 * (d_p * c^p + sum_mi d_s invBeta_p c^s) - applyParInvMap(_g_pSum, parType); + // ====================================================================================// + // solve DG-discretized particle mass balance // + // ====================================================================================// - // ====================================================================================// - // solve DG-discretized particle mass balance // - // ====================================================================================// + /* Solve DG-discretized particle mass balance equation */ - /* solve DG-discretized particle mass balance equation */ - - Eigen::Map, 0, InnerStride> _g_pSum_const(&_g_pSum[0], nPoints, InnerStride(1)); + Eigen::Map, 0, InnerStride> _g_pSum_const(&_g_pSum[0], nPoints, InnerStride(1)); - // Eigen access to particle liquid residual - Eigen::Map, 0, InnerStride> resCp(resC_p + comp, nPoints, InnerStride(idxr.strideParNode(parType))); + // Eigen access to particle liquid residual + Eigen::Map, 0, InnerStride> resCp(resC_p + comp, nPoints, InnerStride(idxr.strideParNode(parType))); - // adds - D_r * (g_sum) to the residual, including metric part. -> res = - D_r * (d_p * c^p + invBeta_p sum_mi d_s c^s) - parVolumeIntegral(parType, false, _g_pSum_const, resCp); + // Add - D_r * (g_sum) to the residual, including metric part. -> res = - D_r * (d_p * c^p + invBeta_p sum_mi d_s c^s) + parVolumeIntegral(parType, false, _g_pSum_const, resCp); - // adds M^-1 B (g_sum - g_sum^*) to the residual -> res = - D_r * (d_p * c^p + invBeta_p sum_mi d_s c^s) + M^-1 B (g_sum - g_sum^*) - parSurfaceIntegral(parType, _g_pSum_const, resCp, strideCell, strideNode, false, comp); + // Add M^-1 B (g_sum - g_sum^*) to the residual -> res = - D_r * (d_p * c^p + invBeta_p sum_mi d_s c^s) + M^-1 B (g_sum - g_sum^*) + parSurfaceIntegral(parType, _g_pSum_const, resCp, strideCell, strideNode, false, comp); + } } return 0; @@ -1991,8 +2069,10 @@ void GeneralRateModelDG::updateRadialDisc() for (unsigned int parType = 0; parType < _disc.nParType; ++parType) { - if (_parDiscType[parType] == ParticleDiscretizationMode::Equidistant) { - for (int cell = 0; cell < _disc.nParCell[parType]; cell++) { + if (_parDiscType[parType] == ParticleDiscretizationMode::Equidistant) + { + for (int cell = 0; cell < _disc.nParCell[parType]; cell++) + { _disc.deltaR[_disc.offsetMetric[parType] + cell] = (_parRadius[parType] - _parCoreRadius[parType]) / _disc.nParCell[parType]; } setEquidistantRadialDisc(parType); @@ -2006,15 +2086,14 @@ void GeneralRateModelDG::updateRadialDisc() /* metrics */ // estimate cell dependent D_r - for (int parType = 0; parType < _disc.nParType; parType++) { - - for (int cell = 0; cell < _disc.nParCell[parType]; cell++) { - + for (int parType = 0; parType < _disc.nParType; parType++) + { + for (int cell = 0; cell < _disc.nParCell[parType]; cell++) + { _disc.Ir[_disc.offsetMetric[parType] + cell] = Vector::Zero(_disc.nParNode[parType]); - for (int node = 0; node < _disc.nParNode[parType]; node++) { + for (int node = 0; node < _disc.nParNode[parType]; node++) _disc.Ir[_disc.offsetMetric[parType] + cell][node] = _disc.deltaR[_disc.offsetMetric[parType] + cell] / 2.0 * (_disc.parNodes[parType][node] + 1.0); - } _disc.Dr[_disc.offsetMetric[parType] + cell].resize(_disc.nParNode[parType], _disc.nParNode[parType]); _disc.Dr[_disc.offsetMetric[parType] + cell].setZero(); @@ -2034,30 +2113,44 @@ void GeneralRateModelDG::updateRadialDisc() _disc.Dr[_disc.offsetMetric[parType] + cell].array().colwise() *= _disc.Ir[_disc.offsetMetric[parType] + cell].array().template cast().cwiseInverse(); // compute mass matrices for exact integration based on particle geometry, via transformation to normalized Jacobi polynomials with weight function w - if (_parGeomSurfToVol[parType] == SurfVolRatioSphere) { // w = (1 + \xi)^2 - - _disc.parInvMM[_disc.offsetMetric[parType] + cell] = parts::dgtoolbox::invMMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 2.0).inverse() * pow((static_cast(_disc.deltaR[_disc.offsetMetric[parType] + cell]) / 2.0), 2.0); - if(cell > 0 || _parCoreRadius[parType] != 0.0) // following contributions are zero for first cell when R_c = 0 (no particle core) - _disc.parInvMM[_disc.offsetMetric[parType] + cell] += parts::dgtoolbox::invMMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 1.0).inverse() * (static_cast(_disc.deltaR[_disc.offsetMetric[parType] + cell]) * static_cast(r_L)) - + parts::dgtoolbox::invMMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 0.0).inverse() * pow(static_cast(r_L), 2.0); + if (_parGeomSurfToVol[parType] == SurfVolRatioSphere) // r^2 = r_i^2 + (1 + \xi) * r_i * DeltaR_i / 2.0 + (1 + \xi)^2 * (DeltaR_i / 2.0)^2 + { + _disc.parInvMM[_disc.offsetMetric[parType] + cell] = parts::dgtoolbox::mMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 2.0) * pow((static_cast(_disc.deltaR[_disc.offsetMetric[parType] + cell]) / 2.0), 2.0); + if (cell > 0 || _parCoreRadius[parType] != 0.0) // following contributions are zero for first cell when R_c = 0 (no particle core) + _disc.parInvMM[_disc.offsetMetric[parType] + cell] += parts::dgtoolbox::mMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 1.0) * (static_cast(_disc.deltaR[_disc.offsetMetric[parType] + cell]) * static_cast(r_L)) + + parts::dgtoolbox::mMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 0.0) * pow(static_cast(r_L), 2.0); _disc.parInvMM[_disc.offsetMetric[parType] + cell] = _disc.parInvMM[_disc.offsetMetric[parType] + cell].inverse(); - _disc.minus_InvMM_ST[_disc.offsetMetric[parType] + cell] = - _disc.parInvMM[_disc.offsetMetric[parType] + cell] * _disc.parPolyDerM[parType].transpose() * _disc.parInvMM[_disc.offsetMetric[parType] + cell].inverse(); - } - else if (_parGeomSurfToVol[parType] == SurfVolRatioCylinder) { // w = (1 + \xi) + _disc.minus_InvMM_ST[_disc.offsetMetric[parType] + cell] = -_disc.parInvMM[_disc.offsetMetric[parType] + cell] * _disc.parPolyDerM[parType].transpose() * _disc.parInvMM[_disc.offsetMetric[parType] + cell].inverse(); - _disc.parInvMM[_disc.offsetMetric[parType] + cell] = parts::dgtoolbox::invMMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 1.0).inverse() * (static_cast(_disc.deltaR[_disc.offsetMetric[parType] + cell]) / 2.0); + // particle GSM specific second order stiffness matrix (single element, i.e. nParCell = 1) + _disc.secondOrderStiffnessM[parType] = std::pow(static_cast(_parCoreRadius[parType]), 2.0) * parts::dgtoolbox::secondOrderStiffnessMatrix(_disc.parPolyDeg[parType], 0.0, 0.0, _disc.parNodes[parType]); + _disc.secondOrderStiffnessM[parType] += static_cast(_disc.deltaR[_disc.offsetMetric[parType]]) * static_cast(_parCoreRadius[parType]) * parts::dgtoolbox::secondOrderStiffnessMatrix(_disc.parPolyDeg[parType], 0.0, 1.0, _disc.parNodes[parType]); + _disc.secondOrderStiffnessM[parType] += std::pow(static_cast(_disc.deltaR[_disc.offsetMetric[parType]]) / 2.0, 2.0) * parts::dgtoolbox::secondOrderStiffnessMatrix(_disc.parPolyDeg[parType], 0.0, 2.0, _disc.parNodes[parType]); + } + else if (_parGeomSurfToVol[parType] == SurfVolRatioCylinder) // r = r_i + (1 + \xi) * DeltaR_i / 2.0 + { + _disc.parInvMM[_disc.offsetMetric[parType] + cell] = parts::dgtoolbox::mMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 1.0) * (static_cast(_disc.deltaR[_disc.offsetMetric[parType] + cell]) / 2.0); if (cell > 0 || _parCoreRadius[parType] != 0.0) // following contribution is zero for first cell when R_c = 0 (no particle core) - _disc.parInvMM[_disc.offsetMetric[parType] + cell] += parts::dgtoolbox::invMMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 0.0).inverse() * static_cast(r_L); + _disc.parInvMM[_disc.offsetMetric[parType] + cell] += parts::dgtoolbox::mMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 0.0) * static_cast(r_L); _disc.parInvMM[_disc.offsetMetric[parType] + cell] = _disc.parInvMM[_disc.offsetMetric[parType] + cell].inverse(); _disc.minus_InvMM_ST[_disc.offsetMetric[parType] + cell] = -_disc.parInvMM[_disc.offsetMetric[parType] + cell] * _disc.parPolyDerM[parType].transpose() * _disc.parInvMM[_disc.offsetMetric[parType] + cell].inverse(); + + // particle GSM specific second order stiffness matrix (single element, i.e. nParCell = 1) + _disc.secondOrderStiffnessM[parType] = static_cast(_parCoreRadius[parType]) * parts::dgtoolbox::secondOrderStiffnessMatrix(_disc.parPolyDeg[parType], 0.0, 0.0, _disc.parNodes[parType]); + _disc.secondOrderStiffnessM[parType] += static_cast(_disc.deltaR[_disc.offsetMetric[parType]]) / 2.0 * parts::dgtoolbox::secondOrderStiffnessMatrix(_disc.parPolyDeg[parType], 0.0, 1.0, _disc.parNodes[parType]); } - else if (_parGeomSurfToVol[parType] == SurfVolRatioSlab) { // w = 1 + else if (_parGeomSurfToVol[parType] == SurfVolRatioSlab) // r = 1 + { + _disc.minus_InvMM_ST[_disc.offsetMetric[parType] + cell] = -_disc.parInvMM[_disc.offsetMetric[parType] + cell] * _disc.parPolyDerM[parType].transpose() * _disc.parInvMM[_disc.offsetMetric[parType] + cell].inverse(); + _disc.secondOrderStiffnessM[parType] = parts::dgtoolbox::secondOrderStiffnessMatrix(_disc.parPolyDeg[parType], 0.0, 0.0, _disc.parNodes[parType]); _disc.parInvMM[_disc.offsetMetric[parType] + cell] = parts::dgtoolbox::invMMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 0.0); } } + + _disc.minus_parInvMM_Ar[parType] = -_disc.parInvMM[_disc.offsetMetric[parType]] * _disc.secondOrderStiffnessM[parType]; } } diff --git a/src/libcadet/model/GeneralRateModelDG.hpp b/src/libcadet/model/GeneralRateModelDG.hpp index fde7b26af..14a36fa0d 100644 --- a/src/libcadet/model/GeneralRateModelDG.hpp +++ b/src/libcadet/model/GeneralRateModelDG.hpp @@ -270,7 +270,8 @@ class GeneralRateModelDG : public UnitOperationBase unsigned int* parPolyDeg; //!< polynomial degree of particle elements unsigned int* nParNode; //!< Array with number of radial nodes per cell in each particle type unsigned int* nParPoints; //!< Array with number of radial nodes per cell in each particle type - bool* parExactInt; //!< 1 for exact integration, 0 for inexact LGL quadrature for each particle type + bool* parExactInt; //!< 1 for exact integration, 0 for inexact LGL quadrature for each particle type + bool* parGSM; //!< specifies whether (single element) Galerkin spectral method should be used in particles unsigned int* parTypeOffset; //!< Array with offsets (in particle block) to particle type, additional last element contains total number of particle DOFs unsigned int* nBound; //!< Array with number of bound states for each component and particle type (particle type major ordering) unsigned int* boundOffset; //!< Array with offset to the first bound state of each component in the solid phase (particle type major ordering) @@ -296,6 +297,8 @@ class GeneralRateModelDG : public UnitOperationBase Eigen::VectorXd* parInvWeights; //!< Array with weights for LGL quadrature of size nNodes for each particle Eigen::MatrixXd* parInvMM; //!< dense inverse mass matrix for exact integration of integrals with metrics, for each particle Eigen::MatrixXd* parInvMM_Leg; //!< dense inverse mass matrix (Legendre) for exact integration of integral without metric, for each particle + Eigen::MatrixXd* secondOrderStiffnessM; //!< specific second order stiffness matrix + Eigen::MatrixXd* minus_parInvMM_Ar; //!< inverse mass matrix times specific second order stiffness matrix Eigen::Vector* Ir; //!< metric part for each particle type and cell, particle type major ordering Eigen::MatrixXd* Dr; //!< derivative matrices including metrics for each particle type and cell, particle type major ordering Eigen::VectorXi offsetMetric; //!< offset required to access metric dependent DG operator storage of Ir, Dr -> summed up nCells of all previous parTypes @@ -356,7 +359,9 @@ class GeneralRateModelDG : public UnitOperationBase Ir = new Vector[offsetMetric[nParType]]; minus_InvMM_ST = new MatrixXd[offsetMetric[nParType]]; parInvMM = new MatrixXd[offsetMetric[nParType]]; - + secondOrderStiffnessM = new MatrixXd[nParType]; + minus_parInvMM_Ar = new MatrixXd[nParType]; + /* compute metric independent DG operators for bulk and particles. Note that metric dependent DG operators are computet in updateRadialDisc(). */ for (int parType = 0; parType < nParType; parType++) @@ -372,9 +377,14 @@ class GeneralRateModelDG : public UnitOperationBase // particle jacobian blocks (each is unique) DGjacParDispBlocks = new MatrixXd[std::accumulate(nParCell, nParCell + nParType, 0)]; - for (unsigned int type = 0; type < nParType; type++) { - for (unsigned int block = 0; block < nParCell[type]; block++) { - DGjacParDispBlocks[offsetMetric[type] + block] = DGjacobianParDispBlock(block + 1u, type, parGeomSurfToVol[type]); + for (unsigned int type = 0; type < nParType; type++) + { + for (unsigned int block = 0; block < nParCell[type]; block++) + { + if (parGSM[type]) + DGjacParDispBlocks[offsetMetric[type] + block] = GSMjacobianParDispBlock(type, parGeomSurfToVol[type]); + else + DGjacParDispBlocks[offsetMetric[type] + block] = DGjacobianParDispBlock(block + 1u, type, parGeomSurfToVol[type]); } } } @@ -414,7 +424,7 @@ class GeneralRateModelDG : public UnitOperationBase gBlock *= 2.0 / static_cast(deltaR[offsetMetric[parType] + (cellIdx - 1)]); } else { - + // inexact integration not maintained due to inferior performance. Code is part of calcParticleCollocationDGSEMJacobian() } return gBlock; @@ -463,32 +473,50 @@ class GeneralRateModelDG : public UnitOperationBase } /** * @brief calculates the dispersion part of the DG jacobian - * @param [in] exInt true if exact integration DG scheme + * @param [in] parType particle type index + * @param [in] parGeomSurfToVol particle geometry + */ + Eigen::MatrixXd GSMjacobianParDispBlock(unsigned int parType, double parGeomSurfToVol) { + + MatrixXd dispBlock; + + // We have to match the DGSEM interface, where the dispersion block [ d RHS_disp / d c ] depends on whole previous and subsequent cell plus first entries of subsubsequent cells + dispBlock = MatrixXd::Zero(nParNode[parType], 3 * nParNode[parType] + 2); + + dispBlock.block(0, nParNode[parType] + 1, nParNode[parType], nParNode[parType]) = minus_parInvMM_Ar[parType]; + dispBlock *= 2.0 / static_cast(deltaR[offsetMetric[parType]]) * 2.0 / static_cast(deltaR[offsetMetric[parType]]); + + return -dispBlock; // *-1 for residual + } + /** + * @brief calculates the dispersion part of the DG jacobian * @param [in] cellIdx cell index + * @param [in] parType particle type index + * @param [in] parGeomSurfToVol particle geometry */ Eigen::MatrixXd DGjacobianParDispBlock(unsigned int cellIdx, unsigned int parType, double parGeomSurfToVol) { MatrixXd dispBlock; + // Inner dispersion block [ d RHS_disp / d c ], depends on whole previous and subsequent cell plus first entries of subsubsequent cells + dispBlock = MatrixXd::Zero(nParNode[parType], 3 * nParNode[parType] + 2); - if (parExactInt[parType]) { - // Inner dispersion block [ d RHS_disp / d c ], depends on whole previous and subsequent cell plus first entries of subsubsequent cells - dispBlock = MatrixXd::Zero(nParNode[parType], 3 * nParNode[parType] + 2); + if (parExactInt[parType]) + { MatrixXd B = getParBMatrix(parType, cellIdx, parGeomSurfToVol); // "Lifting" matrix MatrixXd gBlock = getParGBlock(cellIdx, parType); // current cell auxiliary block matrix MatrixXd gStarDC = parAuxBlockGstar(cellIdx, parType, getParGBlock(cellIdx - 1, parType), gBlock, getParGBlock(cellIdx + 1, parType)); // Numerical flux block - if (parGeomSurfToVol != SurfVolRatioSlab) { + if (parGeomSurfToVol != SurfVolRatioSlab) // weak form DGSEM required dispBlock.block(0, nParNode[parType], nParNode[parType], nParNode[parType] + 2) = minus_InvMM_ST[offsetMetric[parType] + (cellIdx - 1)] * gBlock; - dispBlock += parInvMM[offsetMetric[parType] + (cellIdx - 1)] * B * gStarDC; - } - else { - dispBlock.block(0, nParNode[parType], nParNode[parType], nParNode[parType] + 2) = (parPolyDerM[parType] - parInvMM[parType] * B) * gBlock; - dispBlock += parInvMM[parType] * B * gStarDC; - } + else // strong form DGSEM + dispBlock.block(0, nParNode[parType], nParNode[parType], nParNode[parType] + 2) = (parPolyDerM[parType] - parInvMM[offsetMetric[parType] + (cellIdx - 1)] * B) * gBlock; + + dispBlock += parInvMM[offsetMetric[parType] + (cellIdx - 1)] * B * gStarDC; dispBlock *= 2.0 / static_cast(deltaR[offsetMetric[parType] + (cellIdx - 1)]); } - else { - // inexact integration collocation DGSEM deprecated here + else + { + // inexact integration is not maintained due to inferior performance. Code is in calcParticleCollocationDGSEMJacobian } return -dispBlock; // *-1 for residual @@ -762,6 +790,15 @@ class GeneralRateModelDG : public UnitOperationBase } } + template + void parGSMVolumeIntegral(const int parType, Eigen::Map, 0, InnerStride>& state, Eigen::Map, 0, InnerStride>& stateDer) { + + int nNodes = _disc.nParNode[parType]; + + stateDer.segment(0, nNodes) + -= (_disc.minus_parInvMM_Ar[parType].template cast() * state.segment(0, nNodes)).template cast(); + } + template void parVolumeIntegral(const int parType, const bool aux, Eigen::Map, 0, InnerStride>& state, Eigen::Map, 0, InnerStride>& stateDer) { @@ -792,7 +829,7 @@ class GeneralRateModelDG : public UnitOperationBase if (_parGeomSurfToVol[parType] != _disc.SurfVolRatioSlab && _parCoreRadius[parType] == 0.0) { Cell0 = 1; - // estimate volume integral except for boundary node + // compute volume integral except for boundary node stateDer.segment(1, nNodes - 1) -= (_disc.Dr[_disc.offsetMetric[parType]].block(1, 1, nNodes - 1, nNodes - 1).template cast() * state.segment(1, nNodes - 1)).template cast(); // estimate volume integral for boundary node: sum_{j=1}^N state_j * w_j * D_{j,0} * r_j stateDer[0] += static_cast( @@ -849,7 +886,7 @@ class GeneralRateModelDG : public UnitOperationBase // film diffusion BC _surfFluxPar[_disc.nParCell[parType]] = static_cast(_disc.localFlux[comp]) / (static_cast(_parPorosity[parType]) * static_cast(_poreAccessFactor[parType * _disc.nComp + comp])) - * (2.0 / static_cast(_disc.deltaR[_disc.offsetMetric[parType]])); // inverse squared mapping is also applied, so we apply Map * invMap^2 = invMap + * (2.0 / static_cast(_disc.deltaR[_disc.offsetMetric[parType]])); // inverse squared mapping was also applied, so we apply Map * invMap^2 = invMap // inner particle BC _surfFluxPar[0] = 0.0; diff --git a/src/libcadet/model/parts/DGToolbox.cpp b/src/libcadet/model/parts/DGToolbox.cpp index 3955308ac..aa1850e4a 100644 --- a/src/libcadet/model/parts/DGToolbox.cpp +++ b/src/libcadet/model/parts/DGToolbox.cpp @@ -305,7 +305,7 @@ VectorXd barycentricWeights(const unsigned int polyDeg, const VectorXd nodes) { return baryWeights; } /** - * @brief computation of nodal (lagrange) polynomial derivative matrix + * @brief computation of nodal polynomial derivative matrix * @param [in] polyDeg polynomial degree * @param [in] nodes polynomial interpolation nodes */ @@ -341,7 +341,7 @@ double orthonFactor(const int polyDeg, double a, double b) { * @param [in] a Jacobi polynomial parameter * @param [in] b Jacobi polynomial parameter */ -MatrixXd getVandermonde_JACOBI(const unsigned int polyDeg, const VectorXd nodes, const double a, const double b) { +MatrixXd jacVandermondeMatrix(const unsigned int polyDeg, const VectorXd nodes, const double a, const double b) { const unsigned int nNodes = polyDeg + 1u; MatrixXd V(nNodes, nNodes); @@ -369,22 +369,65 @@ MatrixXd getVandermonde_JACOBI(const unsigned int polyDeg, const VectorXd nodes, return V; } +double jacPDerivativePreFactor(const unsigned int pIndex, const unsigned int derOrder, const double a, const double b) { + + double prefac = std::tgamma(a + b + static_cast(pIndex) + 1.0 + static_cast(derOrder)) / (std::pow(2.0, static_cast(derOrder)) * std::tgamma(a + b + static_cast(pIndex) + 1.0)); + + return prefac * orthonFactor(pIndex - derOrder, a + static_cast(derOrder), b + static_cast(derOrder)) / orthonFactor(pIndex, a, b); +} /** * @brief calculates the Vandermonde matrix of the normalized Legendre polynomials * @param [in] polyDeg polynomial degree * @param [in] nodes polynomial interpolation nodes */ -MatrixXd getVandermonde_LEGENDRE(const unsigned int polyDeg, const VectorXd nodes) { - return getVandermonde_JACOBI(polyDeg, nodes, 0.0, 0.0); +MatrixXd legVandermondeMatrix(const unsigned int polyDeg, const VectorXd nodes) { + return jacVandermondeMatrix(polyDeg, nodes, 0.0, 0.0); } /** - * @brief calculates mass matrix via transformation to orthonormal Jacobi (modal) basis - * @detail exact integration for integrals of the form \int_E \ell_i(\xi) \ell_j(\xi) (1 - \xi)^\alpha (1 + \xi)^\beta d\xi + * @brief calculates the inverse mass matrix via transformation to orthonormal Jacobi (modal) basis + * @detail the mass matrix used to compute integrals of the form \int_E \ell_i(\xi) \ell_j(\xi) (1 - \xi)^\alpha (1 + \xi)^\beta d\xi * @param [in] polyDeg polynomial degree * @param [in] nodes polynomial interpolation nodes */ Eigen::MatrixXd invMMatrix(const unsigned int polyDeg, const Eigen::VectorXd nodes, const double alpha, const double beta) { - return (getVandermonde_JACOBI(polyDeg, nodes, alpha, beta) * (getVandermonde_JACOBI(polyDeg, nodes, alpha, beta).transpose())); + return (jacVandermondeMatrix(polyDeg, nodes, alpha, beta) * (jacVandermondeMatrix(polyDeg, nodes, alpha, beta).transpose())); +} +/** + * @brief calculates the mass matrix via transformation to orthonormal Jacobi (modal) basis + * @detail mass matrix used to compute integrals of the form \int_E \ell_i(\xi) \ell_j(\xi) (1 - \xi)^\alpha (1 + \xi)^\beta d\xi + * @param [in] polyDeg polynomial degree + * @param [in] nodes polynomial interpolation nodes + */ +Eigen::MatrixXd mMatrix(const unsigned int polyDeg, const Eigen::VectorXd nodes, const double alpha, const double beta) { + return invMMatrix(polyDeg, nodes, alpha, beta).inverse(); +} +/** + * @brief calculates the derivative of the Vandermonde matrix of the normalized Legendre polynomials + * @param [in] polyDeg polynomial degree + * @param [in] a Jacobi polynomial parameter + * @param [in] b Jacobi polynomial parameter + * @param [in] nodes polynomial interpolation nodes + */ +MatrixXd jacDerVandermondeMatrix(const unsigned int polyDeg, const double a, const double b, const VectorXd nodes) { + + MatrixXd derVan = MatrixXd::Zero(polyDeg + 1, polyDeg + 1); + derVan.block(0, 1, polyDeg + 1, polyDeg) = jacVandermondeMatrix(polyDeg, nodes, a + 1.0, b + 1.0).block(0, 0, polyDeg + 1, polyDeg); + + for (int hm = 1; hm < polyDeg + 1; hm++) + derVan.block(0, hm, polyDeg + 1, 1) *= jacPDerivativePreFactor(hm, 1, a, b); + //derVan.block(0, hm, polyDeg + 1, 1) *= std::sqrt(hm * (hm + a + b + 1)); // todo + + return derVan; +} +/** + * @brief calculates the second order nodal stiffness matrix (via transformation to normalized Legendre polynomials) + * @param [in] polyDeg polynomial degree + * @param [in] a Jacobi polynomial parameter + * @param [in] b Jacobi polynomial parameter + * @param [in] nodes polynomial interpolation nodes + */ +MatrixXd secondOrderStiffnessMatrix(const unsigned int polyDeg, const double alpha, const double beta, const VectorXd nodes) { + return derivativeMatrix(polyDeg, nodes).transpose() * mMatrix(polyDeg, nodes, alpha, beta) * derivativeMatrix(polyDeg, nodes); } } // namespace dgtoolbox diff --git a/src/libcadet/model/parts/DGToolbox.hpp b/src/libcadet/model/parts/DGToolbox.hpp index c6e482547..6f4fdfc08 100644 --- a/src/libcadet/model/parts/DGToolbox.hpp +++ b/src/libcadet/model/parts/DGToolbox.hpp @@ -69,12 +69,28 @@ Eigen::VectorXd barycentricWeights(const unsigned int polyDeg, const Eigen::Vect */ Eigen::MatrixXd derivativeMatrix(const unsigned int polyDeg, const Eigen::VectorXd nodes); /** - * @brief calculates mass matrix via transformation to orthonormal Jacobi (modal) basis - * @detail exact integration for integrals of the form \int_E \ell_i(\xi) \ell_j(\xi) (1 - \xi)^\alpha (1 + \xi)^\beta d\xi + * @brief calculates the inverse mass matrix via transformation to orthonormal Jacobi (modal) basis + * @detail the mass matrix used to compute integrals of the form \int_E \ell_i(\xi) \ell_j(\xi) (1 - \xi)^\alpha (1 + \xi)^\beta d\xi * @param [in] polyDeg polynomial degree * @param [in] nodes polynomial interpolation nodes */ Eigen::MatrixXd invMMatrix(const unsigned int polyDeg, const Eigen::VectorXd nodes, const double alpha = 0.0, const double beta = 0.0); +/** + * @brief calculates the mass matrix via transformation to orthonormal Jacobi (modal) basis + * @detail mass matrix used to compute integrals of the form \int_E \ell_i(\xi) \ell_j(\xi) (1 - \xi)^\alpha (1 + \xi)^\beta d\xi + * @param [in] polyDeg polynomial degree + * @param [in] nodes polynomial interpolation nodes + */ +Eigen::MatrixXd mMatrix(const unsigned int polyDeg, const Eigen::VectorXd nodes, const double alpha, const double beta); +/** + * @brief calculates a specific second order nodal stiffness matrix + * @detail for integrals including terms of the form (1 - \xi)^\alpha (1 + \xi)^\beta. Computation via transformation to the respective Jacobi polynomial + * @param [in] polyDeg polynomial degree + * @param [in] a Jacobi polynomial parameter + * @param [in] b Jacobi polynomial parameter + * @param [in] nodes polynomial interpolation nodes + */ +Eigen::MatrixXd secondOrderStiffnessMatrix(const unsigned int polyDeg, const double alpha, const double beta, const Eigen::VectorXd nodes); } // namespace dgtoolbox } // namespace parts