Skip to content

Commit a812d54

Browse files
jbreue16jbreue16
committed
Add 1st order FV for subcell blending
1 parent 2fa89de commit a812d54

File tree

2 files changed

+108
-44
lines changed

2 files changed

+108
-44
lines changed

‎src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp

Lines changed: 53 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -141,6 +141,14 @@ bool AxialConvectionDispersionOperatorBaseDG::configureModelDiscretization(IPara
141141

142142
_OSmode = 2;
143143
int order = paramProvider.exists("FV_ORDER") ? paramProvider.getInt("FV_ORDER") : 1;
144+
if (order <= 0)
145+
throw InvalidParameterException("Low order blending scheme must be of order >= 1, but was specified as " + std::to_string(order));
146+
else if(order > 1)
147+
throw InvalidParameterException("Low order blending scheme of order > 1, not implemented yet. Please use order 1.");
148+
149+
_maxBlending = paramProvider.exists("MAX_BLENDING") ? paramProvider.getDouble("MAX_BLENDING") : 1.0;
150+
if(_maxBlending < 0.0)
151+
throw InvalidParameterException("Low order blending factor must be >= 0.0, but was specified as " + std::to_string(_maxBlending));
144152

145153
_modalVanInv.resize(_nNodes, _nNodes);
146154
_modalVanInv.setZero();
@@ -503,56 +511,63 @@ int AxialConvectionDispersionOperatorBaseDG::residualImpl(const IModel& model, d
503511
const ParamType u = static_cast<ParamType>(_curVelocity);
504512
const ParamType d_ax = static_cast<ParamType>(getSectionDependentSlice(_colDispersion, _nComp, secIdx)[comp]);
505513

506-
double FVblending = 0.0;
514+
double FVblending = _maxBlending; // todo use indicator
507515
if (_OSmode == 2) {
508-
FVblending = 1.0;
509516
_boundary[0] = y[comp]; // copy inlet DOFs to ghost node
510-
subcellFVintegral<StateType, ResidualType, ParamType>(FVblending, y + offsetC() + comp, res + offsetC() + comp, _strideNode, _strideNode);
517+
if (FVblending > 0.0 && FVblending < 1.0) // only use lower order scheme for advection, treat dispersion with high order DG.
518+
subcellFVconvectionIntegral<StateType, ResidualType, ParamType>(FVblending, y + offsetC() + comp, res + offsetC() + comp, _strideNode, _strideNode);
519+
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
520+
subcellFVintegral<StateType, ResidualType, ParamType>(FVblending, y + offsetC() + comp, res + offsetC() + comp, d_ax, _strideNode, _strideNode);
521+
continue;
522+
}
511523
}
512-
else
513-
{
514-
// ===================================//
515-
// reset cache //
516-
// ===================================//
517524

518-
_h.setZero();
519-
_g.setZero();
520-
_boundary[0] = y[comp]; // copy inlet DOFs to ghost node
525+
double DGblending = 1.0 - FVblending;
521526

522-
// ======================================//
523-
// solve auxiliary system g = d c / d x //
524-
// ======================================//
527+
// ===================================//
528+
// reset cache //
529+
// ===================================//
525530

526-
// DG volume integral in strong form
527-
volumeIntegral<StateType, StateType>(_C, _g);
531+
_h.setZero();
532+
_g.setZero();
533+
_boundary[0] = y[comp]; // copy inlet DOFs to ghost node
528534

529-
// calculate numerical flux values c*
530-
InterfaceFluxAuxiliary<StateType>(y + offsetC() + comp, _strideNode, _strideCell);
535+
// ======================================//
536+
// solve auxiliary system g = d c / d x //
537+
// ======================================//
531538

532-
// DG surface integral in strong form
533-
surfaceIntegral<StateType, StateType>(y + offsetC() + comp, reinterpret_cast<StateType*>(&_g[0]),
534-
_strideNode, _strideCell, 1u, _nNodes);
539+
// DG volume integral in strong form
540+
volumeIntegral<StateType, StateType>(_C, _g);
535541

536-
// ======================================//
537-
// solve main equation RHS d h / d x //
538-
// ======================================//
542+
// calculate numerical flux values c*
543+
InterfaceFluxAuxiliary<StateType>(y + offsetC() + comp, _strideNode, _strideCell);
539544

540-
// calculate the substitute h(g(c), c) and apply inverse mapping jacobian (reference space)
541-
_h = 2.0 / static_cast<ParamType>(_deltaZ) * (-u * _C + d_ax * (-2.0 / static_cast<ParamType>(_deltaZ)) * _g).template cast<ResidualType>();
545+
// DG surface integral in strong form
546+
surfaceIntegral<StateType, StateType>(y + offsetC() + comp, reinterpret_cast<StateType*>(&_g[0]),
547+
_strideNode, _strideCell, 1u, _nNodes);
542548

543-
// DG volume integral in strong form
544-
volumeIntegral<ResidualType, ResidualType>(_h, _resC);
549+
// ======================================//
550+
// solve main equation RHS d h / d x //
551+
// ======================================//
545552

546-
// update boundary values for auxiliary variable g (solid wall)
547-
calcBoundaryValues<StateType>();
553+
// calculate the substitute h(g(c), c) and apply inverse mapping jacobian (reference space)
554+
_h = 2.0 / static_cast<ParamType>(_deltaZ) * (-u * _C * DGblending + d_ax * (-2.0 / static_cast<ParamType>(_deltaZ)) * _g).template cast<ResidualType>();
548555

549-
// calculate numerical flux values h*
550-
InterfaceFlux<StateType, ParamType>(y + offsetC() + comp, d_ax);
556+
// DG volume integral in strong form
557+
volumeIntegral<ResidualType, ResidualType>(_h, _resC);
551558

552-
// DG surface integral in strong form
553-
surfaceIntegral<ResidualType, ResidualType>(&_h[0], res + offsetC() + comp,
554-
1u, _nNodes, _strideNode, _strideCell);
555-
}
559+
// update boundary values for auxiliary variable g (solid wall)
560+
calcBoundaryValues<StateType>();
561+
562+
// calculate numerical flux values h*
563+
InterfaceFlux<StateType, ParamType>(y + offsetC() + comp, d_ax);
564+
565+
// Not needed, when re-accounted for in FV subcell residual (i.e. partly reverted FV subcell DG weak form formulation)
566+
//_h += 2.0 / static_cast<ParamType>(_deltaZ) * (-u * _C * FVblending).template cast<ResidualType>();
567+
568+
// DG surface integral in strong form
569+
surfaceIntegral<ResidualType, ResidualType>(&_h[0], res + offsetC() + comp,
570+
1u, _nNodes, _strideNode, _strideCell);
556571
}
557572

