LifeV
HyperbolicSolver.hpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
5  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
6  This file is part of LifeV.
7  LifeV is free software; you can redistribute it and/or modify
8  it under the terms of the GNU Lesser General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11  LifeV is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  Lesser General Public License for more details.
15  You should have received a copy of the GNU Lesser General Public License
16  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
17 *******************************************************************************
18 */
19 //@HEADER
20 /*!
21  * @file
22  * @brief Solver class for hyperbolic scalar equations.
23  *
24  *
25  * @date 30-09-2010
26  *
27  * @author Alessio Fumagalli <alessio.fumagalli@mail.polimi.it>
28  * @author Michel Kern <michel.kern@inria.fr>
29  *
30  * @contributor
31  *
32  * @mantainer Alessio Fumagalli <alessio.fumagalli@mail.polimi.it>
33  *
34  */
35 
36 #ifndef _HYPERBOLICSOLVER_H_
37 #define _HYPERBOLICSOLVER_H_ 1
38 
39 #include <Epetra_LAPACK.h>
40 
41 #include <lifev/core/algorithm/SolverAztecOO.hpp>
42 
43 #include <lifev/core/fem/AssemblyElemental.hpp>
44 #include <lifev/core/fem/BCManage.hpp>
45 #include <lifev/core/fem/HyperbolicFluxNumerical.hpp>
46 
47 #include <lifev/core/solver/HyperbolicData.hpp>
48 
49 // LifeV namespace.
50 namespace LifeV
51 {
52 //! HyperbolicSolver Implements an hyperbolic solver.
53 /*!
54 
55  @author Alessio Fumagalli <alessio.fumagalli@mail.polimi.it>
56  @author Michel Kern <michel.kern@inria.fr>
57  @see For applications related to two-phase flow see \cite Fumagalli2011a
58 
59  This class implements an hyperbolic solver.
60  <br>
61  <br>
62  This class solves a general hyperbolic scalar equation in the conservative form: find \f$ u \in V \f$ such that
63  \f[
64  \left\{
65  \begin{array}{l l}
66  \displaystyle \frac{\partial u}{\partial t} + \nabla \cdot \mathbf{F} ( u ) = f & \mathrm{ in } \quad \Omega \times [0,T]\,, \\
67  u(\mathbf{x}, 0) = u_0( \mathbf{x} ) & \mathrm{ in } \quad \Omega\,,
68  \end{array}
69  \right.
70  \f]
71  with prescribed inflow boundary conditions on \f$ \partial \Omega_{in} \f$
72  \f[
73  u = \bar{u} \quad \mathrm{ on } \quad \partial \Omega_{in}\,.
74  \f]
75  We have \f$ \partial \Omega = \partial \Omega_{in} \cap \partial \Omega_{out}\f$.
76  Since this solver can be used coupled with an elliptic solver, the Neumann and Robin boundary
77  are trated as outflow boundary, while the Dirichlet boundary is automatically dedived into
78  outflow boundary and inflow boundary. This choice is due to the general form of the flux function \f$ \mathbf{F} \f$.
79  <br>
80  Introducing a conforming triangolation \f$ \mathcal{T}_h = \cup K \f$ of the domain \f$ \Omega \f$, we write the strong formulation
81  in a weak form in each element \f$ K \f$: give a test function \f$ v \f$
82  \f[
83  \displaystyle \int_K \frac{\partial u}{\partial t} v + \int_K \nabla \cdot \mathbf{F} ( u ) v = \int_K f v \quad \forall v \in V\,,
84  \f]
85  integrating by part the integral with the divergence we find
86  \f[
87  \displaystyle \int_K \frac{\partial u}{\partial t} v - \int_K \mathbf{F} ( u ) \cdot \nabla v +
88  \int_{\partial K} \mathbf{F}( u ) \cdot \mathbf{n} v = \int_K f v \quad \forall v \in V\,.
89  \f]
90  The Discontinuous Galerkin formulation of the problem is given approximating \f$ V \f$ with a discontinuous finite element
91  space \f$ V_h \f$, while the numerical flux at the boundaries is approximated using a suitable numerical scheme obtaining
92  \f$ \hat{\mathbf{F}} \f$.
93  <br>
94  Summing on all the elements \f$ K \in \mathcal{T}_h \f$, the finite element problem reads: find \f$ u_h \in V_h \f$ such that
95  \f[
96  \displaystyle \sum_{K \in \mathcal{T}_h} \int_K \frac{\partial u_h}{\partial t} v_h - \int_K \mathbf{F} ( u_h ) \cdot \nabla v_h +
97  \int_{\partial K} \hat{\mathbf{F}}( u_h ) \cdot \mathbf{n} v_h = \sum_{K \in \mathcal{T}_h} \int_K f v_h \quad \forall v_h \in V_h\,.
98  \f]
99  Each time step \f$ \Delta t^n \f$ and mesh size \f$ h \f$ are coupled via \f$ CFL \f$ condition at step \f$ n \f$, using esplicit Euler
100  scheme we find the relation between \f$ \Delta t^n \f$ and \f$ h \f$ satysfying \f$ CFL^n < 1 \f$.
101  <br>
102  A general form for computing the \f$ CFL \f$ condition is
103  \f[
104  CFL^n = \displaystyle \sup_{e \in \partial K, K \in \mathcal{T}_h } \Delta t^n \frac{\vert e \vert}{ \vert K \vert }
105  \Vert \mathbf{F}^\prime \cdot \mathbf{n}_{e, K} \Vert_{L^\infty(a_0, b_0) }\,,
106  \f]
107  where
108  \f[
109  \begin{array}{l}
110  \displaystyle a_0 = \inf_{x \in \Omega, t \in (t^{n-1}, t^n), y \in \partial \Omega_{in} } \left\{ u_0, \bar{u}(y, t) \right\}\,, \\
111  \displaystyle b_0 = \sup_{x \in \Omega, t \in (t^{n-1}, t^n), y \in \partial \Omega_{in} } \left\{ u_0, \bar{u}(y, t) \right\}\,.
112  \end{array}
113  \f]
114  Local approximation for the flux function \f$ \hat{\mathbf{F}} \cdot \mathbf{n} \f$ and the computation of
115  \f$ \mathbf{F}^\prime \cdot \mathbf{n}_{e, K} \f$ are in NumericalFlux.hpp file.
116  @note The implementation is given just for lowest order discontinuous finite elements.
117  @todo Implement the forcing term \f$ f \f$ and implement high order finite elements.
118  @todo When we will pass to Trilinos >= 10.6 use Epetra wrapper for LAPACK functions.
119 */
120 template < typename Mesh,
121  typename SolverType = LifeV::SolverAztecOO >
123 {
124 
125 public:
126 
127  //! @name Public Types
128  //@{
129 
130  typedef std::function < Real ( const Real&, const Real&, const Real&,
131  const Real&, const UInt& ) >
133 
134  typedef HyperbolicData< Mesh > data_Type;
135 
136  typedef BCHandler bchandler_Type;
138 
139  typedef typename SolverType::vector_type vector_Type;
141 
144 
145  typedef AbstractNumericalFlux<Mesh, SolverType> flux_Type;
147 
152 
153  //@}
154 
155  //! @name Constructors & Destructor
156  //@{
157 
158  //! Full constructor for the class.
159  /*!
160  @param dataFile Data for the problem.
161  @param fESpace Discontinuous finite element space.
162  @param bcHandler Boundary conditions for the problem.
163  @param comm Shared pointer of the Epetra communicator.
164  */
165  HyperbolicSolver ( const data_Type& dataFile,
166  FESpace<Mesh, MapEpetra>& fESpace,
167  MapEpetra& ghostMap,
168  bchandler_Type& bcHandler,
169  commPtr_Type& comm );
170 
171  //! Constructor for the class without the definition of the boundary handler.
172  /*!
173  @param dataFile Data for the problem.
174  @param fESpace Discontinuous finite element space.
175  @param comm Shared pointer of the Epetra communicator.
176  */
177  HyperbolicSolver ( const data_Type& dataFile,
178  FESpace<Mesh, MapEpetra>& fESpace,
179  MapEpetra& ghostMap,
180  commPtr_Type& comm );
181 
182  //! Virtual destructor.
183  virtual ~HyperbolicSolver ();
184 
185  //@}
186 
187  //! @name Methods
188  //@{
189 
190  //! Setup the local mass matrices.
191  void setup ();
192 
193  //! Solve one time step of the hyperbolic problem.
194  void solveOneTimeStep();
195 
196  //! Compute the global CFL condition.
197  Real CFL();
198 
199  //@}
200 
201  //! @name Set Methos
202  //@{
203 
204  //! Set the initial solution for the computation.
205  /*!
206  Compute the initial solution as an interpolation on the analytical
207  initial solution.
208  @param initialSolution The initial solution function.
209  */
210  void setInitialSolution ( const Function_Type& initialSolution );
211 
212  //! Set the boundary conditions.
213  /*!
214  @param bcHandler Boundary condition handler for the problem.
215  */
216  inline void setBoundaryCondition ( bchandler_Type& bcHandler )
217  {
218  M_BCh.reset ( new bchandler_Type ( bcHandler ) );
219  M_setBC = true;
220  }
221 
222  //! Set the source term.
223  /*
224  The default setted source term is \f$ f( \mathbf{x}, t ) \equiv 0 \f$.
225  @param source Source term for the problem.
226  */
227  inline void setSourceTerm ( const Function_Type& source )
228  {
229  M_source = source;
230  }
231 
232  //! Set the mass function term.
233  /*!
234  It does not depend on time. The default setted source term is \f$ \phi( \mathbf{x} ) \equiv 1 \f$.
235  @param mass Mass term for the problem.
236  */
237  inline void setMassTerm ( const Function_Type& mass )
238  {
239  M_mass = mass;
240  }
241 
242  //! Set the numerical flux.
243  /*!
244  @param flux The numerical flux class
245  */
246  inline void setNumericalFlux ( const flux_Type& flux )
247  {
248  typedef GodunovNumericalFlux<Mesh, SolverType> godunov_Type;
249  M_numericalFlux.reset ( new godunov_Type ( dynamic_cast<const godunov_Type&> ( flux ) ) );
250  }
251 
252  //! Set the solution vector.
253  /*!
254  @param solution Constant vector_type reference of the solution.
255  */
256  inline void setSolution ( const vectorPtr_Type& solution )
257  {
258  // Set both the final step solution and beginning step solution.
259  *M_u = *solution;
260  *M_uOld = *solution;
261  }
262 
263  //@}
264 
265  //! @name Get Methods
266  //@{
267 
268  //! Returns the solution vector.
269  /*!
270  @return Constant vector_type reference of the solution vector.
271  */
272  inline const vectorPtr_Type& solution () const
273  {
274  return M_u;
275  }
276 
277  //! Return if the bounday conditions is setted or not.
278  /*!
279  @return Constant boolean with value true if the boundary condition is setted,
280  false otherwise
281  */
282  inline bool isBoundaryConditionSet() const
283  {
284  return M_setBC;
285  }
286 
287  //! Return the boundary conditions handler.
288  /*!
289  @return Reference of boundary conditions handler.
290  */
292  {
293  return M_BCh;
294  }
295 
296 
297  //! Return the Epetra local map.
298  /*!
299  @return Constant MapEpetra reference of the problem.
300  */
301  inline MapEpetra const& map () const
302  {
303  return M_localMap;
304  }
305 
306  //! Return the displayer.
307  /*!
308  Useful for parallel print in programs.
309  @return Constant reference of the displayer of the problem.
310  */
311  inline Displayer const& getDisplayer() const
312  {
313  return M_displayer;
314  }
315 
316  //! Returns displayer.
317  /*!
318  @return Reference of the displayer of the problem.
319  */
320  inline Displayer& getDisplayer()
321  {
322  return M_displayer;
323  }
324 
325  //@}
326 
327 protected:
328 
329  //! @name Protected Methods
330  //@{
331 
332  //! Reconstruct locally the solution.
333  void localReconstruct ( const UInt& Elem );
334 
335  //! Compute the local contribute
336  void localEvolve ( const UInt& iElem );
337 
338  //! Apply the flux limiters locally.
339  void localAverage ( const UInt& iElem );
340 
341  //@}
342 
343  //! MPI process identifier.
344  const UInt M_me;
345 
346  //! Local map.
347  MapEpetra M_localMap;
348 
349  //! Parallel displayer
350  Displayer M_displayer;
351 
352  //! Data for Darcy solvers.
354 
355  //! Source function.
357 
358  //! Mass function, it does not depend on time
360 
361  //! Bondary conditions handler.
363 
364  //! Flag if the boundary conditions are setted or not.
365  bool M_setBC;
366 
367  //! Function of initial solution.
369 
370  //! Class of type AbstractNumericalFlux for local flux computations.
372 
373  //! Finite element space.
374  FESpace<Mesh, MapEpetra>& M_FESpace;
375 
376  //! Right hand side.
378 
379  //! Solution at current time step.
381 
382  //! Solution at previous time step.
384 
385  //! Computed numerical flux.
387 
388  //! Auxiliary vector for local fluxes.
389  VectorElemental M_localFlux;
390 
391  //! Vector of all local mass matrices, possibly with mass function.
393 
394 private:
395 
396  //! @name Private Constructors
397  //@{
398 
399  //! Inhibited copy constructor.
400  HyperbolicSolver ( const HyperbolicSolver<Mesh, SolverType>& );
401 
402  //@}
403 
404  //! @name Private Operators
405  //@{
406 
407  //! Inhibited assign operator.
408  HyperbolicSolver& operator= ( const HyperbolicSolver<Mesh, SolverType>& );
409 
410  //@}
411 
412 }; // class HyperbolicSolver
413 
414 // ===================================================
415 // Constructors & Destructor
416 // ===================================================
417 
418 // Complete constructor.
419 template< typename Mesh, typename SolverType >
420 HyperbolicSolver< Mesh, SolverType >::
421 HyperbolicSolver ( const data_Type& dataFile,
422  FESpace<Mesh, MapEpetra>& fESpace,
423  MapEpetra& ghostMap,
424  bchandler_Type& bcHandler,
425  commPtr_Type& comm ) :
426  // Parallel stuff.
427  M_me ( comm->MyPID() ),
428  M_localMap ( fESpace.map() ),
429  M_displayer ( comm ),
430  // Data of the problem.
431  M_data ( dataFile ),
432  M_source ( NULL ),
433  M_mass ( NULL ),
434  M_BCh ( new bchandler_Type ( bcHandler ) ),
435  M_setBC ( true ),
436  M_initialSolution ( NULL ),
437  M_numericalFlux ( ),
438  // Finite element spaces.
439  M_FESpace ( fESpace ),
440  // Algebraic stuff.
441  M_rhs ( new vector_Type ( M_localMap ) ),
442  M_u ( new vector_Type ( M_FESpace.map(), Repeated ) ),
443  M_uOld ( new vector_Type ( ghostMap, Repeated ) ),
444  M_globalFlux ( new vector_Type ( ghostMap, Repeated ) ),
445  // Local matrices and vectors.
446  M_localFlux ( M_FESpace.refFE().nbDof(), 1 ),
447  M_elmatMass ( )
448 {
449 
450  M_elmatMass.reserve ( M_FESpace.mesh()->numElements() );
451 
452 } // Constructor
453 
454 
455 // Constructor without boundary condition handler.
456 template< typename Mesh, typename SolverType >
457 HyperbolicSolver< Mesh, SolverType >::
458 HyperbolicSolver ( const data_Type& dataFile,
459  FESpace<Mesh, MapEpetra>& fESpace,
460  MapEpetra& ghostMap,
461  commPtr_Type& comm ) :
462  // Parallel stuff.
463  M_me ( comm->MyPID() ),
464  M_localMap ( fESpace.map() ),
465  M_displayer ( comm ),
466  // Data of the problem.
467  M_data ( dataFile ),
468  M_source ( NULL ),
469  M_mass ( NULL ),
470  M_BCh ( ),
471  M_setBC ( false ),
472  M_initialSolution ( NULL ),
473  M_numericalFlux ( ),
474  // Finite element spaces.
475  M_FESpace ( fESpace ),
476  // Algebraic stuff.
477  M_rhs ( new vector_Type ( M_localMap ) ),
478  M_u ( new vector_Type ( M_FESpace.map(), Repeated ) ),
479  M_uOld ( new vector_Type ( ghostMap, Repeated ) ),
480  M_globalFlux ( new vector_Type ( ghostMap, Repeated ) ),
481  // Local matrices and vectors.
482  M_localFlux ( M_FESpace.refFE().nbDof(), 1 ),
483  M_elmatMass ( )
484 {
485 
486  M_elmatMass.reserve ( M_FESpace.mesh()->numElements() );
487 
488 } // Constructor
489 
490 // Virtual destructor.
491 template< typename Mesh, typename SolverType >
492 HyperbolicSolver< Mesh, SolverType >::
494 {
495 
496 } // Destructor
497 
498 // ===================================================
499 // Methods
500 // ===================================================
501 
502 // Setup the local mass matrices.
503 template<typename Mesh, typename SolverType>
504 void
505 HyperbolicSolver<Mesh, SolverType>::
507 {
508  // LAPACK wrapper of Epetra
509  Epetra_LAPACK lapack;
510 
511  // Flags for LAPACK routines.
512  Int INFO[1] = {0};
513  Int NB = M_FESpace.refFE().nbDof();
514 
515  // Parameter that indicate the Lower storage of matrices.
516  char param_L = 'L';
517 
518  // Total number of elements.
519  UInt meshNumberOfElements = M_FESpace.mesh()->numElements();
520 
521  // Vector of interpolation of the mass function
522  vector_Type vectorMass ( M_FESpace.map(), Repeated );
523 
524  // If the mass function is given take it, otherwise use one as a value.
525  if ( M_mass != NULL )
526  {
527  // Interpolate the mass function on the finite element space.
528  M_FESpace.interpolate ( M_mass, vectorMass );
529  }
530  else
531  {
532  vectorMass = 1.;
533  }
534 
535  // For each element it creates the mass matrix and factorize it using Cholesky.
536  for ( UInt iElem (0); iElem < meshNumberOfElements; ++iElem )
537  {
538 
539  // Update the element.
540  M_FESpace.fe().update ( M_FESpace.mesh()->element ( iElem),
542 
543  // Local mass matrix
544  MatrixElemental matElem (M_FESpace.refFE().nbDof(), 1, 1);
545  matElem.zero();
546 
547  // Compute the mass matrix for the current element
548  VectorElemental massValue ( M_FESpace.refFE().nbDof(), 1 );
549  extract_vec ( vectorMass, massValue, M_FESpace.refFE(), M_FESpace.dof(), iElem, 0 );
550  // TODO: this works only for P0
551  AssemblyElemental::mass ( massValue[ 0 ], matElem, M_FESpace.fe(), 0, 0);
552 
553  /* Put in M the matrix L and L^T, where L and L^T is the Cholesky factorization of M.
554  For more details see http://www.netlib.org/lapack/double/dpotrf.f */
555  lapack.POTRF ( param_L, NB, matElem.mat(), NB, INFO );
556  ASSERT_PRE ( !INFO[0], "Lapack factorization of M is not achieved." );
557 
558  // Save the local mass matrix in the global vector of mass matrices
559  M_elmatMass.push_back ( matElem );
560 
561  }
562 
563  //make sure mesh facets are updated
564  if (! M_FESpace.mesh()->hasLocalFacets() )
565  {
566  M_FESpace.mesh()->updateElementFacets();
567  }
568 
569 } // setup
570 
571 // Solve one time step of the hyperbolic problem.
572 template< typename Mesh, typename SolverType >
573 void
574 HyperbolicSolver< Mesh, SolverType >::
576 {
577  // Total number of elements in the mesh
578  const UInt meshNumberOfElements ( M_FESpace.mesh()->numElements() );
579 
580  // Loop on all the elements to perform the fluxes
581  for ( UInt iElem (0); iElem < meshNumberOfElements; ++iElem )
582  {
583 
584  // Update the property of the current element
585  M_FESpace.fe().update ( M_FESpace.mesh()->element (iElem), UPDATE_QUAD_NODES | UPDATE_WDET);
586 
587  // Reconstruct step of the current element
588  localReconstruct ( iElem );
589 
590  // Evolve step of the current element
591  localEvolve ( iElem );
592 
593  // Put the total flux of the current element in the global vector of fluxes
594  assembleVector ( *M_globalFlux,
595  M_FESpace.fe().currentLocalId(),
596  M_localFlux,
597  M_FESpace.refFE().nbDof(),
598  M_FESpace.dof(), 0 );
599 
600  // Average step of the current element
601  localAverage ( iElem );
602 
603  }
604 
605  // Assemble the global hybrid vector.
606  M_globalFlux->globalAssemble();
607 
608  // alternative: instead of modifying M_globalFlux.map, we can make a local copy with the correct map
609  // // this is needed since M_uOld.map != M_globalFlux.map
610  // vector_Type fluxCopy ( M_uOld->map() );
611  // fluxCopy = *M_globalFlux;
612 
613  // Update the value of the solution
614  (*M_u) = (*M_uOld) - M_data.dataTime()->timeStep() * (*M_globalFlux);
615  // *M_u = *M_uOld - M_data.dataTime()->timeStep() * fluxCopy;
616 
617  // Clean the vector of fluxes
618  M_globalFlux.reset ( new vector_Type ( M_uOld->map(), Repeated ) );
619 
620  // Update the solution at previous time step
621  *M_uOld = *M_u;
622 
623 } // solveOneStep
624 
625 // Compute the global CFL condition.
626 template< typename Mesh, typename SolverType >
627 Real
628 HyperbolicSolver< Mesh, SolverType >::
630 {
631  // Total number of elements in the mesh
632  const UInt meshNumberOfElements ( M_FESpace.mesh()->numElements() );
633 
634  // The local value for the CFL condition, without the time step
635  Real localCFL (0.), localCFLOld ( - 1. );
636 
637  // Loop on all the elements to perform the fluxes
638  for ( UInt iElem (0); iElem < meshNumberOfElements; ++iElem )
639  {
640  // Update the property of the current element
641  M_FESpace.fe().update ( M_FESpace.mesh()->element (iElem), UPDATE_QUAD_NODES | UPDATE_WDET);
642 
643  // Volumetric measure of the current element
644  const Real K ( M_FESpace.fe().measure() );
645 
646  // Loop on the faces of the element iElem and compute the local contribution
647  for ( UInt iFace (0); iFace < M_FESpace.mesh()->numLocalFaces(); ++iFace )
648  {
649 
650  const UInt iGlobalFace ( M_FESpace.mesh()->localFacetId ( iElem, iFace ) );
651 
652  // Update the normal vector of the current face in each quadrature point
653  M_FESpace.feBd().update ( M_FESpace.mesh()->boundaryFacet ( iGlobalFace ), UPDATE_W_ROOT_DET_METRIC | UPDATE_NORMALS | UPDATE_QUAD_NODES );
654 
655  // Take the left element to the face, see regionMesh for the meaning of left element
656  const UInt leftElement ( M_FESpace.mesh()->faceElement ( iGlobalFace, 0 ) );
657 
658  // Take the right element to the face, see regionMesh for the meaning of right element
659  const UInt rightElement ( M_FESpace.mesh()->faceElement ( iGlobalFace, 1 ) );
660 
661  // Solution in the left element
662  VectorElemental leftValue ( M_FESpace.refFE().nbDof(), 1 );
663 
664  // Solution in the right element
665  VectorElemental rightValue ( M_FESpace.refFE().nbDof(), 1 );
666 
667  // Extract the solution in the current element, now is the leftElement
668  extract_vec ( *M_uOld,
669  leftValue,
670  M_FESpace.refFE(),
671  M_FESpace.dof(),
672  leftElement , 0 );
673 
674  if ( !Flag::testOneSet ( M_FESpace.mesh()->face ( iGlobalFace ).flag(),
676  {
677  // Extract the solution in the current element, now is the leftElement
678  extract_vec ( *M_uOld,
679  rightValue,
680  M_FESpace.refFE(),
681  M_FESpace.dof(),
682  rightElement , 0 );
683  }
684  else if ( Flag::testOneSet ( M_FESpace.mesh()->face ( iGlobalFace ).flag(), EntityFlags::SUBDOMAIN_INTERFACE ) )
685  {
686  // TODO: this works only for P0 elements
687  // but extract_vec works only with lids while RightElement is a gid
688  rightValue[ 0 ] = (*M_uOld) [ rightElement ];
689  }
690  else // Flag::testOneSet ( M_FESpace.mesh()->face ( iGlobalFace ).flag(), PHYSICAL_BOUNDARY )
691  {
692  rightValue = leftValue;
693  }
694 
695  // Area of the current face
696  const Real e ( M_FESpace.feBd().measure() );
697 
698  // Loop on all the quadrature points
699  for ( UInt ig (0); ig < M_FESpace.feBd().nbQuadPt(); ++ig )
700  {
701 
702  // Cuurent quadrature point
703  KN<Real> quadPoint (3);
704  // normal vector
705  KN<Real> normal (3);
706 
707  for (UInt icoor (0); icoor < 3; ++icoor)
708  {
709  quadPoint (icoor) = M_FESpace.feBd().quadPt ( ig, icoor );
710  normal (icoor) = M_FESpace.feBd().normal ( icoor, ig ) ;
711  }
712 
713  // Compute the local CFL without the time step
714  localCFL = e / K * M_numericalFlux->normInfinity ( leftValue[0],
715  rightValue[0],
716  normal,
717  iElem,
718  M_data.dataTime()->time(),
719  quadPoint (0),
720  quadPoint (1),
721  quadPoint (2) );
722 
723  // Select the maximum between the old CFL condition and the new CFL condition
724  if ( localCFL > localCFLOld )
725  {
726  localCFLOld = localCFL;
727  }
728  else
729  {
730  localCFL = localCFLOld;
731  }
732 
733  }
734 
735  }
736 
737  }
738 
739 
740  // Compute the time step according to CLF for the current process
741  Real timeStepLocal[] = { M_data.getCFLRelaxParameter() / localCFLOld };
742  Real timeStepGlobal[] = { 0. };
743 
744  // Compute the minimum of the computed time step for all the processes
745  M_displayer.comm()->MinAll ( timeStepLocal, timeStepGlobal, 1 );
746 
747  // Return the computed value
748  return *timeStepGlobal;
749 
750 } //CFL
751 
752 // ===================================================
753 // Set Methods
754 // ===================================================
755 
756 // Set the initial solution for the computation.
757 template< typename Mesh, typename SolverType >
758 void
759 HyperbolicSolver< Mesh, SolverType >::
760 setInitialSolution ( const Function_Type& initialSolution )
761 {
762 
763  // interpolation must be done on a Unique map
764  vector_Type uUnique ( M_u->map(), Unique );
765 
766  // Interpolate the initial solution.
767  M_FESpace.interpolate ( initialSolution,
768  uUnique,
769  M_data.dataTime()->initialTime() );
770 
771  // Update the solutions
772  *M_uOld = uUnique;
773  *M_u = uUnique;
774 } // setInitialSolution
775 
776 // ===================================================
777 // Protected Methods
778 // ===================================================
779 
780 // Reconstruct locally the solution.
781 template< typename Mesh, typename SolverType >
782 void
783 HyperbolicSolver< Mesh, SolverType >::
784 localReconstruct ( const UInt& /* iElem */ )
785 {
786 
787  ;
788 
789 } // localReconstruct
790 
791 // Compute the local contribute
792 template< typename Mesh, typename SolverType >
793 void
794 HyperbolicSolver< Mesh, SolverType >::
795 localEvolve ( const UInt& iElem )
796 {
797 
798  // LAPACK wrapper of Epetra
799  Epetra_LAPACK lapack;
800 
801  // Flags for LAPACK routines.
802  Int INFO[1] = { 0 };
803  Int NB = M_FESpace.refFE().nbDof();
804 
805  // Parameter that indicate the Lower storage of matrices.
806  char param_L = 'L';
807  char param_N = 'N';
808 
809  // Paramater that indicate the Transpose of matrices.
810  char param_T = 'T';
811 
812  // Numbers of columns of the right hand side := 1.
813  Int NBRHS = 1;
814 
815  // Clean the local flux
816  M_localFlux.zero();
817 
818  // Check if the boundary conditions were updated.
819  if ( !M_BCh->bcUpdateDone() )
820  {
821  // Update the boundary conditions handler. We use the finite element of the boundary of the dual variable.
822  M_BCh->bcUpdate ( *M_FESpace.mesh(), M_FESpace.feBd(), M_FESpace.dof() );
823  }
824 
825  // Loop on the faces of the element iElem and compute the local contribution
826  for ( UInt iFace (0); iFace < M_FESpace.mesh()->numLocalFaces(); ++iFace )
827  {
828  // Id mapping
829  const UInt iGlobalFace ( M_FESpace.mesh()->localFacetId ( iElem, iFace ) );
830 
831  // Take the left element to the face, see regionMesh for the meaning of left element
832  const UInt leftElement ( M_FESpace.mesh()->faceElement ( iGlobalFace, 0 ) );
833 
834  // Take the right element to the face, see regionMesh for the meaning of right element
835  const UInt rightElement ( M_FESpace.mesh()->faceElement ( iGlobalFace, 1 ) );
836 
837  // Update the normal vector of the current face in each quadrature point
838  M_FESpace.feBd().update ( M_FESpace.mesh()->boundaryFacet ( iGlobalFace ), UPDATE_W_ROOT_DET_METRIC | UPDATE_NORMALS | UPDATE_QUAD_NODES );
839 
840  // Local flux of a face times the integration weight
841  VectorElemental localFaceFluxWeight ( M_FESpace.refFE().nbDof(), 1 );
842 
843  // Solution in the left element
844  VectorElemental leftValue ( M_FESpace.refFE().nbDof(), 1 );
845 
846  // Solution in the right element
847  VectorElemental rightValue ( M_FESpace.refFE().nbDof(), 1 );
848 
849  // Extract the solution in the current element, now is the leftElement
850  extract_vec ( *M_uOld,
851  leftValue,
852  M_FESpace.refFE(),
853  M_FESpace.dof(),
854  leftElement , 0 );
855 
856  // Check if the current face is a boundary face, that is rightElement == NotAnId
858  {
859  // Extract the solution in the current element, now is the leftElement
860  extract_vec ( *M_uOld,
861  rightValue,
862  M_FESpace.refFE(),
863  M_FESpace.dof(),
864  rightElement , 0 );
865  }
866  else if ( Flag::testOneSet ( M_FESpace.mesh()->face ( iGlobalFace ).flag(), EntityFlags::SUBDOMAIN_INTERFACE ) )
867  {
868  // TODO: this works only for P0 elements
869  // rightValue[ 0 ] = M_ghostDataMap[ iGlobalFace ];
870  rightValue[ 0 ] = (*M_uOld) [ rightElement ];
871  }
872  else // Flag::testOneSet ( M_FESpace.mesh()->face ( iGlobalFace ).flag(), PHYSICAL_BOUNDARY )
873  {
874 
875  // Clean the value of the right element
876  rightValue.zero();
877 
878  // Take the boundary marker for the current boundary face
879  const ID faceMarker ( M_FESpace.mesh()->boundaryFacet ( iGlobalFace ).markerID() );
880 
881  // Take the corrispective boundary function
882  const BCBase& bcBase ( M_BCh->findBCWithFlag ( faceMarker ) );
883 
884  // Check if the bounday condition is of type Essential, useful for operator splitting strategies
885  if ( bcBase.type() == Essential )
886  {
887 
888  // Loop on all the quadrature points
889  for ( UInt ig (0); ig < M_FESpace.feBd().nbQuadPt(); ++ig)
890  {
891 
892  // Current quadrature point
893  KN<Real> quadPoint (3);
894 
895  // normal vector
896  KN<Real> normal (3);
897 
898  for (UInt icoor (0); icoor < 3; ++icoor)
899  {
900  quadPoint (icoor) = M_FESpace.feBd().quadPt ( ig, icoor );
901  normal (icoor) = M_FESpace.feBd().normal ( icoor, ig ) ;
902  }
903 
904  // Compute the boundary contribution
905  rightValue[0] = bcBase ( M_data.dataTime()->time(), quadPoint (0), quadPoint (1), quadPoint (2), 0 );
906 
907  const Real localFaceFlux = M_numericalFlux->firstDerivativePhysicalFluxDotNormal ( normal,
908  iElem,
909  M_data.dataTime()->time(),
910  quadPoint (0),
911  quadPoint (1),
912  quadPoint (2),
913  rightValue[ 0 ] );
914  // Update the local flux of the current face with the quadrature weight
915  localFaceFluxWeight[0] += localFaceFlux * M_FESpace.feBd().wRootDetMetric ( ig );
916  }
917 
918  }
919  else
920  {
921  /* If the boundary flag is not Essential then is automatically an outflow boundary.
922  We impose to localFaceFluxWeight a positive value. */
923  localFaceFluxWeight[0] = 1.;
924  }
925 
926  // It is an outflow face, we use a ghost cell
927  if ( localFaceFluxWeight[0] > 1e-4 )
928  {
929  rightValue = leftValue;
930  }
931 
932  // Clean the localFaceFluxWeight
933  localFaceFluxWeight.zero();
934 
935  }
936 
937  // Clean the localFaceFluxWeight
938  localFaceFluxWeight.zero();
939 
940  // Loop on all the quadrature points
941  for ( UInt ig (0); ig < M_FESpace.feBd().nbQuadPt(); ++ig )
942  {
943 
944  // Current quadrature point
945  KN<Real> quadPoint (3);
946 
947  // normal vector
948  KN<Real> normal (3);
949 
950  for (UInt icoor (0); icoor < 3; ++icoor)
951  {
952  quadPoint (icoor) = M_FESpace.feBd().quadPt ( ig, icoor );
953  normal (icoor) = M_FESpace.feBd().normal ( icoor, ig ) ;
954  }
955 
956  // If the normal is orientated inward, we change its sign and swap the left and value of the solution
957  if ( iElem == rightElement )
958  {
959  normal *= -1.;
960  std::swap ( leftValue, rightValue );
961  }
962 
963  const Real localFaceFlux = (*M_numericalFlux) ( leftValue[ 0 ],
964  rightValue[ 0 ],
965  normal,
966  iElem,
967  M_data.dataTime()->time(),
968  quadPoint (0),
969  quadPoint (1),
970  quadPoint (2) );
971 
972  // Update the local flux of the current face with the quadrature weight
973  localFaceFluxWeight[0] += localFaceFlux * M_FESpace.feBd().wRootDetMetric ( ig );
974 
975  }
976 
977  /* Put in localFlux the vector L^{-1} * localFlux
978  For more details see http://www.netlib.org/lapack/lapack-3.1.1/SRC/dtrtrs.f */
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." );
981 
982  /* Put in localFlux the vector L^{-T} * localFlux
983  For more details see http://www.netlib.org/lapack/lapack-3.1.1/SRC/dtrtrs.f */
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." );
986 
987  // Add to the local flux the local flux of the current face
988  M_localFlux += localFaceFluxWeight;
989 
990  }
991 
992 } // localEvolve
993 
994 // Apply the flux limiters locally.
995 template< typename Mesh, typename SolverType >
996 void
997 HyperbolicSolver< Mesh, SolverType >::
998 localAverage ( const UInt& /* iElem */ )
999 {
1000 
1001  ;
1002 
1003 } // localAverage
1004 
1005 } // namespace LifeV
1006 
1007 #endif /*_HYPERBOLICSOLVER_H_ */
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!
Definition: FESpace.hpp:78
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.
Definition: LifeV.hpp:188
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
Definition: LifeV.hpp:216
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&#39;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.
#define ASSERT_PRE(X, A)
Definition: LifeAssert.hpp:96
uint32_type ID
IDs.
Definition: LifeV.hpp:194
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.
Definition: LifeV.hpp:175
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
Definition: LifeV.hpp:203
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
Definition: MeshEntity.hpp:48
GodunovNumericalFlux Gives an implementation for Godunov solver for hyperbolic&#39;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)
Definition: LifeV.hpp:191
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.