LifeV
AssemblyPolicyNavierStokesNewton.hpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4 
5  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
6  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7 
8  This file is part of LifeV.
9 
10  LifeV is free software; you can redistribute it and/or modify
11  it under the terms of the GNU Lesser General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  LifeV is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 
27 /*!
28  @file AssemblyPolicyNavierStokesSemiImplicit class
29  @brief This class contains all the informations necessary
30  to assemble a Navier-Stokes problem using a
31  Newton scheme
32 
33  @author Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
34  @date 06-12-2012
35  */
36 
37 #ifndef ASSEMBLYPOLICYNAVIERSTOKESNEWTON_HPP
38 #define ASSEMBLYPOLICYNAVIERSTOKESNEWTON_HPP
39 
40 #include <iostream>
41 #include <string>
42 #include <boost/shared_ptr.hpp>
43 
44 
45 #ifdef HAVE_MPI
46 #include <Epetra_MpiComm.h>
47 #else
48 #include <Epetra_SerialComm.h>
49 #endif
50 
51 #include <Teuchos_ParameterList.hpp>
52 #include <Teuchos_XMLParameterListHelpers.hpp>
53 #include <Teuchos_RCP.hpp>
54 
55 
56 #include <lifev/core/LifeV.hpp>
57 #include <lifev/core/array/MatrixEpetra.hpp>
58 #include <lifev/core/array/VectorEpetra.hpp>
59 #include <lifev/core/util/Displayer.hpp>
60 #include <lifev/core/mesh/RegionMesh.hpp>
61 #include <lifev/core/fem/FESpace.hpp>
62 #include <lifev/core/fem/TimeAdvanceBDF.hpp>
63 #include <lifev/core/algorithm/Preconditioner.hpp>
64 #include <lifev/core/util/LifeChrono.hpp>
65 
66 #include <lifev/navier_stokes/algorithm/PreconditionerPCD.hpp>
67 #include <lifev/navier_stokes/solver/NavierStokesSolver/AssemblyPolicyStokes.hpp>
68 #include <lifev/navier_stokes/solver/NavierStokesSolver/NavierStokesProblem.hpp>
69 
70 
71 namespace LifeV
72 {
73 
74 template< class mesh_Type >
76 {
82  typedef MeshPartitioner< mesh_Type > meshPartitioner_Type;
85  typedef FESpace< mesh_Type, map_Type > fespace_Type;
87  typedef TimeAdvanceBDF<vector_Type> bdf_Type;
91 
92  enum { BDFOrder = 1 };
93 
94  void initAssembly ( Teuchos::ParameterList& list );
95 
96  void assembleSystem ( matrixPtr_Type systemMatrix,
97  vectorPtr_Type rhs,
98  vectorPtr_Type solution,
99  preconditionerPtr_Type preconditioner );
100 
102 
103  virtual Displayer displayer() = 0;
104  virtual fespacePtr_Type uFESpace() const = 0;
105  virtual fespacePtr_Type pFESpace() const = 0;
106  virtual NSProblemPtr_Type problem() const = 0;
107  virtual Real timestep() const = 0;
108  virtual bdfPtr_Type bdf() const = 0;
109 };
110 
111 template< class mesh_Type >
112 void
113 AssemblyPolicyNavierStokesNewton< mesh_Type >::initAssembly ( Teuchos::ParameterList& list )
114 {
115  AssemblyPolicyStokes< mesh_Type >::initAssembly ( list );
116 
117  LifeChrono assemblyChrono;
118  assemblyChrono.start();
119 
120  displayer().leaderPrint ( "Creating the mass matrix... " );
121  map_Type solutionMap ( uFESpace()->map() + pFESpace()->map() );
122  M_massMatrix.reset ( new matrix_Type ( solutionMap ) );
123  M_massMatrix->zero();
124  AssemblyPolicyStokes< mesh_Type >::M_assembler->addMass ( *M_massMatrix, 1.0 );
125  M_massMatrix->globalAssemble();
126  displayer().leaderPrint ( "done\n" );
127 
128  assemblyChrono.stop();
129  displayer().leaderPrintMax ("Matrices assembly time: ", assemblyChrono.diff(), " s.\n");
130 }
131 
132 template< class mesh_Type >
133 void
135  vectorPtr_Type rhs,
136  vectorPtr_Type solution,
137  preconditionerPtr_Type preconditioner )
138 {
139  AssemblyPolicyStokes< mesh_Type >::assembleSystem ( systemMatrix, rhs, solution, preconditioner );
140 
141  bdf()->updateRHSContribution ( timestep() );
142  *rhs += *M_massMatrix * bdf()->rhsContributionFirstDerivative();
143 
144  double alpha = bdf()->coefficientFirstDerivative ( 0 ) / timestep();
145  *systemMatrix += *M_massMatrix * alpha;
146 
147  vector_Type beta ( systemMatrix->map(), Repeated );
148  beta += *solution;
149  AssemblyPolicyStokes< mesh_Type >::M_assembler->addConvection ( *systemMatrix, 1.0, beta );
150  AssemblyPolicyStokes< mesh_Type >::M_assembler->addNewtonConvection ( *systemMatrix, beta );
151  AssemblyPolicyStokes< mesh_Type >::M_assembler->addConvectionRhs ( *rhs, 1.0, beta );
152 
153  if ( preconditioner->preconditionerType() == "PCD" )
154  {
155  PreconditionerPCD* pcdPtr = dynamic_cast<PreconditionerPCD*> ( preconditioner.get() );
156  pcdPtr->updateBeta ( beta );
157  }
158 }
159 
160 } // namespace LifeV
161 
162 #endif /* ASSEMBLYPOLICYNAVIERSTOKESSEMIIMPLICIT_HPP */
VectorEpetra - The Epetra Vector format Wrapper.
std::shared_ptr< NavierStokesProblem< mesh_Type > > NSProblemPtr_Type
PreconditionerPCD(std::shared_ptr< Epetra_Comm > comm=std::shared_ptr< Epetra_Comm >(new Epetra_MpiComm(MPI_COMM_WORLD)))
default constructor
void start()
Start the timer.
Definition: LifeChrono.hpp:93
void importFromHDF5(std::string const &fileName, std::string const &matrixName="matrix")
Read a matrix from a HDF5 (.h5) file.
void updateInverseJacobian(const UInt &iQuadPt)
virtual fespacePtr_Type pFESpace() const =0
Epetra_Import const & importer()
Getter for the Epetra_Import.
Definition: MapEpetra.cpp:394
Real diff()
Compute the difference in time between start and stop.
Definition: LifeChrono.hpp:111
double Real
Generic real data.
Definition: LifeV.hpp:175
Preconditioner - Abstract preconditioner class.
void stop()
Stop the timer.
Definition: LifeChrono.hpp:100
virtual bdfPtr_Type bdf() const =0
virtual fespacePtr_Type uFESpace() const =0
std::shared_ptr< preconditioner_Type > preconditionerPtr_Type
void updateBeta(const vector_Type &beta)
Update the vector beta of the convective term in Fp.
virtual NSProblemPtr_Type problem() const =0
Displayer - This class is used to display messages in parallel simulations.
Definition: Displayer.hpp:62
void assembleSystem(matrixPtr_Type systemMatrix, vectorPtr_Type rhs, vectorPtr_Type solution, preconditionerPtr_Type preconditioner)