LifeV
test_bdf_2d.cpp
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
29  @brief
30 
31  @author Umberto Villa <uvilla@emory.edu>
32  @date 14-04-2010
33  */
34 
35 
36 // ===================================================
37 //! Includes
38 // ===================================================
39 
40 #include <Epetra_ConfigDefs.h>
41 #ifdef EPETRA_MPI
42 #include <mpi.h>
43 #include <Epetra_MpiComm.h>
44 #else
45 #include <Epetra_SerialComm.h>
46 #endif
47 
48 
49 #include <lifev/core/array/MapEpetra.hpp>
50 #include <lifev/core/mesh/MeshData.hpp>
51 #include <lifev/core/mesh/MeshPartitioner.hpp>
52 #include <lifev/core/algorithm/SolverAztecOO.hpp>
53 #include <lifev/core/algorithm/LinearSolver.hpp>
54 #include <lifev/core/fem/Assembly.hpp>
55 #include <lifev/core/fem/AssemblyElemental.hpp>
56 #include <lifev/core/fem/BCManage.hpp>
57 #include <lifev/core/fem/FESpace.hpp>
58 #include <lifev/core/fem/TimeAdvanceBDFVariableStep.hpp>
59 
60 #include <lifev/core/algorithm/PreconditionerIfpack.hpp>
61 #include <lifev/core/algorithm/PreconditionerML.hpp>
62 
63 #include <lifev/core/filter/Exporter.hpp>
64 #include <lifev/core/filter/ExporterEmpty.hpp>
65 #include <lifev/core/filter/ExporterEnsight.hpp>
66 #include <lifev/core/filter/ExporterHDF5.hpp>
67 #include <lifev/core/mesh/RegionMesh.hpp>
68 
69 #include <Teuchos_ParameterList.hpp>
70 #include <Teuchos_XMLParameterListHelpers.hpp>
71 #include <Teuchos_RCP.hpp>
72 
73 #include "ud_functions.hpp"
74 #include "test_bdf.hpp"
75 
76 // ===================================================
77 //! Namespaces & Define
78 // ===================================================
79 using namespace LifeV;
80 
81 const int TOP = 3;
82 const int BOTTOM = 1;
83 const int LEFT = 1;
84 const int RIGHT = 2;
85 
87 
88 const std::string discretization_section = "space_discretization";
89 
90 // ===================================================
91 //! Private Members
92 // ===================================================
93 struct test_bdf::Private
94 {
96  {
97  }
98 
99  typedef std::function < Real (Real const&, Real const&, Real const&,
100  Real const&, ID const&) > fct_type;
101 
104  Real errorNorm;
105 
106 };
107 
108 // ===================================================
109 //! Constructors
110 // ===================================================
111 test_bdf::test_bdf (int argc, char** argv) :
112  Members (new Private)
113 {
114  GetPot command_line (argc, argv);
115  const std::string data_file_name =
116  command_line.follow ("data_2d", 2, "-f", "--file");
117  GetPot dataFile (data_file_name);
118 
119  Members->data_file_name = data_file_name;
120 
121 #ifdef EPETRA_MPI
122  std::cout << "Epetra Initialization" << std::endl;
123  Members->comm.reset (new Epetra_MpiComm (MPI_COMM_WORLD) );
124 #else
125  Members->comm.reset (new Epetra_SerialComm() );
126 #endif
127 }
128 
129 // ===================================================
130 //! Methods
131 // ===================================================
132 void test_bdf::run()
133 {
134  //Useful typedef
135  typedef VectorEpetra vector_Type;
136  typedef std::shared_ptr<vector_Type> vectorPtr_Type;
137 
138  typedef LifeV::Preconditioner basePrec_Type;
139  typedef std::shared_ptr<basePrec_Type> basePrecPtr_Type;
140  typedef LifeV::PreconditionerIfpack prec_Type;
141  typedef std::shared_ptr<prec_Type> precPtr_Type;
142 
143  // Reading from data file
144  GetPot dataFile (Members->data_file_name.c_str() );
145  // Select Pid(0) to write on std output
146  bool verbose = (Members->comm->MyPID() == 0);
147 
148  if (verbose)
149  {
150  std::cout << "The BDF Solver" << std::flush;
151  }
152 
153  //the forcing term
154  SourceFct_2d sf;
155  // the boundary conditions
156  BCFunctionBase g_Ess (AnalyticalSol_2d::u);
157 
158  BCHandler bc;
159 
160  bc.addBC ("Top", TOP, Essential, Full, g_Ess, 1);
161  bc.addBC ("Bottom", BOTTOM, Essential, Full, g_Ess, 1);
162  bc.addBC ("Left", LEFT, Essential, Full, g_Ess, 1);
163  bc.addBC ("Right", RIGHT, Essential, Full, g_Ess, 1);
164  //=============================================================================
165  //Mesh stuff
166  Members->comm->Barrier();
167  MeshData meshData (dataFile, ("bdf/" + discretization_section).c_str() );
168  std::shared_ptr<regionMesh> fullMeshPtr ( new regionMesh ( Members->comm ) );
169  readMesh (*fullMeshPtr, meshData);
170  std::shared_ptr<regionMesh> meshPtr;
171  {
172  MeshPartitioner<regionMesh> meshPart ( fullMeshPtr, Members->comm );
173  meshPtr = meshPart.meshPartition();
174  }
175 
176  //=============================================================================
177  //finite element space of the solution
178  std::shared_ptr<FESpace<regionMesh, MapEpetra> > feSpacePtr (
179  new FESpace<regionMesh, MapEpetra> (
180  meshPtr, dataFile ( ("bdf/" + discretization_section + "/order").c_str(), "P2"), 1, Members->comm) );
181 
182  if (verbose)
183  std::cout << " Number of unknowns : "
184  << feSpacePtr->map().map (Unique)->NumGlobalElements()
185  << std::endl;
186 
187  bc.bcUpdate (* (feSpacePtr->mesh() ), feSpacePtr->feBd(), feSpacePtr->dof() );
188 
189  //=============================================================================
190  //Fe Matrices and vectors
191  MatrixElemental elmat (feSpacePtr->fe().nbFEDof(), 1, 1); //local matrix
192  MatrixEpetra<double> matM (feSpacePtr->map() ); //mass matrix
193  std::shared_ptr<MatrixEpetra<double> > matA_ptr (
194  new MatrixEpetra<double> (feSpacePtr->map() ) ); //stiff matrix
195  std::shared_ptr<VectorEpetra> u ( new VectorEpetra (feSpacePtr->map(), Unique) ); // solution vector
196  std::shared_ptr<VectorEpetra> f ( new VectorEpetra (feSpacePtr->map(), Unique) ); // forcing term vector
197 
198  LifeChrono chrono;
199  //Assembling Matrix M
200  Members->comm->Barrier();
201  chrono.start();
202  for (UInt iVol = 0; iVol < feSpacePtr->mesh()->numElements(); iVol++)
203  {
204  feSpacePtr->fe().updateJac (feSpacePtr->mesh()->element (iVol) );
205  elmat.zero();
206  AssemblyElemental::mass (1., elmat, feSpacePtr->fe(), 0, 0);
207  assembleMatrix (matM, elmat, feSpacePtr->fe(), feSpacePtr->fe(), feSpacePtr->dof(),
208  feSpacePtr->dof(), 0, 0, 0, 0);
209  }
210  matM.globalAssemble();
211  Members->comm->Barrier();
212  chrono.stop();
213  if (verbose)
214  std::cout << "\n \n -- Mass matrix assembling time = " << chrono.diff() << std::endl
215  << std::endl;
216 
217  // ==========================================
218  // Definition of the time integration stuff
219  // ==========================================
220  Real Tfin = dataFile ("bdf/endtime", 10.0);
221  Real delta_t = dataFile ("bdf/timestep", 0.5);
222  Real t0 = 1.;
223  UInt ord_bdf = dataFile ("bdf/order", 3);
224  TimeAdvanceBDFVariableStep<VectorEpetra> bdf;
225  bdf.setup (ord_bdf);
226 
227  //Initialization
228  bdf.setInitialCondition<Real (*) (Real, Real, Real, Real, UInt), FESpace<regionMesh, MapEpetra> > (AnalyticalSol_2d::u, *u, *feSpacePtr, t0, delta_t);
229 
230  if (verbose)
231  {
232  bdf.showMe();
233  }
234  Members->comm->Barrier();
235 
236  //===================================================
237  // post processing setup
238  std::shared_ptr<Exporter<regionMesh> > exporter;
239  std::string const exporterType = dataFile ( "exporter/type", "hdf5");
240 
241 #ifdef HAVE_HDF5
242  if (exporterType.compare ("hdf5") == 0)
243  {
244  exporter.reset ( new ExporterHDF5<regionMesh > ( dataFile, "bdf_test" ) );
245  }
246  else
247 #endif
248  {
249  if (exporterType.compare ("none") == 0)
250  {
251  exporter.reset ( new ExporterEmpty<regionMesh > ( dataFile, meshPtr, "bdf_test", Members->comm->MyPID() ) );
252  }
253  else
254  {
255  exporter.reset ( new ExporterEnsight<regionMesh > ( dataFile, meshPtr, "bdf_test", Members->comm->MyPID() ) );
256  }
257  }
258 
259  exporter->setPostDir ( "./" );
260  exporter->setMeshProcId ( meshPtr, Members->comm->MyPID() );
261 
262  std::shared_ptr<VectorEpetra> u_display_ptr (new VectorEpetra (
263  feSpacePtr->map(), exporter->mapType() ) );
264  exporter->addVariable (ExporterData<regionMesh >::ScalarField, "u", feSpacePtr,
265  u_display_ptr, UInt (0) );
266  *u_display_ptr = *u;
267  exporter->postProcess (0);
268 
269 
270  //===================================================
271  //Definition of the linear solver
272  LinearSolver linearSolver;
273 
274  Teuchos::RCP< Teuchos::ParameterList > belosList = Teuchos::rcp ( new Teuchos::ParameterList );
275  belosList = Teuchos::getParametersFromXmlFile ( "SolverParamList.xml" );
276 
277  linearSolver.setCommunicator ( Members->comm );
278  linearSolver.setParameters ( *belosList );
279 
280  prec_Type* precRawPtr;
281  basePrecPtr_Type precPtr;
282  precRawPtr = new prec_Type;
283  precRawPtr->setDataFromGetPot ( dataFile, "bdf/prec" );
284  precPtr.reset ( precRawPtr );
285  linearSolver.setPreconditioner ( precPtr );
286 
287  //===================================================
288  // TIME LOOP
289  //===================================================
290 
291  matA_ptr.reset (new MatrixEpetra<double> (feSpacePtr->map() ) );
292 
293  for (Real t = t0 + delta_t; t <= Tfin; t += delta_t)
294  {
295  Members->comm->Barrier();
296  if (verbose)
297  {
298  std::cout << "Now we are at time " << t << std::endl;
299  }
300 
301  chrono.start();
302  //Assemble A
303  Real coeff = bdf.coefficientDerivative (0) / delta_t;
304  Real visc = nu (t);
305  Real s = sigma (t);
306  for (UInt i = 0; i < feSpacePtr->mesh()->numElements(); i++)
307  {
308  feSpacePtr->fe().updateFirstDerivQuadPt (feSpacePtr->mesh()->element (i) );
309  elmat.zero();
310  AssemblyElemental::mass (coeff + s, elmat, feSpacePtr->fe() );
311  AssemblyElemental::stiff (visc, elmat, feSpacePtr->fe() );
312  assembleMatrix (*matA_ptr, elmat, feSpacePtr->fe(), feSpacePtr->fe(),
313  feSpacePtr->dof(), feSpacePtr->dof(), 0, 0, 0, 0);
314  }
315 
316  chrono.stop();
317  if (verbose)
318  {
319  std::cout << "A has been constructed in " << chrono.diff() << "s." << std::endl;
320  }
321 
322  // Handling of the right hand side
323  *f = (matM * bdf.rhsContribution() ); //f = M*\sum_{i=1}^{orderBdf} \alpha_i u_{n-i}
324  feSpacePtr->l2ScalarProduct (sf, *f, t); //f +=\int_\Omega{ volumeForces *v dV}
325  Members->comm->Barrier();
326 
327  // Treatment of the Boundary conditions
328  if (verbose)
329  {
330  std::cout << "*** BC Management: " << std::endl;
331  }
332  Real tgv = 1.;
333  chrono.start();
334  bcManage (*matA_ptr, *f, *feSpacePtr->mesh(), feSpacePtr->dof(), bc, feSpacePtr->feBd(), tgv, t);
335  matA_ptr->globalAssemble();
336  chrono.stop();
337  if (verbose)
338  {
339  std::cout << chrono.diff() << "s." << std::endl;
340  }
341 
342  //Set Up the linear system
343  Members->comm->Barrier();
344  chrono.start();
345 
346  linearSolver.setOperator ( matA_ptr );
347  linearSolver.setReusePreconditioner (false);
348  linearSolver.setRightHandSide ( f );
349  linearSolver.solve ( u );
350 
351  bdf.shiftRight (*u);
352 
353  if (verbose)
354  std::cout << "*** Solution computed in " << chrono.diff() << "s."
355  << std::endl;
356 
357  Members->comm->Barrier();
358 
359  // Error in the L2
360  vector_Type uComputed (*u, Repeated);
361 
362  Real L2_Error, L2_RelError;
363 
364  L2_Error = feSpacePtr->l2Error (AnalyticalSol_2d::u, uComputed, t, &L2_RelError);
365 
366  if (verbose)
367  std::cout << "Error Norm L2: " << L2_Error
368  << "\nRelative Error Norm L2: " << L2_RelError << std::endl;
369 
370  Members->errorNorm = L2_Error;
371 
372  //transfer the solution at time t.
373  *u_display_ptr = *u;
374 
375  matA_ptr->zero();
376 
377  exporter->postProcess (t);
378  }
379 
380 
381 
382 }
383 
384 bool test_bdf::check()
385 {
386  // Reading from data file
387  GetPot dataFile (Members->data_file_name.c_str() );
388  return Members->errorNorm < dataFile ("errorNorms/l2Error", -10e10);
389 }
const int BOTTOM
Definition: test_bdf_2d.cpp:82
double operator()(const char *VarName, const double &Default) const
Definition: GetPot.hpp:2033
const std::string discretization_section
Definition: test_bdf_2d.cpp:88
GetPot(const int argc_, char **argv_, const char *FieldSeparator=0x0)
Definition: GetPot.hpp:507
std::shared_ptr< Epetra_Comm > comm
Definition: test_bdf.cpp:105
const int LEFT
Definition: test_bdf_2d.cpp:83
TimeAdvanceBDFVariableStep - Backward differencing formula time discretization with non-uniform time ...
std::shared_ptr< std::vector< Int > > M_isOnProc
const int RIGHT
Definition: test_bdf_2d.cpp:84
MeshData - class for handling spatial discretization.
Definition: MeshData.hpp:72
double Real
Generic real data.
Definition: LifeV.hpp:175
std::string data_file_name
Definition: test_bdf.cpp:104
int operator()(const char *VarName, int Default) const
Definition: GetPot.hpp:2021
std::function< Real(Real const &, Real const &, Real const &, Real const &, ID const &) > fct_type
const int TOP
Definition: test_bdf_2d.cpp:81
RegionMesh< LinearTriangle > regionMesh
Definition: test_bdf_2d.cpp:86
Private Members.
Definition: test_bdf.cpp:95
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191