diff --git a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp index 99e04cb28..210e78f76 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp @@ -141,6 +141,14 @@ bool AxialConvectionDispersionOperatorBaseDG::configureModelDiscretization(IPara _OSmode = 2; int order = paramProvider.exists("FV_ORDER") ? paramProvider.getInt("FV_ORDER") : 1; + if (order <= 0) + throw InvalidParameterException("Low order blending scheme must be of order >= 1, but was specified as " + std::to_string(order)); + else if(order > 1) + throw InvalidParameterException("Low order blending scheme of order > 1, not implemented yet. Please use order 1."); + + _maxBlending = paramProvider.exists("MAX_BLENDING") ? paramProvider.getDouble("MAX_BLENDING") : 1.0; + if(_maxBlending < 0.0) + throw InvalidParameterException("Low order blending factor must be >= 0.0, but was specified as " + std::to_string(_maxBlending)); _modalVanInv.resize(_nNodes, _nNodes); _modalVanInv.setZero(); @@ -503,56 +511,63 @@ 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 FVblending = 0.0; + double FVblending = _maxBlending; // todo use indicator if (_OSmode == 2) { - FVblending = 1.0; _boundary[0] = y[comp]; // copy inlet DOFs to ghost node - subcellFVintegral(FVblending, y + offsetC() + comp, res + offsetC() + comp, _strideNode, _strideNode); + if (FVblending > 0.0 && FVblending < 1.0) // only use lower order scheme for advection, treat dispersion with high order DG. + subcellFVconvectionIntegral(FVblending, y + offsetC() + comp, res + offsetC() + comp, _strideNode, _strideNode); + else if (FVblending == 1.0) { // use lower order scheme for both, advection and dispersion. // TODO: either special mode where this is done or deprecated when troubled cell identification works + subcellFVintegral(FVblending, y + offsetC() + comp, res + offsetC() + comp, d_ax, _strideNode, _strideNode); + continue; + } } - else - { - // ===================================// - // reset cache // - // ===================================// - _h.setZero(); - _g.setZero(); - _boundary[0] = y[comp]; // copy inlet DOFs to ghost node + double DGblending = 1.0 - FVblending; - // ======================================// - // solve auxiliary system g = d c / d x // - // ======================================// + // ===================================// + // reset cache // + // ===================================// - // DG volume integral in strong form - volumeIntegral(_C, _g); + _h.setZero(); + _g.setZero(); + _boundary[0] = y[comp]; // copy inlet DOFs to ghost node - // calculate numerical flux values c* - InterfaceFluxAuxiliary(y + offsetC() + comp, _strideNode, _strideCell); + // ======================================// + // solve auxiliary system g = d c / d x // + // ======================================// - // DG surface integral in strong form - surfaceIntegral(y + offsetC() + comp, reinterpret_cast(&_g[0]), - _strideNode, _strideCell, 1u, _nNodes); + // DG volume integral in strong form + volumeIntegral(_C, _g); - // ======================================// - // solve main equation RHS d h / d x // - // ======================================// + // calculate numerical flux values c* + InterfaceFluxAuxiliary(y + offsetC() + comp, _strideNode, _strideCell); - // 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(); + // DG surface integral in strong form + surfaceIntegral(y + offsetC() + comp, reinterpret_cast(&_g[0]), + _strideNode, _strideCell, 1u, _nNodes); - // DG volume integral in strong form - volumeIntegral(_h, _resC); + // ======================================// + // solve main equation RHS d h / d x // + // ======================================// - // update boundary values for auxiliary variable g (solid wall) - calcBoundaryValues(); + // calculate the substitute h(g(c), c) and apply inverse mapping jacobian (reference space) + _h = 2.0 / static_cast(_deltaZ) * (-u * _C * DGblending + d_ax * (-2.0 / static_cast(_deltaZ)) * _g).template cast(); - // calculate numerical flux values h* - InterfaceFlux(y + offsetC() + comp, d_ax); + // DG volume integral in strong form + volumeIntegral(_h, _resC); - // DG surface integral in strong form - surfaceIntegral(&_h[0], res + offsetC() + comp, - 1u, _nNodes, _strideNode, _strideCell); - } + // update boundary values for auxiliary variable g (solid wall) + calcBoundaryValues(); + + // calculate numerical flux values h* + InterfaceFlux(y + offsetC() + comp, d_ax); + + // Not needed, when re-accounted for in FV subcell residual (i.e. partly reverted FV subcell DG weak form formulation) + //_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); } return 0; @@ -598,9 +613,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 8245342fd..1f1c7984a 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp @@ -189,6 +189,7 @@ namespace cadet Eigen::VectorXd _modalCoeff; //!< Orthonormal modal polynomial coefficients (used by smoothness indicator) Eigen::MatrixXd _modalVanInv; //!< Inverse Vandermonde matrix of orthonormal modal polynomial coefficients (used by smoothness indicator) Eigen::VectorXd _weights; //!< LGL quadrature weights + double _maxBlending; Eigen::Vector _g; //!< auxiliary variable Eigen::Vector _h; //!< auxiliary substitute @@ -769,7 +770,7 @@ namespace cadet } template - void subcellFVintegral(double blending, const StateType* _C, ResidualType* stateDer, + void subcellFVintegral(double blending, const StateType* _C, ResidualType* stateDer, ParamType d_ax, unsigned int strideNodeC, unsigned int strideNode_stateDer) { const StateType* localC = _C; @@ -789,7 +790,8 @@ namespace cadet // subcell faces for (int interface = 1; interface < _nNodes; interface++, localC += strideNodeC, localStateDer += strideNode_stateDer) { - flux = blending * static_cast(_curVelocity * localC[0]); // left value (upwind) + flux = blending * (static_cast(_curVelocity * localC[0]) // upwind flux for convection + + 0.5 * d_ax * (localC[0] + localC[1]) * (2 / (_invWeights[interface - 1] + _invWeights[interface]))); // central flux for diffusion with first order FD approximation of c_x localStateDer[0] += 2.0 / static_cast(_deltaZ) * _invWeights[interface - 1] * flux; localStateDer[strideNode_stateDer] -= 2.0 / static_cast(_deltaZ) * _invWeights[interface] * flux; @@ -807,15 +809,62 @@ namespace cadet for (unsigned int elemFace = 1; elemFace < _nCells; elemFace++, localC += _nNodes * strideNodeC, localStateDer += _nNodes * strideNode_stateDer) { - flux = blending * static_cast(_curVelocity * localC[0]); + flux = blending * (static_cast(_curVelocity * localC[0]) // upwind flux for convection + + 0.5 * d_ax * (localC[0] + localC[1]) * (2 / (_invWeights[_polyDeg] + _invWeights[0]))); // central flux for diffusion with first order FD approximation of c_x + localStateDer[0] += 2.0 / static_cast(_deltaZ) * _invWeights[_polyDeg] * flux; localStateDer[strideNode_stateDer] -= 2.0 / static_cast(_deltaZ) * _invWeights[0] * flux; } // outlet boundary BC // todo backwards flow! (will be solved once DG operation is being used) - flux = blending * static_cast(_curVelocity * _C[_nCells * (strideNodeC * _nNodes) - strideNodeC]); + flux = blending * (static_cast(_curVelocity * _C[_nCells * (strideNodeC * _nNodes) - strideNodeC]) // upwind flux for convection + + 0.5 * d_ax * (localC[0] + localC[1]) * (2 / (_invWeights[_polyDeg] + _invWeights[0]))); // central flux for diffusion with first order FD approximation of c_x + stateDer[_nCells * (strideNode_stateDer * _nNodes) - strideNode_stateDer] += 2.0 / static_cast(_deltaZ) * _invWeights[_polyDeg] * flux; } + template + void subcellFVconvectionIntegral(double blending, const StateType* _C, ResidualType* stateDer, + unsigned int strideNodeC, unsigned int strideNode_stateDer) { + + const StateType* localC = _C; + ResidualType* localStateDer = stateDer; + ResidualType flux = 0.0; + + // TODO backwards flow! + + /* subcell FV faces */ + + for (unsigned int cell = 0; cell < _nCells; cell++, localC += strideNodeC, localStateDer += strideNode_stateDer) { + + //FVintegral(blending, cell, &_C[cell * (strideNodeC * _nNodes)], &stateDer[cell * (strideNode_stateDer * _nNodes)], strideNodeC, strideNode_stateDer); + + //if (!indicator) { // todo use indicator + // localC += + // localStateDer += + // continue; + //} + + for (int interface = 1; interface < _nNodes; interface++, localC += strideNodeC, localStateDer += strideNode_stateDer) { + + flux = blending * static_cast(_curVelocity * localC[0]); // upwind flux for convection + + localStateDer[0] += 2.0 / static_cast(_deltaZ) * _invWeights[interface - 1] * flux; + localStateDer[strideNode_stateDer] -= 2.0 / static_cast(_deltaZ) * _invWeights[interface] * flux; + } + } + + /* DG element interfaces : account DG strong - form like formulation of FV subcell scheme */ + // Not needed, when re-accounted for in DG surface integral (i.e. partly reverted FV subcell DG weak form formulation) + //localC = _C; + //localStateDer = stateDer; + + //for (unsigned int DGcell = 0; DGcell < _nCells; DGcell++, localC += _nNodes * strideNodeC, localStateDer += _nNodes * strideNode_stateDer) { + // // todo use indicator + + // localStateDer[0] -= 2.0 / static_cast(_deltaZ) * _invWeights[0] * blending * static_cast(_curVelocity) * localC[0]; + // localStateDer[_polyDeg * strideNode_stateDer] += 2.0 / static_cast(_deltaZ) * _invWeights[_polyDeg] * blending * static_cast(_curVelocity) * localC[_polyDeg * strideNodeC]; + //} + } /** * @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 @@ -835,7 +884,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 */ @@ -963,7 +1012,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 */