Skip to content

Commit

Permalink
Add 1st order FV for subcell blending
Browse files Browse the repository at this point in the history
  • Loading branch information
@jbreue16
jbreue16 committed Jun 18, 2023
1 parent 276500e commit 4473fa7
Show file tree
Hide file tree
Showing 2 changed files with 108 additions and 44 deletions.
91 changes: 53 additions & 38 deletions src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -503,56 +511,63 @@ int AxialConvectionDispersionOperatorBaseDG::residualImpl(const IModel& model, d
const ParamType u = static_cast<ParamType>(_curVelocity);
const ParamType d_ax = static_cast<ParamType>(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<StateType, ResidualType, ParamType>(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<StateType, ResidualType, ParamType>(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<StateType, ResidualType, ParamType>(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<StateType, StateType>(_C, _g);
_h.setZero();
_g.setZero();
_boundary[0] = y[comp]; // copy inlet DOFs to ghost node

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

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

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

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

// DG volume integral in strong form
volumeIntegral<ResidualType, ResidualType>(_h, _resC);
// ======================================//
// solve main equation RHS d h / d x //
// ======================================//

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

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

// DG surface integral in strong form
surfaceIntegral<ResidualType, ResidualType>(&_h[0], res + offsetC() + comp,
1u, _nNodes, _strideNode, _strideCell);
}
// update boundary values for auxiliary variable g (solid wall)
calcBoundaryValues<StateType>();

// calculate numerical flux values h*
InterfaceFlux<StateType, ParamType>(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<ParamType>(_deltaZ) * (-u * _C * FVblending).template cast<ResidualType>();

// DG surface integral in strong form
surfaceIntegral<ResidualType, ResidualType>(&_h[0], res + offsetC() + comp,
1u, _nNodes, _strideNode, _strideCell);
}

return 0;
Expand Down Expand Up @@ -598,9 +613,9 @@ unsigned int AxialConvectionDispersionOperatorBaseDG::nConvDispEntries(bool pure
void model::parts::AxialConvectionDispersionOperatorBaseDG::convDispJacPattern(std::vector<T>& 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
Expand Down
61 changes: 55 additions & 6 deletions src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<active, Eigen::Dynamic> _g; //!< auxiliary variable
Eigen::Vector<active, Eigen::Dynamic> _h; //!< auxiliary substitute
Expand Down Expand Up @@ -769,7 +770,7 @@ namespace cadet
}

template<typename StateType, typename ResidualType, typename ParamType>
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;
Expand All @@ -789,7 +790,8 @@ namespace cadet
// subcell faces
for (int interface = 1; interface < _nNodes; interface++, localC += strideNodeC, localStateDer += strideNode_stateDer) {

flux = blending * static_cast<ResidualType>(_curVelocity * localC[0]); // left value (upwind)
flux = blending * (static_cast<ResidualType>(_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<ParamType>(_deltaZ) * _invWeights[interface - 1] * flux;
localStateDer[strideNode_stateDer] -= 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[interface] * flux;
Expand All @@ -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<ResidualType>(_curVelocity * localC[0]);
flux = blending * (static_cast<ResidualType>(_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<ParamType>(_deltaZ) * _invWeights[_polyDeg] * flux;
localStateDer[strideNode_stateDer] -= 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[0] * flux;
}
// outlet boundary BC // todo backwards flow! (will be solved once DG operation is being used)
flux = blending * static_cast<ResidualType>(_curVelocity * _C[_nCells * (strideNodeC * _nNodes) - strideNodeC]);
flux = blending * (static_cast<ResidualType>(_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<ParamType>(_deltaZ) * _invWeights[_polyDeg] * flux;

}
template<typename StateType, typename ResidualType, typename ParamType>
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<StateType, ResidualType, ParamType>(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<ResidualType>(_curVelocity * localC[0]); // upwind flux for convection

localStateDer[0] += 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[interface - 1] * flux;
localStateDer[strideNode_stateDer] -= 2.0 / static_cast<ParamType>(_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<ParamType>(_deltaZ) * _invWeights[0] * blending * static_cast<ResidualType>(_curVelocity) * localC[0];
// localStateDer[_polyDeg * strideNode_stateDer] += 2.0 / static_cast<ParamType>(_deltaZ) * _invWeights[_polyDeg] * blending * static_cast<ResidualType>(_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
Expand All @@ -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<T>& tripletList, const int offC = 0) {
int ConvDispCollocationDGPattern(std::vector<T>& tripletList, const int offC = 0) {

/*======================================================*/
/* Define Convection Jacobian Block */
Expand Down Expand Up @@ -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<T>& tripletList, const int offC = 0) {
int ConvDispExIntDGPattern(std::vector<T>& tripletList, const int offC = 0) {

/*======================================================*/
/* Define Convection Jacobian Block */
Expand Down

0 comments on commit 4473fa7

Please sign in to comment.