39 #ifndef _SDSTABILIZATION_HPP_ 40 #define _SDSTABILIZATION_HPP_ 43 #include <lifev/core/filter/GetPot.hpp> 44 #include <lifev/core/util/LifeDebug.hpp> 58 template<
typename MeshType,
typename DofType>
83 const ReferenceFE& refFE,
84 const QuadratureRule& quadRule,
117 template<
typename MatrixType,
typename VectorType>
118 void applySUPG (
const Real dt, MatrixType& matrix,
const VectorType& state );
128 template<
typename MatrixType,
typename VectorType>
129 void applySD (
const Real dt, MatrixType& matrix,
const VectorType& state );
145 template <
typename VectorType,
typename SourceType >
146 void applyRHS (
const Real dt, VectorType& vector,
const VectorType& state,
147 const SourceType& source,
const Real& time);
150 void showMe (std::ostream& output = std::cout)
const;
181 template <
typename VectorType>
183 VectorElemental& beta,
Real& coeffBeta,
Real& coeffDiv)
const;
186 void bgradu_bgradv (
const Real& coef, VectorElemental& vel, ElemMat& elmat,
const CurrentFE& fe,
190 void lapu_bgradv (
const Real& coef, VectorElemental& vel, MatrixElemental& elmat,
const CurrentFE& fe,
194 void gradp_bgradv (
const Real& coef, VectorElemental& vel, MatrixElemental& elmat,
const CurrentFE& fe)
const;
197 void lapu_gradq (
const Real& coef, MatrixElemental& elmat,
const CurrentFE& fe)
const;
200 template <
typename SourceType>
201 void f_bgradv (
const Real& coef, SourceType& source, VectorElemental& vel,
202 VectorElemental& elvec,
const CurrentFE& fe,
UInt iblock,
const Real& time)
const;
205 template<
typename SourceType>
206 void f_gradq (
const Real& coef, SourceType& source, VectorElemental& elvec,
207 const CurrentFE& fe,
UInt iblock,
const Real& time)
const;
235 template<
typename MeshType,
typename DofType>
239 const ReferenceFE& refFE,
240 const QuadratureRule& quadRule,
256 template <
typename MatrixType,
typename VectorType>
266 Chrono chronoElemComp;
267 Chrono chronoAssembly;
270 Real coeffBeta, coeffDiv;
273 VectorElemental beta ( M_fe.nbNode, nDimensions );
276 for (
UInt iVol = 0; iVol <
M_mesh.numVolumes(); iVol++ )
278 chronoUpdate.start();
280 M_fe.updateFirstSecondDeriv ( M_mesh.volumeList ( iVol ) );
284 this->computeParameters (dt, iVol, state, beta, coeffBeta, coeffDiv);
287 chronoElemComp.start();
292 this->bgradu_bgradv (coeffBeta, beta, M_elMat, M_fe, 0, 0, nDimensions);
296 this->gradp_bgradv (coeffBeta, beta, M_elMat, M_fe);
300 stiff ( coeffBeta, M_elMat, M_fe, nDimensions, nDimensions );
304 this->lapu_bgradv (-coeffBeta * M_viscosity, beta, M_elMat, M_fe, 0, 0, nDimensions);
308 this->lapu_gradq (-coeffBeta * M_viscosity, M_elMat, M_fe);
312 stiff_div (coeffDiv, M_elMat, M_fe );
314 chronoElemComp.stop();
316 chronoAssembly.start();
317 for ( UInt iComp = 0; iComp <= nDimensions; ++iComp )
318 for ( UInt jComp = 0; jComp <= nDimensions; ++jComp )
320 assembleMatrix ( matrix, M_elMat, M_fe, M_dof,
322 iComp * M_dof.numTotalDof(), jComp * M_dof.numTotalDof() );
323 chronoAssembly.stop();
329 debugStream (7100) <<
" Updating of element done in " 330 << chronoUpdate.diffCumul() <<
"s." << std::endl;
331 debugStream (7100) <<
" Determination of parameters done in " 332 << chronoBeta.diffCumul() <<
"s." << std::endl;
333 debugStream (7100) <<
" Element computations done in " 334 << chronoElemComp.diffCumul() <<
"s." << std::endl;
336 << chronoAssembly.diffCumul() <<
"s." << std::endl;
343 template <
typename MatrixType,
typename VectorType>
353 Chrono chronoElemComp;
354 Chrono chronoAssembly;
360 VectorElemental beta ( M_fe.nbNode, nDimensions );
363 for (
UInt iVol = 0; iVol <
M_mesh.numVolumes(); iVol++ )
365 chronoUpdate.start();
367 M_fe.updateFirstSecondDeriv ( M_mesh.volumeList ( iVol ) );
371 this->computeParameters (dt, iVol, state, beta, coeffBeta, coeffDiv);
374 chronoElemComp.start();
379 this->bgradu_bgradv (coeffBeta, beta, M_elMat, M_fe, 0, 0, nDimensions);
389 chronoElemComp.stop();
391 chronoAssembly.start();
395 assemb_mat ( matrix, M_elMat, M_fe, M_dof, iComp, jComp );
397 chronoAssembly.stop();
403 debugStream (7100) <<
" Updating of element done in " 404 << chronoUpdate.diffCumul() <<
"s." << std::endl;
405 debugStream (7100) <<
" Determination of parameters done in " 406 << chronoBeta.diffCumul() <<
"s." << std::endl;
407 debugStream (7100) <<
" Element computations done in " 408 << chronoElemComp.diffCumul() <<
"s." << std::endl;
410 << chronoAssembly.diffCumul() <<
"s." << std::endl;
416 template <
typename VectorType,
typename SourceType >
418 const SourceType& source,
const Real& time)
428 Chrono chronoElemComp;
429 Chrono chronoAssembly;
432 Real coeffBeta, coeffDiv;
435 VectorElemental beta ( M_fe.nbNode, nDimensions );
438 for (
UInt iVol = 0; iVol <
M_mesh.numVolumes(); iVol++ )
440 chronoUpdate.start();
442 M_fe.updateFirstDeriv ( M_mesh.volumeList ( iVol ) );
447 this->computeParameters (dt, iVol, state, beta, coeffBeta, coeffDiv);
450 chronoElemComp.start();
455 this->f_bgradv (coeffBeta, source, beta, M_elVec, M_fe, 0, time);
459 this->f_gradq (coeffBeta, source, M_elVec, M_fe, nDimensions, time);
460 chronoElemComp.stop();
462 chronoAssembly.start();
466 assembleVector ( vector, M_elVec, M_fe, M_dof, iComp, iComp * M_dof.numTotalDof() );
468 chronoAssembly.stop();
474 debugStream (7100) <<
" Updating of element done in " 475 << chronoUpdate.diffCumul() <<
"s." << std::endl;
476 debugStream (7100) <<
" Determination of parameters done in " 477 << chronoBeta.diffCumul() <<
"s." << std::endl;
478 debugStream (7100) <<
" Element computations done in " 479 << chronoElemComp.diffCumul() <<
"s." << std::endl;
481 << chronoAssembly.diffCumul() <<
"s." << std::endl;
486 template<
typename MeshType,
typename DofType>
489 output <<
"StabilizationSD::showMe() " << std::endl;
490 output <<
"Fluid Viscosity: " << M_viscosity << std::endl;
491 output <<
"Stabilization coefficient SUPG/SD: " << M_gammaBeta << std::endl;
492 output <<
"Stabilization coefficient div grad: " << M_gammaDiv << std::endl;
494 M_dof.showMe (output);
501 template<
typename VectorType>
503 VectorElemental& beta,
Real& coeffBeta,
Real& coeffDiv)
const 511 hK = M_fe.diameter();
517 for ( UInt iNode = 0; iNode < M_fe.nbNode; ++iNode )
519 for ( UInt iCoor = 0; iCoor < M_fe.nbLocalCoor(); ++iCoor )
521 UInt ig = M_dof.localToGlobalMap ( iVol, iNode ) + iCoor * nDof;
522 beta.vec() [ iCoor * M_fe.nbNode + iNode ] = state[ig];
527 Real bmax = std::fabs ( beta.vec() [ 0 ] );
528 for ( UInt l = 1; l < UInt ( M_fe.nbLocalCoor() *M_fe.nbNode ); ++l )
530 if ( bmax < std::fabs ( beta.vec() [ l ] ) )
532 bmax = std::fabs ( beta.vec() [ l ] );
536 coeffBeta =
M_gammaBeta / std::sqrt ( 4. / ( dt * dt)
537 + 4.*bmax * bmax / hK2
544 template<
typename MeshType,
typename DofType>
546 MatrixElemental& elmat,
const CurrentFE& fe)
const 549 "advection_grad matrix needs at least the first derivatives");
551 MatrixElemental::matrix_type v (fe.nbLocalCoor(), fe.nbQuadPt() );
556 for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
558 for (UInt icoor (0); icoor < fe.nbLocalCoor(); ++icoor)
560 VectorElemental::vector_view velicoor = vel.block (icoor);
562 for (UInt k (0); k < fe.nbNode; ++k)
564 v (icoor, ig) += velicoor (k) * fe.phi (k, ig);
569 for (UInt ic (0); ic < fe.nbLocalCoor(); ++ic)
571 MatrixElemental::matrix_view mat_ic3 = elmat.block (ic, fe.nbLocalCoor() );
572 MatrixElemental::matrix_view mat_3ic = elmat.block (fe.nbLocalCoor(), ic);
573 for (UInt i = 0; i < fe.nbNode; ++i)
575 for (UInt j = 0; j < fe.nbNode; ++j)
578 for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
579 for (UInt jcoor (0); jcoor < fe.nbLocalCoor(); ++jcoor)
581 s += fe.phiDer (j, ic, ig) * v (jcoor, ig) * fe.phiDer (i, jcoor, ig) * fe.weightDet (ig);
583 mat_ic3 (i, j) += coef * s;
584 mat_3ic (j, i) += coef * s;
591 template<
typename MeshType,
typename DofType>
596 "advection (vect) matrix needs at least the first derivatives");
599 MatrixElemental::matrix_type mat_tmp (fe.nbNode, fe.nbNode);
600 MatrixElemental::matrix_type v ( fe.nbLocalCoor(), fe.nbQuadPt() );
605 for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
607 for (UInt icoor (0); icoor < fe.nbLocalCoor(); ++icoor)
609 VectorElemental::vector_view velicoor = vel.block (icoor);
611 for (UInt k (0); k < fe.nbNode; k++)
613 v (icoor, ig) += velicoor (k) * fe.phi (k, ig);
618 for (UInt i (0); i < fe.nbNode; ++i)
620 for (UInt j (0); j < fe.nbNode; ++j)
624 for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
625 for (UInt icoor (0); icoor < fe.nbLocalCoor(); ++icoor)
626 for (UInt jcoor (0); jcoor < fe.nbLocalCoor(); ++jcoor)
628 s += fe.phiDer (i, jcoor, ig) * v (jcoor, ig) * v (icoor, ig) * fe.phiDer (j, icoor, ig) * fe.weightDet (ig);
630 mat_tmp (i, j) = coef * s;
635 for (
UInt icomp (0); icomp < nb; icomp++)
637 MatrixElemental::matrix_view mat_icomp = elmat.block (iblock + icomp, jblock + icomp);
638 for (UInt i (0); i < fe.nbDiag(); ++i)
640 for (UInt j (0); j < fe.nbDiag(); ++j)
642 mat_icomp (i, j) += mat_tmp (i, j);
650 template<
typename MeshType,
typename DofType>
657 "lapu_bgradv matrix needs first derivatives");
660 "lapu_bgradv matrix needs second derivatives");
662 MatrixElemental::matrix_type mat_tmp (fe.nbNode, fe.nbNode);
663 MatrixElemental::matrix_type v ( fe.nbLocalCoor(), fe.nbQuadPt() );
668 for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
670 for (UInt icoor (0); icoor < fe.nbLocalCoor(); ++icoor)
672 VectorElemental::vector_view velicoor = vel.block (icoor);
674 for (UInt k (0); k < fe.nbNode; ++k)
676 v (icoor, ig) += velicoor (k) * fe.phi (k, ig);
683 for (UInt i (0); i < fe.nbNode; ++i)
685 for (UInt j (0); j < fe.nbNode; ++j)
688 for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
689 for (UInt icoor (0); icoor < fe.nbLocalCoor(); ++icoor)
690 for (UInt jcoor (0); jcoor < fe.nbLocalCoor(); ++jcoor)
692 s += fe.phiDer2 (j, icoor, icoor, ig) * v (jcoor, ig) * fe.phiDer (i, jcoor, ig) * fe.weightDet (ig);
694 mat_tmp (i, j) = coef * s;
699 for (
UInt icomp (0); icomp < nb; ++icomp)
701 MatrixElemental::matrix_view mat_icomp = elmat.block (iblock + icomp, jblock + icomp);
702 for (UInt i (0); i < fe.nbDiag(); ++i)
704 for (UInt j (0); j < fe.nbDiag(); ++j)
706 mat_icomp (i, j) += mat_tmp (i, j);
713 template<
typename MeshType,
typename DofType>
718 "lapu_gradq matrix needs first derivatives");
721 "lapu_gradq matrix needs second derivatives");
725 for (UInt jc (0); jc < fe.nbLocalCoor(); ++jc)
727 MatrixElemental::matrix_view mat_view = elmat.block (fe.nbLocalCoor(), jc);
728 for (UInt i (0); i < fe.nbNode; ++i)
730 for (UInt j (0); j < fe.nbNode; ++j)
734 for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
735 for (UInt jcoor (0); jcoor < fe.nbLocalCoor(); ++jcoor)
737 s += fe.phiDer2 (j, jcoor, jcoor, ig) * fe.phiDer (i, jc, ig) * fe.weightDet (ig);
739 mat_view (i, j) += coef * s;
747 template<
typename SourceType>
749 VectorElemental& elvec,
const CurrentFE& fe,
UInt iblock,
const Real& time)
const 753 "f_bgradv vector needs at least the first derivatives");
755 MatrixElemental::matrix_type v (fe.nbLocalCoor(), fe.nbQuadPt() );
760 for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
762 for (UInt icoor (0); icoor < fe.nbLocalCoor(); ++icoor)
764 VectorElemental::vector_view velicoor = vel.block (icoor);
766 for (UInt k (0); k < fe.nbNode; ++k)
768 v (icoor, ig) += velicoor (k) * fe.phi (k, ig);
774 for (UInt ic (0); ic < fe.nbLocalCoor(); ++ic)
776 VectorElemental::vector_view vec_ic = elvec.block (ic + iblock);
777 for (UInt i (0); i < fe.nbNode; ++i)
780 for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
781 for (UInt jcoor (0); jcoor < fe.nbLocalCoor(); ++jcoor)
782 s += source (time, fe.quadPt (ig, 0), fe.quadPt (ig, 1), fe.quadPt (ig, 2), ic + 1)
783 * fe.phiDer (i, jcoor, ig) * v (jcoor, ig) * fe.weightDet (ig);
784 vec_ic (i) += coef * s;
792 template<
typename SourceType>
797 "f_gradq vector needs at least the first derivatives");
801 VectorElemental::vector_view vec_ic = elvec.block (iblock);
802 for (UInt i (0); i < fe.nbNode; ++i)
805 for (UInt ig = 0; ig < fe.nbQuadPt(); ++ig)
806 for (UInt jcoor (0); jcoor < fe.nbLocalCoor(); ++jcoor)
807 s += source (time, fe.quadPt (ig, 0), fe.quadPt (ig, 1), fe.quadPt (ig, 2), jcoor + 1)
808 * fe.phiDer (i, jcoor, ig) * fe.weightDet (ig);
809 vec_ic (i) += coef * s;
CurrentFE M_fe
current fe for the assembling of stabilization terms.
double operator()(const char *VarName, const double &Default) const
Real M_viscosity
fluid viscosity
Real M_gammaBeta
Stabilization coefficient of .
void bgradu_bgradv(const Real &coef, VectorElemental &vel, ElemMat &elmat, const CurrentFE &fe, UInt iblock, UInt jblock, UInt nb) const
Evaluate the varf .
void applySUPG(const Real dt, MatrixType &matrix, const VectorType &state)
compute the SUPG stabilization terms and add them into the monolithic N.S. matrix ...
StabilizationSD()
Default Constructor.
void f_gradq(const Real &coef, SourceType &source, VectorElemental &elvec, const CurrentFE &fe, UInt iblock, const Real &time) const
Evaluate the varf .
void lapu_gradq(const Real &coef, MatrixElemental &elmat, const CurrentFE &fe) const
Evaluate the varf .
void applySD(const Real dt, MatrixType &matrix, const VectorType &state)
compute the SD stabilization terms and add them into the Momentum matrix
void applyRHS(const Real dt, VectorType &vector, const VectorType &state, const SourceType &source, const Real &time)
compute the SUPG stabilization terms and add them into the right and side
void updateInverseJacobian(const UInt &iQuadPt)
StabilizationSD(const GetPot &dataFile, const mesh_Type &mesh, const dof_Type &dof, const ReferenceFE &refFE, const QuadratureRule &quadRule, Real viscosity)
Constructor.
void gradp_bgradv(const Real &coef, VectorElemental &vel, MatrixElemental &elmat, const CurrentFE &fe) const
Evaluate the varf .
Real M_gammaDiv
Stabilization coefficient of .
VectorElemental M_elVec
Elementary Vector for assembling the stabilization terms.
void setGammaBeta(const Real &gammaBeta)
Set the stabilization parameter for .
void lapu_bgradv(const Real &coef, VectorElemental &vel, MatrixElemental &elmat, const CurrentFE &fe, UInt iblock, UInt jblock, UInt nb) const
Evaluate the varf .
void computeParameters(const Real dt, const UInt iVol, const VectorType &state, VectorElemental &beta, Real &coeffBeta, Real &coeffDiv) const
Compute the stabilization coefficients for each element.
void setGammaDiv(const Real &gammaDiv)
Set the stabilization parameter for .
const dof_Type & M_dof
the dof object
StabilizationSD(const StabilizationSD< mesh_Type, dof_Type > &original)
Copy Constructor.
double Real
Generic real data.
const UInt nDimensions(NDIM)
virtual ~StabilizationSD()
~Destructor
const mesh_Type & M_mesh
the mesh object
void showMe(std::ostream &output=std::cout) const
Display class informations.
void f_bgradv(const Real &coef, SourceType &source, VectorElemental &vel, VectorElemental &elvec, const CurrentFE &fe, UInt iblock, const Real &time) const
Evaluate the varf .
MatrixElemental M_elMat
Elementary Matrix for assembling the stabilization terms.
uint32_type UInt
generic unsigned integer (used mainly for addressing)