36 #ifndef _HYPERBOLICSOLVER_H_ 37 #define _HYPERBOLICSOLVER_H_ 1
39 #include <Epetra_LAPACK.h> 41 #include <lifev/core/algorithm/SolverAztecOO.hpp> 43 #include <lifev/core/fem/AssemblyElemental.hpp> 44 #include <lifev/core/fem/BCManage.hpp> 45 #include <lifev/core/fem/HyperbolicFluxNumerical.hpp> 47 #include <lifev/core/solver/HyperbolicData.hpp> 120 template <
typename Mesh,
121 typename SolverType = LifeV::SolverAztecOO >
166 FESpace<Mesh, MapEpetra>& fESpace,
178 FESpace<Mesh, MapEpetra>& fESpace,
218 M_BCh.reset (
new bchandler_Type ( bcHandler ) );
249 M_numericalFlux.reset (
new godunov_Type (
dynamic_cast<
const godunov_Type&> ( flux ) ) );
301 inline MapEpetra
const&
map ()
const 419 template<
typename Mesh,
typename SolverType >
422 FESpace<Mesh, MapEpetra>& fESpace,
450 M_elmatMass.reserve ( M_FESpace.mesh()->numElements() );
456 template<
typename Mesh,
typename SolverType >
459 FESpace<Mesh, MapEpetra>& fESpace,
486 M_elmatMass.reserve ( M_FESpace.mesh()->numElements() );
491 template<
typename Mesh,
typename SolverType >
503 template<
typename Mesh,
typename SolverType>
509 Epetra_LAPACK lapack;
536 for (
UInt iElem (0); iElem < meshNumberOfElements; ++iElem )
544 MatrixElemental matElem (
M_FESpace.refFE().nbDof(), 1, 1);
548 VectorElemental massValue (
M_FESpace.refFE().nbDof(), 1 );
551 AssemblyElemental::mass ( massValue[ 0 ], matElem, M_FESpace.fe(), 0, 0);
555 lapack.POTRF ( param_L, NB, matElem.mat(), NB, INFO );
556 ASSERT_PRE ( !INFO[0],
"Lapack factorization of M is not achieved." );
559 M_elmatMass.push_back ( matElem );
564 if (!
M_FESpace.mesh()->hasLocalFacets() )
572 template<
typename Mesh,
typename SolverType >
578 const UInt meshNumberOfElements (
M_FESpace.mesh()->numElements() );
581 for (
UInt iElem (0); iElem < meshNumberOfElements; ++iElem )
594 assembleVector ( *M_globalFlux,
595 M_FESpace.fe().currentLocalId(),
597 M_FESpace.refFE().nbDof(),
598 M_FESpace.dof(), 0 );
606 M_globalFlux->globalAssemble();
614 (*M_u) = (*M_uOld) - M_data.dataTime()->timeStep() * (*M_globalFlux);
618 M_globalFlux.reset (
new vector_Type ( M_uOld->map(), Repeated ) );
626 template<
typename Mesh,
typename SolverType >
632 const UInt meshNumberOfElements (
M_FESpace.mesh()->numElements() );
635 Real localCFL (0.), localCFLOld ( - 1. );
638 for (
UInt iElem (0); iElem < meshNumberOfElements; ++iElem )
647 for (
UInt iFace (0); iFace <
M_FESpace.mesh()->numLocalFaces(); ++iFace )
650 const UInt iGlobalFace (
M_FESpace.mesh()->localFacetId ( iElem, iFace ) );
656 const UInt leftElement (
M_FESpace.mesh()->faceElement ( iGlobalFace, 0 ) );
659 const UInt rightElement (
M_FESpace.mesh()->faceElement ( iGlobalFace, 1 ) );
662 VectorElemental leftValue (
M_FESpace.refFE().nbDof(), 1 );
665 VectorElemental rightValue (
M_FESpace.refFE().nbDof(), 1 );
668 extract_vec ( *M_uOld,
678 extract_vec ( *M_uOld,
688 rightValue[ 0 ] = (*M_uOld) [ rightElement ];
692 rightValue = leftValue;
703 KN<
Real> quadPoint (3);
707 for (
UInt icoor (0); icoor < 3; ++icoor)
709 quadPoint (icoor) = M_FESpace.feBd().quadPt ( ig, icoor );
710 normal (icoor) = M_FESpace.feBd().normal ( icoor, ig ) ;
714 localCFL = e / K * M_numericalFlux->normInfinity ( leftValue[0],
718 M_data.dataTime()->time(),
724 if ( localCFL > localCFLOld )
726 localCFLOld = localCFL;
730 localCFL = localCFLOld;
741 Real timeStepLocal[] = {
M_data.getCFLRelaxParameter() / localCFLOld };
742 Real timeStepGlobal[] = { 0. };
745 M_displayer.comm()->MinAll ( timeStepLocal, timeStepGlobal, 1 );
748 return *timeStepGlobal;
757 template<
typename Mesh,
typename SolverType >
769 M_data.dataTime()->initialTime() );
781 template<
typename Mesh,
typename SolverType >
792 template<
typename Mesh,
typename SolverType >
799 Epetra_LAPACK lapack;
819 if ( !M_BCh->bcUpdateDone() )
822 M_BCh->bcUpdate ( *M_FESpace.mesh(), M_FESpace.feBd(), M_FESpace.dof() );
826 for (
UInt iFace (0); iFace <
M_FESpace.mesh()->numLocalFaces(); ++iFace )
829 const UInt iGlobalFace (
M_FESpace.mesh()->localFacetId ( iElem, iFace ) );
832 const UInt leftElement (
M_FESpace.mesh()->faceElement ( iGlobalFace, 0 ) );
835 const UInt rightElement (
M_FESpace.mesh()->faceElement ( iGlobalFace, 1 ) );
841 VectorElemental localFaceFluxWeight (
M_FESpace.refFE().nbDof(), 1 );
844 VectorElemental leftValue (
M_FESpace.refFE().nbDof(), 1 );
847 VectorElemental rightValue (
M_FESpace.refFE().nbDof(), 1 );
850 extract_vec ( *M_uOld,
860 extract_vec ( *M_uOld,
870 rightValue[ 0 ] = (*M_uOld) [ rightElement ];
879 const ID faceMarker (
M_FESpace.mesh()->boundaryFacet ( iGlobalFace ).markerID() );
882 const BCBase& bcBase ( M_BCh->findBCWithFlag ( faceMarker ) );
885 if ( bcBase.type() == Essential )
893 KN<
Real> quadPoint (3);
898 for (
UInt icoor (0); icoor < 3; ++icoor)
900 quadPoint (icoor) = M_FESpace.feBd().quadPt ( ig, icoor );
901 normal (icoor) = M_FESpace.feBd().normal ( icoor, ig ) ;
905 rightValue[0] = bcBase ( M_data.dataTime()->time(), quadPoint (0), quadPoint (1), quadPoint (2), 0 );
907 const Real localFaceFlux = M_numericalFlux->firstDerivativePhysicalFluxDotNormal ( normal,
909 M_data.dataTime()->time(),
915 localFaceFluxWeight[0] += localFaceFlux * M_FESpace.feBd().wRootDetMetric ( ig );
923 localFaceFluxWeight[0] = 1.;
927 if ( localFaceFluxWeight[0] > 1e-4 )
929 rightValue = leftValue;
933 localFaceFluxWeight.zero();
938 localFaceFluxWeight.zero();
945 KN<
Real> quadPoint (3);
950 for (
UInt icoor (0); icoor < 3; ++icoor)
952 quadPoint (icoor) = M_FESpace.feBd().quadPt ( ig, icoor );
953 normal (icoor) = M_FESpace.feBd().normal ( icoor, ig ) ;
957 if ( iElem == rightElement )
960 std::swap ( leftValue, rightValue );
963 const Real localFaceFlux = (*M_numericalFlux) ( leftValue[ 0 ],
967 M_data.dataTime()->time(),
973 localFaceFluxWeight[0] += localFaceFlux * M_FESpace.feBd().wRootDetMetric ( ig );
979 lapack.TRTRS ( param_L, param_N, param_N, NB, NBRHS, M_elmatMass[ iElem ].mat(), NB, localFaceFluxWeight, NB, INFO);
980 ASSERT_PRE ( !INFO[0],
"Lapack Computation M_elvecSource = LB^{-1} rhs is not achieved." );
984 lapack.TRTRS ( param_L, param_T, param_N, NB, NBRHS, M_elmatMass[ iElem ].mat(), NB, localFaceFluxWeight, NB, INFO);
985 ASSERT_PRE ( !INFO[0],
"Lapack Computation M_elvecSource = LB^{-1} rhs is not achieved." );
988 M_localFlux += localFaceFluxWeight;
995 template<
typename Mesh,
typename SolverType >
const flag_Type UPDATE_NORMALS(UPDATE_ONLY_NORMALS|UPDATE_ONLY_TANGENTS|UPDATE_ONLY_CELL_NODES)
std::map< UInt, ghostDataContainer_Type > buffer_Type
AbstractNumericalFlux< Mesh, SolverType > flux_Type
bool isBoundaryConditionSet() const
Return if the bounday conditions is setted or not.
std::shared_ptr< bchandler_Type > bchandlerPtr_Type
std::shared_ptr< vector_Type > vectorPtr_Type
bool M_setBC
Flag if the boundary conditions are setted or not.
FESpace - Short description here please!
virtual ~HyperbolicSolver()
Virtual destructor.
HyperbolicSolver(const data_Type &dataFile, FESpace< Mesh, MapEpetra > &fESpace, MapEpetra &ghostMap, commPtr_Type &comm)
Constructor for the class without the definition of the boundary handler.
HyperbolicSolver(const HyperbolicSolver< Mesh, SolverType > &)
Inhibited copy constructor.
void setNumericalFlux(const flux_Type &flux)
Set the numerical flux.
int32_type Int
Generic integer data.
Displayer & getDisplayer()
Returns displayer.
VectorElemental M_localFlux
Auxiliary vector for local fluxes.
void setInitialSolution(const Function_Type &initialSolution)
Set the initial solution for the computation.
const flag_Type UPDATE_WDET(UPDATE_ONLY_CELL_NODES|UPDATE_ONLY_DPHI_GEO_MAP|UPDATE_ONLY_JACOBIAN|UPDATE_ONLY_DET_JACOBIAN|UPDATE_ONLY_W_DET_JACOBIAN)
HyperbolicData< Mesh > data_Type
bool testOneSet(flag_Type const &inputFlag, flag_Type const &refFlag)
returns true if at least one flag set in refFlag is set in inputFlag
vectorPtr_Type M_rhs
Right hand side.
const vectorPtr_Type & solution() const
Returns the solution vector.
void localAverage(const UInt &iElem)
Apply the flux limiters locally.
AbstractNumericalFlux Gives a common interface for hyperbolic's flux function.
void setMassTerm(const Function_Type &mass)
Set the mass function term.
const flag_Type SUBDOMAIN_INTERFACE(0x04)
void setSolution(const vectorPtr_Type &solution)
Set the solution vector.
const data_Type & M_data
Data for Darcy solvers.
bchandlerPtr_Type & boundaryConditionHandler()
Return the boundary conditions handler.
MapEpetra M_localMap
Local map.
std::shared_ptr< flux_Type > fluxPtr_Type
const flag_Type PHYSICAL_BOUNDARY(0x01)
void localEvolve(const UInt &iElem)
Compute the local contribute.
std::shared_ptr< comm_Type > commPtr_Type
HyperbolicSolver Implements an hyperbolic solver.
std::vector< MatrixElemental > M_elmatMass
Vector of all local mass matrices, possibly with mass function.
std::map< ID, ghostData_Type > ghostDataMap_Type
vectorPtr_Type M_u
Solution at current time step.
double Real
Generic real data.
Displayer M_displayer
Parallel displayer.
Real CFL()
Compute the global CFL condition.
void solveOneTimeStep()
Solve one time step of the hyperbolic problem.
Function_Type M_source
Source function.
flag related free functions and functors
void localReconstruct(const UInt &Elem)
Reconstruct locally the solution.
const flag_Type UPDATE_W_ROOT_DET_METRIC(UPDATE_ONLY_W_ROOT_DET_METRIC|UPDATE_METRIC|UPDATE_ONLY_DET_METRIC)
available bit-flags for different geometric properties
GodunovNumericalFlux Gives an implementation for Godunov solver for hyperbolic's flux function...
const flag_Type UPDATE_QUAD_NODES(UPDATE_ONLY_CELL_NODES|UPDATE_ONLY_QUAD_NODES)
void setBoundaryCondition(bchandler_Type &bcHandler)
Set the boundary conditions.
MapEpetra const & map() const
Return the Epetra local map.
FESpace< Mesh, MapEpetra > & M_FESpace
Finite element space.
Displayer const & getDisplayer() const
Return the displayer.
fluxPtr_Type M_numericalFlux
Class of type AbstractNumericalFlux for local flux computations.
HyperbolicSolver & operator=(const HyperbolicSolver< Mesh, SolverType > &)
Inhibited assign operator.
bchandlerPtr_Type M_BCh
Bondary conditions handler.
vectorPtr_Type M_globalFlux
Computed numerical flux.
std::vector< ghostData_Type > ghostDataContainer_Type
const UInt M_me
MPI process identifier.
void setup()
Setup the local mass matrices.
std::function< Real(const Real &, const Real &, const Real &, const Real &, const UInt &) > Function_Type
void setSourceTerm(const Function_Type &source)
Set the source term.
Function_Type M_mass
Mass function, it does not depend on time.
SolverType::vector_type vector_Type
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Function_Type M_initialSolution
Function of initial solution.
HyperbolicSolver(const data_Type &dataFile, FESpace< Mesh, MapEpetra > &fESpace, MapEpetra &ghostMap, bchandler_Type &bcHandler, commPtr_Type &comm)
Full constructor for the class.
vectorPtr_Type M_uOld
Solution at previous time step.