558573
return 0;
@@ -598,9 +613,9 @@ unsigned int AxialConvectionDispersionOperatorBaseDG::nConvDispEntries(bool pure
598613
void model::parts::AxialConvectionDispersionOperatorBaseDG::convDispJacPattern(std::vector<T>& tripletList, const int bulkOffset)
599614
{
600615
if (_exactInt)
601-
ConvDispModalPattern(tripletList, bulkOffset);
616+
ConvDispExIntDGPattern(tripletList, bulkOffset);
602617
else
603-
ConvDispNodalPattern(tripletList, bulkOffset);
618+
ConvDispCollocationDGPattern(tripletList, bulkOffset);
604619
}
605620
/**
606621
* @brief Multiplies the time derivative Jacobian @f$ \frac{\partial F}{\partial \dot{y}}\left(t, y, \dot{y}\right) @f$ with a given vector

‎src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp

Lines changed: 55 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -189,6 +189,7 @@ namespace cadet
189189
Eigen::VectorXd _modalCoeff; //!< Orthonormal modal polynomial coefficients (used by smoothness indicator)
190190
Eigen::MatrixXd _modalVanInv; //!< Inverse Vandermonde matrix of orthonormal modal polynomial coefficients (used by smoothness indicator)
191191
Eigen::VectorXd _weights; //!< LGL quadrature weights
192+
double _maxBlending;
192193

193194
Eigen::Vector<active, Eigen::Dynamic> _g; //!< auxiliary variable
194195
Eigen::Vector<active, Eigen::Dynamic> _h; //!< auxiliary substitute
@@ -769,7 +770,7 @@ namespace cadet
769770
}
770771

771772
template<typename StateType, typename ResidualType, typename ParamType>
772-
void subcellFVintegral(double blending, const StateType* _C, ResidualType* stateDer,
773+
void subcellFVintegral(double blending, const StateType* _C, ResidualType* stateDer, ParamType d_ax,
773774
unsigned int strideNodeC, unsigned int strideNode_stateDer) {
774775

775776
const StateType* localC = _C;
@@ -789,7 +790,8 @@ namespace cadet
789790
// subcell faces
790791
for (int interface = 1; interface < _nNodes; interface++, localC += strideNodeC, localStateDer += strideNode_stateDer) {
791792

792-
flux = blending * static_cast<ResidualType>(_curVelocity * localC[0]); // left value (upwind)
793+
flux = blending * (static_cast<ResidualType>(_curVelocity * localC[0]) // upwind flux for convection
794+
+ 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
793795

794796
localStateDer[0] += 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[interface - 1] * flux;
795797
localStateDer[strideNode_stateDer] -= 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[interface] * flux;
@@ -807,15 +809,62 @@ namespace cadet
807809

808810
for (unsigned int elemFace = 1; elemFace < _nCells; elemFace++, localC += _nNodes * strideNodeC, localStateDer += _nNodes * strideNode_stateDer) {
809811

810-
flux = blending * static_cast<ResidualType>(_curVelocity * localC[0]);
812+
flux = blending * (static_cast<ResidualType>(_curVelocity * localC[0]) // upwind flux for convection
813+
+ 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
814+
811815
localStateDer[0] += 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[_polyDeg] * flux;
812816
localStateDer[strideNode_stateDer] -= 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[0] * flux;
813817
}
814818
// outlet boundary BC // todo backwards flow! (will be solved once DG operation is being used)
815-
flux = blending * static_cast<ResidualType>(_curVelocity * _C[_nCells * (strideNodeC * _nNodes) - strideNodeC]);
819+
flux = blending * (static_cast<ResidualType>(_curVelocity * _C[_nCells * (strideNodeC * _nNodes) - strideNodeC]) // upwind flux for convection
820+
+ 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
821+
816822
stateDer[_nCells * (strideNode_stateDer * _nNodes) - strideNode_stateDer] += 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[_polyDeg] * flux;
817823

818824
}
825+
template<typename StateType, typename ResidualType, typename ParamType>
826+
void subcellFVconvectionIntegral(double blending, const StateType* _C, ResidualType* stateDer,
827+
unsigned int strideNodeC, unsigned int strideNode_stateDer) {
828+
829+
const StateType* localC = _C;
830+
ResidualType* localStateDer = stateDer;
831+
ResidualType flux = 0.0;
832+
833+
// TODO backwards flow!
834+
835+
/* subcell FV faces */
836+
837+
for (unsigned int cell = 0; cell < _nCells; cell++, localC += strideNodeC, localStateDer += strideNode_stateDer) {
838+
839+
//FVintegral<StateType, ResidualType, ParamType>(blending, cell, &_C[cell * (strideNodeC * _nNodes)], &stateDer[cell * (strideNode_stateDer * _nNodes)], strideNodeC, strideNode_stateDer);
840+
841+
//if (!indicator) { // todo use indicator
842+
// localC +=
843+
// localStateDer +=
844+
// continue;
845+
//}
846+
847+
for (int interface = 1; interface < _nNodes; interface++, localC += strideNodeC, localStateDer += strideNode_stateDer) {
848+
849+
flux = blending * static_cast<ResidualType>(_curVelocity * localC[0]); // upwind flux for convection
850+
851+
localStateDer[0] += 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[interface - 1] * flux;
852+
localStateDer[strideNode_stateDer] -= 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[interface] * flux;
853+
}
854+
}
855+
856+
/* DG element interfaces : account DG strong - form like formulation of FV subcell scheme */
857+
// Not needed, when re-accounted for in DG surface integral (i.e. partly reverted FV subcell DG weak form formulation)
858+
//localC = _C;
859+
//localStateDer = stateDer;
860+
861+
//for (unsigned int DGcell = 0; DGcell < _nCells; DGcell++, localC += _nNodes * strideNodeC, localStateDer += _nNodes * strideNode_stateDer) {
862+
// // todo use indicator
863+
864+
// localStateDer[0] -= 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[0] * blending * static_cast<ResidualType>(_curVelocity) * localC[0];
865+
// localStateDer[_polyDeg * strideNode_stateDer] += 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[_polyDeg] * blending * static_cast<ResidualType>(_curVelocity) * localC[_polyDeg * strideNodeC];
866+
//}
867+
}
819868
/**
820869
* @brief computes ghost nodes to implement boundary conditions
821870
* @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
835884
/**
836885
* @brief sets the sparsity pattern of the convection dispersion Jacobian for the nodal DG scheme
837886
*/
838-
int ConvDispNodalPattern(std::vector<T>& tripletList, const int offC = 0) {
887+
int ConvDispCollocationDGPattern(std::vector<T>& tripletList, const int offC = 0) {
839888

840889
/*======================================================*/
841890
/* Define Convection Jacobian Block */
@@ -963,7 +1012,7 @@ namespace cadet
9631012
/**
9641013
* @brief sets the sparsity pattern of the convection dispersion Jacobian for the exact integration (her: modal) DG scheme
9651014
*/
966-
int ConvDispModalPattern(std::vector<T>& tripletList, const int offC = 0) {
1015+
int ConvDispExIntDGPattern(std::vector<T>& tripletList, const int offC = 0) {
9671016

9681017
/*======================================================*/
9691018
/* Define Convection Jacobian Block */

0 commit comments

Comments
 (0)