LifeV
lifev/fsi/examples/benchmark_GreenshieldsWeller/main.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  * @include fluidstructure.dox
29  * @file
30  * @brief for testing the benchmark contained in the
31  *
32  * @date 2009-04-09
33  * @author Paolo Crosetto <crosetto@iacspc70.epfl.ch>
34  *
35  * @maintainer Paolo Crosetto <crosetto@iacspc70.epfl.ch>
36  *
37  * Monolithic problem. Features:
38  * - fullMonolithic (CE):
39  * -# solution with exact Newton (semiImplicit = false, useShapeDerivatives = true, domainVelImplicit = false, convectiveTermDer = false)
40  * -# solution with quasi Newton (semiImplicit = false, useShapeDerivatives = false, domainVelImplicit = false, convectiveTermDer = false)
41  * -# preconditioner choice: see the classes Monolithic and fullMonolithic
42  * - Monolithic (GCE):
43  * -# solution extrapolating the fluid domain (semiImplicit = true, useShapeDerivatives = false, domainVelImplicit = false, convectiveTermDer = false)
44  * -# preconditioner choice: see the classes Monolithic and fullMonolithic
45  * - boundary conditions:
46  * -# Neumann
47  * -# Dirichlet
48  * -# Robin
49  * -# Fluxes (defective)
50  * -# absorbing \cite BadiaNobileVergara2008 :
51  * through the class flowConditions.
52  * - optional: computation of wall shear stress (not properly tested in parallel)
53  * - optional: computation of the largest singular values of the preconditioned matrix
54  *
55  * \b Features:
56  * This test by default solves the FSI probem discretized in time using the GCE or CE methods, implemented respectively
57  * in the files monolithicGE.hpp and monolithicGI.hpp . The geometry is that of a tube (benchmark test introduced in \cite Gerbeau2003).
58  * In this test the boundary conditions assigned are of type:
59  * - flux (defective b.c.) at the inlet
60  * - absorbing (see \cite BadiaNobileVergara2008) at the outlet
61  * - Robin b.c. on the solid external wall
62  * - Dirichlet homogeneous at the solid rings on the inlet-outlet (clamped tube).
63  *
64  * The output is written at every timestep, in HDF5 (if available) or ensight format.
65  * This test implements an inlet flux bundary condition for the first three time steps, then at the fourth time step
66  * the inlet boundary condition is replaced by a Neumann one (this mechanism is useful to implement rudimental valves).
67  * The outflow boundary condition is of absorbing type. At the outer wall for the structure a Robin condition is imposed.
68  */
69 
70 
71 #include <cassert>
72 #include <cstdlib>
73 
74 #include <boost/timer.hpp>
75 
76 #include <Epetra_ConfigDefs.h>
77 #ifdef EPETRA_MPI
78 #include <mpi.h>
79 #include <Epetra_MpiComm.h>
80 #else
81 #include <Epetra_SerialComm.h>
82 #endif
83 
84 
85 // LifeV includes
86 #include <lifev/core/fem/FESpace.hpp>
87 #include <lifev/core/array/MapEpetra.hpp>
88 #include <lifev/core/fem/BCHandler.hpp>
89 #include <lifev/core/LifeV.hpp>
90 
91 #include <lifev/core/algorithm/PreconditionerIfpack.hpp>
92 #include <lifev/core/algorithm/PreconditionerML.hpp>
93 #ifdef HAVE_NS_PREC
94 #include <lifev/core/algorithm/PreconditionerPCD.hpp>
95 #endif
96 
97 #include <lifev/fsi/solver/FSISolver.hpp>
98 #include <lifev/structure/solver/StructuralOperator.hpp>
99 #include <lifev/fsi/solver/FSIMonolithicGI.hpp>
100 
101 #include <lifev/core/filter/ExporterEnsight.hpp>
102 #include <lifev/core/filter/ExporterEmpty.hpp>
103 #ifdef HAVE_HDF5
104 #include <lifev/core/filter/ExporterHDF5.hpp>
105 #endif
106 
108 
109 class Problem
110 {
111 public:
112 
113  typedef std::shared_ptr<LifeV::FSISolver> FSISolverPtr_Type;
114 
115  typedef LifeV::FSIOperator::data_Type data_Type;
116  typedef LifeV::FSIOperator::dataPtr_Type dataPtr_Type;
117 
118  typedef LifeV::FSIOperator::vector_Type vector_Type;
119  typedef LifeV::FSIOperator::vectorPtr_Type vectorPtr_Type;
120 
121  typedef std::shared_ptr< LifeV::Exporter<LifeV::RegionMesh<LifeV::LinearTetra> > > filterPtr_Type;
122 
123  typedef LifeV::ExporterEnsight<LifeV::FSIOperator::mesh_Type> ensightFilter_Type;
124  typedef std::shared_ptr<ensightFilter_Type> ensightFilterPtr_Type;
125 #ifdef HAVE_HDF5
126  typedef LifeV::ExporterHDF5<LifeV::FSIOperator::mesh_Type> hdf5Filter_Type;
127  typedef std::shared_ptr<hdf5Filter_Type> hdf5FilterPtr_Type;
128 #endif
129  typedef LifeV::FactorySingleton<LifeV::Factory<LifeV::FSIOperator, std::string> > FSIFactory_Type;
130  typedef LifeV::RegionMesh<LifeV::LinearTetra> mesh_Type;
131  typedef LifeV::MapEpetra map_Type;
132  typedef LifeV::FESpace<mesh_Type, map_Type> fespace_Type;
133  /*!
134  This routine sets up the problem:
135 
136  -# create the standard boundary conditions for the fluid and
137  structure problems.
138 
139  -# initialize and setup the FSIsolver
140  */
141 
142  Problem ( GetPot const& data_file ) :
143  M_Tstart (0.),
144  M_saveEvery (1)
145  {
146  using namespace LifeV;
147 
148  FSIOperator::solid_Type::material_Type::isotropicLaw_Type::StructureIsotropicMaterialFactory::instance().registerProduct ( "linearVenantKirchhoff", &FSIOperator::createVenantKirchhoffLinear );
149  FSIOperator::solid_Type::material_Type::isotropicLaw_Type::StructureIsotropicMaterialFactory::instance().registerProduct ( "exponential", &FSIOperator::createExponentialMaterialNonLinear );
150  FSIOperator::solid_Type::material_Type::isotropicLaw_Type::StructureIsotropicMaterialFactory::instance().registerProduct ( "neoHookean", &FSIOperator::createNeoHookeanMaterialNonLinear );
151  FSIOperator::solid_Type::material_Type::isotropicLaw_Type::StructureIsotropicMaterialFactory::instance().registerProduct ( "nonLinearVenantKirchhoff", &FSIOperator::createVenantKirchhoffNonLinear );
152 
153  std::cout << "register MonolithicGE : " << FSIMonolithicGE::S_register << std::endl;
154  std::cout << "register MonolithicGI : " << FSIMonolithicGI::S_register << std::endl;
155 #ifdef HAVE_NS_PREC
156  std::cout << "register PCD : " << PreconditionerPCD::S_register << std::endl;
157 #endif
158 
159  //bool reg=FSIMonolithicGI::S_register&&FSIMonolithicGE::S_register;
160 
161  M_data = dataPtr_Type ( new data_Type() );
162  M_data->setup ( data_file );
163 
164  M_fsi = FSISolverPtr_Type ( new FSISolver( ) );
165  MPI_Barrier ( MPI_COMM_WORLD );
166 
167  M_fsi->setData ( M_data );
168  M_fsi->FSIOper()->setDataFile ( data_file ); //TO BE REMOVED!
169 
170  // Setting FESpace and DOF
171 
172  std::string fluidMeshPartitioned = data_file ( "problem/fluidMeshPartitioned", "none" );
173  std::string solidMeshPartitioned = data_file ( "problem/solidMeshPartitioned", "none" );
174 #ifdef HAVE_HDF5
175  if ( fluidMeshPartitioned.compare ( "none" ) )
176  {
177  FSIOperator::meshFilter_Type fluidMeshFilter ( data_file, fluidMeshPartitioned );
178  fluidMeshFilter.setComm ( M_fsi->FSIOper()->worldComm() );
179  FSIOperator::meshFilter_Type solidMeshFilter ( data_file, solidMeshPartitioned );
180  solidMeshFilter.setComm ( M_fsi->FSIOper( )->worldComm( ) );
181  M_fsi->FSIOper( )->partitionMeshes ( fluidMeshFilter, solidMeshFilter );
182  M_fsi->FSIOper( )->setupFEspace( );
183  M_fsi->FSIOper( )->setupDOF ( fluidMeshFilter );
184  fluidMeshFilter.closeFile( );
185  solidMeshFilter.closeFile( );
186  }
187  else
188 #endif
189  {
190  M_fsi->FSIOper( )->partitionMeshes( );
191  M_fsi->FSIOper( )->setupFEspace( );
192  M_fsi->FSIOper( )->setupDOF( );
193  }
194 
195  MPI_Barrier ( MPI_COMM_WORLD );
196 
197  M_fsi->setFluidBC ( BCh_monolithicFlux ( true ) );
198  M_fsi->setSolidBC ( BCh_monolithicSolid ( *M_fsi->FSIOper( ) ) );
199 
200  M_fsi->setup();
201 
202  M_fsi->setFluidBC ( BCh_monolithicFluid ( *M_fsi->FSIOper( ) ) );
203  M_fsi->setHarmonicExtensionBC ( BCh_harmonicExtension ( *M_fsi->FSIOper( ) ) );
204 
205  dynamic_cast<LifeV::FSIMonolithic*> (M_fsi->FSIOper().get() )->mergeBCHandlers();
206 
207  std::string const exporterType = data_file ( "exporter/type", "ensight" );
208  std::string const fluidName = data_file ( "exporter/fluid/filename", "fluid" );
209  std::string const solidName = data_file ( "exporter/solid/filename", "solid" );
210 
211 #ifdef HAVE_HDF5
212  if (exporterType.compare ("hdf5") == 0)
213  {
214  M_exporterFluid.reset ( new hdf5Filter_Type ( data_file, fluidName) );
215  M_exporterSolid.reset ( new hdf5Filter_Type ( data_file, solidName) );
216  }
217  else
218 #endif
219  {
220  if (exporterType.compare ("none") == 0)
221  {
222  M_exporterFluid.reset ( new ExporterEmpty<RegionMesh<LinearTetra> > ( data_file, M_fsi->FSIOper()->uFESpace().mesh(), fluidName, M_fsi->FSIOper()->uFESpace().map().comm().MyPID() ) );
223  M_exporterSolid.reset ( new ExporterEmpty<RegionMesh<LinearTetra> > ( data_file, M_fsi->FSIOper()->dFESpace().mesh(), solidName, M_fsi->FSIOper()->uFESpace().map().comm().MyPID() ) );
224  }
225  else
226  {
227  M_exporterFluid.reset ( new ensightFilter_Type ( data_file, fluidName) );
228  M_exporterSolid.reset ( new ensightFilter_Type ( data_file, solidName) );
229  }
230  }
231 
232 
233  // load using ensight/hdf5
234  M_saveEvery = data_file ("exporter/saveEvery", 1);
235 
236  // load using ensight/hdf5
237  std::string restartType (data_file ("importer/restartType", "newSimulation") );
238  std::cout << "The load state is: " << restartType << std::endl;
239 
240  if (!restartType.compare ("Stokes") )
241  {
242  initializeStokes (data_file);
243  }
244  else if (!restartType.compare ("restartFSI") )
245  {
246  restartFSI (data_file);
247  }
248  else
249  {
250  M_fsi->initialize();
251 
252  M_velAndPressure.reset ( new vector_Type ( M_fsi->FSIOper()->fluid().getMap(), M_exporterFluid->mapType() ) );
253 
254  M_fluidDisp.reset ( new vector_Type ( M_fsi->FSIOper()->mmFESpace().map(), M_exporterFluid->mapType() ) );
255 
256  M_solidDisp.reset ( new vector_Type ( M_fsi->FSIOper()->dFESpace().map(), M_exporterSolid->mapType() ) );
257  }
258 
259 
260  M_exporterFluid->setMeshProcId (M_fsi->FSIOper()->uFESpace().mesh(), M_fsi->FSIOper()->uFESpace().map().comm().MyPID() );
261  M_exporterSolid->setMeshProcId (M_fsi->FSIOper()->dFESpace().mesh(), M_fsi->FSIOper()->dFESpace().map().comm().MyPID() );
262  M_exporterFluid->addVariable ( ExporterData<FSIOperator::mesh_Type>::VectorField, "f-velocity",
263  M_fsi->FSIOper()->uFESpacePtr(), M_velAndPressure, UInt (0) );
264  M_exporterFluid->addVariable ( ExporterData<FSIOperator::mesh_Type>::ScalarField, "f-pressure",
265  M_fsi->FSIOper()->pFESpacePtr(), M_velAndPressure,
266  UInt (3 * M_fsi->FSIOper()->uFESpace().dof().numTotalDof() ) );
267 
268  M_exporterFluid->addVariable ( ExporterData<FSIOperator::mesh_Type>::VectorField, "f-displacement",
269  M_fsi->FSIOper()->mmFESpacePtr(), M_fluidDisp, UInt (0) );
270 
271 
272 
273  // M_solidVel.reset ( new vector_Type( M_fsi->FSIOper()->dFESpace().map(), M_exporterSolid->mapType() ));
274  // M_WS.reset ( new vector_Type( M_fsi->FSIOper()->dFESpace().map(), M_exporterSolid->mapType() ));
275 
276  M_exporterSolid->addVariable ( ExporterData<FSIOperator::mesh_Type>::VectorField, "s-displacement",
277  M_fsi->FSIOper()->dFESpacePtr(), M_solidDisp, UInt (0) );
278  // M_exporterSolid->addVariable( ExporterData<FSIOperator::mesh_Type>::VectorField, "s-velocity",
279  // M_fsi->FSIOper()->dFESpacePtr(), M_solidVel, UInt(0) );
280  // M_exporterSolid->addVariable( ExporterData<FSIOperator::mesh_Type>::VectorField, "s-ws",
281  // M_fsi->FSIOper()->dFESpacePtr(), M_WS, UInt(0) );
282 
283 
284 
285  M_data->dataFluid()->dataTime()->setInitialTime ( M_Tstart );
286  M_data->dataFluid()->dataTime()->setTime ( M_data->dataFluid()->dataTime()->initialTime() );
287  M_data->dataSolid()->dataTime()->setInitialTime ( M_Tstart );
288  M_data->dataSolid()->dataTime()->setTime ( M_data->dataFluid()->dataTime()->initialTime() );
289  M_data->timeDataALE()->setInitialTime ( M_Tstart );
290  M_data->timeDataALE()->setTime ( M_data->dataFluid()->dataTime()->initialTime() );
291  }
292 
293  /*!
294  This routine runs the temporal loop
295  */
296  void
297  run()
298  {
299  boost::timer _overall_timer;
300 
301  LifeV::UInt iter = 1;
302  //LifeV::UInt offset=dynamic_cast<LifeV::FSIMonolithic*>(M_fsi->FSIOper().get())->offset();
303 
304  dynamic_cast<LifeV::FSIMonolithic*> (M_fsi->FSIOper().get() )->enableStressComputation (1);
305 
306  // #ifdef HAVE_HDF5
307  // if (M_exporterFluid->mapType() == LifeV::Unique)
308  // {
309  // M_exporterFluid->postProcess( M_Tstart-M_data->dataFluid()->dataTime()->getTimeStep() );//ugly way to avoid that hdf5 starts with a deformed mesh
310  // M_exporterSolid->postProcess( M_Tstart-M_data->dataSolid()->dataTime()->getTimeStep() );//ugly way to avoid that hdf5 starts with a deformed mesh
311  // }
312  // #endif
313 
314  for ( ; M_data->dataFluid()->dataTime()->canAdvance(); M_data->dataFluid()->dataTime()->updateTime(), M_data->dataSolid()->dataTime()->updateTime(), ++iter)
315  {
316  boost::timer _timer;
317 
318  if (iter % M_saveEvery == 0)
319  {
320  M_fsi->FSIOper()->exportSolidDisplacement (*M_solidDisp); // displacement(), M_offset);
321  //M_fsi->FSIOper()->exportSolidVelocity(*M_solidVel);// displacement(), M_offset);
322 
323  M_fsi->FSIOper()->exportFluidVelocityAndPressure (*M_velAndPressure);
324  M_exporterSolid->postProcess ( M_data->dataFluid()->dataTime()->time() );
325 
326  //*M_WS= *(dynamic_cast<LifeV::FSIMonolithic*>(M_fsi->FSIOper().get())->/*WS());//*/computeStress());
327 
328  *M_fluidDisp = M_fsi->FSIOper()->meshDisp();
329  M_exporterFluid->postProcess ( M_data->dataFluid()->dataTime()->time() );
330  }
331 
332  M_fsi->iterate();
333 
334 
335  M_fsi->FSIOper()->displayer().leaderPrintMax ("[fsi_run] Iteration ", iter);
336  M_fsi->FSIOper()->displayer().leaderPrintMax (" was done in : ", _timer.elapsed() );
337 
338  // std::cout << "solution norm " << iter << " : "
339  // << M_fsi->displacement().norm2() << "\n";
340 
341  }
342  // if (M_data->method().compare("monolithicGI"))
343  // {
344  // M_fsi->FSIOper()->iterateMesh(M_fsi->displacement());
345 
346  // M_solidDisp->subset(M_fsi->displacement(), offset);
347 
348  // *M_velAndPressure = M_fsi->displacement();
349  // M_exporterSolid->postProcess( M_data->dataFluid()->dataTime()->time() );
350  // *M_fluidDisp = M_fsi->FSIOper()->meshMotion().disp();
351  // M_exporterFluid->postProcess( M_data->dataFluid()->dataTime()->time() );
352  // }
353 
354  std::cout << "Total computation time = "
355  << _overall_timer.elapsed() << "s" << "\n";
356 
357  }
358 
359 private:
360 
361  void initializeStokes ( GetPot const& data_file);
362  void restartFSI ( GetPot const& data_file);
363 
364  FSISolverPtr_Type M_fsi;
365  dataPtr_Type M_data;
366 
367  filterPtr_Type M_exporterSolid;
368  filterPtr_Type M_exporterFluid;
369  filterPtr_Type M_importerSolid;
370  filterPtr_Type M_importerFluid;
371  vectorPtr_Type M_velAndPressure;
372  vectorPtr_Type M_fluidDisp;
373  vectorPtr_Type M_solidDisp;
374 
375  std::vector<vectorPtr_Type> M_solidStencil;
376  std::vector<vectorPtr_Type> M_fluidStencil;
377  std::vector<vectorPtr_Type> M_ALEStencil;
378 
379  //vectorPtr_Type M_solidVel;
380  LifeV::Real M_Tstart;
381  // vectorPtr_Type M_WS;
382  LifeV::UInt M_saveEvery;
383 
384  struct RESULT_CHANGED_EXCEPTION
385  {
386  public:
387  RESULT_CHANGED_EXCEPTION (LifeV::Real time)
388  {
389  std::cout << "Some modifications led to changes in the l2 norm of the solution at time" << time << std::endl;
390  }
391  };
392 };
393 
394 struct FSIChecker
395 {
396  FSIChecker ( GetPot const& _data_file ) :
398  {}
399 
400  void operator() ()
401  {
402  std::shared_ptr<Problem> fsip;
403 
404  try
405  {
406  fsip = std::shared_ptr<Problem> ( new Problem ( data_file ) );
407 
408  fsip->run();
409  }
410  catch ( std::exception const& _ex )
411  {
412  std::cout << "caught exception : " << _ex.what() << "\n";
413  }
414 
415  }
416 
418  LifeV::Vector disp;
419 };
420 
421 
422 
423 namespace LifeV
424 {
425 namespace
426 {
427 static bool regIF = (PRECFactory::instance().registerProduct ( "Ifpack", &createIfpack ) );
428 static bool regML = (PRECFactory::instance().registerProduct ( "ML", &createML ) );
429 }
430 }
431 
432 
433 int main (int argc, char** argv)
434 {
435 #ifdef HAVE_MPI
436  MPI_Init (&argc, &argv);
437 #else
438  std::cout << "% using serial Version" << std::endl;
439 #endif
440 
441  GetPot command_line (argc, argv);
442  const bool check = command_line.search (2, "-c", "--check");
443 
444  if (check)
445  {
446  GetPot data_fileGCE ("dataGCE");
447  FSIChecker _GCE_check ( data_fileGCE );
448  _GCE_check();
449 
450  GetPot data_fileCE ("dataCE");
451  FSIChecker _CE_check (data_fileCE);
452  _CE_check();
453 
454 #ifdef HAVE_MPI
455  MPI_Finalize();
456 #endif
457  return 0;
458  }
459  else
460  {
461  const std::string data_file_name = command_line.follow ("data", 2, "-f", "--file");
462  GetPot data_file (data_file_name);
463  FSIChecker _sp_check ( data_file );
464  _sp_check();
465  }
466 #ifdef HAVE_MPI
467  MPI_Finalize();
468 #endif
469 
470  return 0;
471 
472 }
473 
474 void Problem::restartFSI ( GetPot const& data_file)
475 {
476 
477  using namespace LifeV;
478 
479  typedef FSIOperator::mesh_Type mesh_Type;
480 
481  //Creating the pointer to the filter
482  std::string const importerType = data_file ( "importer/type", "ensight");
483  std::string const fluidName = data_file ( "importer/fluid/filename", "fluid");
484  std::string const solidName = data_file ( "importer/solid/filename", "solid");
485 
486  std::string const loadInitSol = data_file ( "importer/initSol", "00000");
487  std::string const loadInitSolFD = data_file ("importer/initSolFD", "-1");
488  std::string iterationString;
489 
490  M_Tstart = data_file ( "fluid/time_discretization/initialtime", 0.);
491 
492  std::cout << "The file for fluid is : " << fluidName << std::endl;
493  std::cout << "The file for solid is : " << solidName << std::endl;
494  std::cout << "The importerType is : " << importerType << std::endl;
495  std::cout << "The iteration is : " << loadInitSol << std::endl;
496  std::cout << "For the fluid disp is : " << loadInitSolFD << std::endl;
497  std::cout << "Starting time : " << M_Tstart << std::endl;
498 
499 #ifdef HAVE_HDF5
500  if (importerType.compare ("hdf5") == 0)
501  {
502  M_importerFluid.reset ( new hdf5Filter_Type ( data_file, fluidName) );
503  M_importerSolid.reset ( new hdf5Filter_Type ( data_file, solidName) );
504  }
505  else
506 #endif
507  {
508  if (importerType.compare ("none") == 0)
509  {
510  M_importerFluid.reset ( new ExporterEmpty<mesh_Type > ( data_file, M_fsi->FSIOper()->uFESpace().mesh(), "fluid", M_fsi->FSIOper()->uFESpace().map().comm().MyPID() ) );
511  M_importerSolid.reset ( new ExporterEmpty<mesh_Type > ( data_file, M_fsi->FSIOper()->dFESpace().mesh(), "solid", M_fsi->FSIOper()->uFESpace().map().comm().MyPID() ) );
512  }
513  else
514  {
515  M_importerFluid.reset ( new ensightFilter_Type ( data_file, fluidName) );
516  M_importerSolid.reset ( new ensightFilter_Type ( data_file, solidName) );
517  }
518  }
519 
520  M_importerFluid->setMeshProcId (M_fsi->FSIOper()->uFESpace().mesh(), M_fsi->FSIOper()->uFESpace().map().comm().MyPID() );
521  M_importerSolid->setMeshProcId (M_fsi->FSIOper()->dFESpace().mesh(), M_fsi->FSIOper()->dFESpace().map().comm().MyPID() );
522 
523  //Each of the vectors of the stencils has the dimension of the big vector
524  //Performing a cycle of the size of the timeAdvance classes for each problem
525  //The size of TimeAdvanceClass depends on the order of the BDF (data file)
526 
527  //It should work just initializing the timeAdvance classes
528  //Three stencils are needed (Fluid-Structure-Geometric)
529 
530 
531  vectorPtr_Type fluidSol (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique) );
532  vectorPtr_Type initFluid (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
533  vectorPtr_Type HarmonicSol (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
534  vectorPtr_Type structureSol (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
535  vectorPtr_Type temporarySol (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
536 
537  vectorPtr_Type vel (new vector_Type (M_fsi->FSIOper()->uFESpace().map(), M_importerFluid->mapType() ) );
538  vectorPtr_Type pressure (new vector_Type (M_fsi->FSIOper()->pFESpace().map(), M_importerFluid->mapType() ) );
539  vectorPtr_Type solidDisp (new vector_Type (M_fsi->FSIOper()->dFESpace().map(), M_importerSolid->mapType() ) );
540  vectorPtr_Type fluidDisp (new vector_Type (M_fsi->FSIOper()->mmFESpace().map(), M_importerFluid->mapType() ) );
541 
542  vectorPtr_Type firstFluidDisp (new vector_Type (M_fsi->FSIOper()->mmFESpace().map(), M_importerFluid->mapType() ) );
543 
544  UInt offset = dynamic_cast<LifeV::FSIMonolithic*> (M_fsi->FSIOper().get() )->offset();
545 
546  iterationString = loadInitSol;
547  std::cout << "Fluid size TimeAdvance:" << M_fsi->FSIOper()->fluidTimeAdvance()->size() << std::endl;
548  for (UInt iterInit = 0; iterInit < M_fsi->FSIOper()->fluidTimeAdvance()->size(); iterInit++ )
549  {
550 
551  /*!
552  definition of the vector to fill with the initialization.
553  */
554  temporarySol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
555  fluidSol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique) );
556  initFluid.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
557 
558 
559  LifeV::ExporterData<mesh_Type> initSolFluidVel (LifeV::ExporterData<mesh_Type>::VectorField, std::string ("f-velocity." + iterationString), M_fsi->FSIOper()->uFESpacePtr(), vel, UInt (0), LifeV::ExporterData<mesh_Type>::UnsteadyRegime );
560  LifeV::ExporterData<mesh_Type> initSolFluidPress (LifeV::ExporterData<mesh_Type>::ScalarField, std::string ("f-pressure." + iterationString), M_fsi->FSIOper()->pFESpacePtr(), pressure, UInt (0), LifeV::ExporterData<mesh_Type>::UnsteadyRegime );
561 
562  /*!load of the solutions*/
563  M_importerFluid->readVariable (initSolFluidVel);
564  M_importerFluid->readVariable (initSolFluidPress);
565 
566  fluidSol.reset ( new vector_Type (*pressure, Unique, Zero) );
567  initFluid->subset (*fluidSol, fluidSol->map(), UInt (0), (UInt) 3 * M_fsi->FSIOper()->uFESpace().dof().numTotalDof() );
568  vector_Type tmpVec (*initFluid);
569  //temporarySol->subset(*fluidSol, fluidSol->map(), UInt(0), (UInt)3*M_fsi->FSIOper()->uFESpace().dof().numTotalDof());
570  fluidSol.reset ( new vector_Type (*vel, Unique, Zero) );
571  tmpVec = *fluidSol;
572  *initFluid += tmpVec;
573 
574  // std::string firstFl="firstFluid";
575  // initFluid->spy(firstFl);
576 
577  std::cout << "Norm of first Fluid sol: " << initFluid->norm2() << std::endl;
578  *temporarySol = *initFluid;
579 
580  // if(!iterInit)
581  // {
582  // *un += *temporarySol;
583  // }
584 
585  M_fluidStencil.push_back (temporarySol);
586  //Updating string name
587  int iterations = std::atoi (iterationString.c_str() );
588  iterations--;
589 
590  std::ostringstream iter;
591  iter.fill ( '0' );
592  iter << std::setw (5) << ( iterations );
593  iterationString = iter.str();
594 
595  }
596 
597  iterationString = loadInitSol;
598  for (UInt iterInit = 0; iterInit < M_fsi->FSIOper()->solidTimeAdvance()->size(); iterInit++ )
599  {
600  /*!
601  definition of the vector to fill with the initialization.
602  */
603  temporarySol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
604  structureSol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
605 
606 
607  /*!Definition of the ExporterData, used to load the solution inside the previously defined vectors*/
608  LifeV::ExporterData<mesh_Type> initSolSolidDisp (LifeV::ExporterData<mesh_Type>::VectorField, "s-displacement." + iterationString, M_fsi->FSIOper()->dFESpacePtr(), solidDisp, UInt (0) /*offset*/, LifeV::ExporterData<mesh_Type>::UnsteadyRegime );
609 
610  /*!load of the solutions*/
611  M_importerSolid->readVariable (initSolSolidDisp);
612 
613  structureSol->subset (*solidDisp, solidDisp->map(), (UInt) 0, offset);
614 
615  *temporarySol = *structureSol / (M_fsi->FSIOper()->solid().rescaleFactor() );
616 
617  M_solidStencil.push_back (temporarySol);
618 
619  //Updating the string name for the next iteration
620  UInt iterations = std::atoi (iterationString.c_str() );
621  iterations--;
622 
623  std::ostringstream iter;
624  iter.fill ( '0' );
625  iter << std::setw (5) << ( iterations );
626  iterationString = iter.str();
627 
628  }
629 
630  Int convectiveTerm = 0; //loadInitSol;
631  if (!M_data->dataFluid()->domainVelImplicit() )
632  {
633  convectiveTerm = 1;
634  }
635  HarmonicSol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
636 
637  iterationString = loadInitSolFD;
638 
639  for (UInt iterInit = 0; iterInit < M_fsi->FSIOper()->ALETimeAdvance()->size() + convectiveTerm + 1; iterInit++ )
640  {
641  temporarySol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV:: Unique, Zero) );
642  /*!
643  definition of the vector to fill with the initialization.
644  */
645 
646  LifeV::ExporterData<mesh_Type> initSolFluidDisp (LifeV::ExporterData<mesh_Type>::VectorField, "f-displacement." + iterationString, M_fsi->FSIOper()->mmFESpacePtr(), fluidDisp, UInt (0), LifeV::ExporterData<mesh_Type>::UnsteadyRegime );
647 
648  /*!load of the solutions*/
649  M_importerFluid->readVariable (initSolFluidDisp);
650 
651  std::cout << "Reloaded Harmonic sol norm: " << fluidDisp->norm2() << std::endl;
652  if (iterInit == 0)
653  {
654  HarmonicSol->subset (*fluidDisp, fluidDisp->map(), (UInt) 0, dynamic_cast<LifeV::FSIMonolithicGI*> (M_fsi->FSIOper().get() )->mapWithoutMesh().map (Unique)->NumGlobalElements() );
655  if (!convectiveTerm)
656  {
657  *firstFluidDisp = *fluidDisp;
658  }
659  }
660  else if (convectiveTerm && iterInit == 1)
661  {
662  *firstFluidDisp = *fluidDisp;
663  }
664  else
665  {
666  M_ALEStencil.push_back (/*temporarySol*/fluidDisp);
667  }
668  //Updating the iteration String name
669  int iterations = std::atoi (iterationString.c_str() );
670  iterations--;
671 
672  std::ostringstream iter;
673  iter.fill ( '0' );
674  iter << std::setw (5) << ( iterations );
675  iterationString = iter.str();
676  }
677 
678  *M_fluidStencil[0] += *M_solidStencil[0];
679  *M_fluidStencil[0] += *HarmonicSol;
680  //this is going to be the global solutions returned by the method solution()
681 
682  M_fsi->initialize (M_fluidStencil, M_solidStencil, M_ALEStencil);
683 
684  if (!M_data->dataFluid()->domainVelImplicit() )
685  {
686  //The following is needed because (and if) of the extrapolation of the fluid domain velocity is used, i.e. M_domainVelImplicit
687  M_fsi->FSIOper()->ALETimeAdvance()->updateRHSFirstDerivative ( M_data->dataSolid()->dataTime()->timeStep() );
688  M_fsi->FSIOper()->ALETimeAdvance()->shiftRight (*firstFluidDisp);
689  }
690 
691  M_velAndPressure.reset ( new vector_Type ( M_fsi->FSIOper()->fluid().getMap(), M_importerFluid->mapType() ) );
692  M_velAndPressure->subset (*pressure, pressure->map(), UInt (0), (UInt) 3 * M_fsi->FSIOper()->uFESpace().dof().numTotalDof() );
693  *M_velAndPressure += *vel;
694 
695  M_fluidDisp.reset ( new vector_Type ( *fluidDisp, M_importerFluid->mapType() ) );
696 
697  M_solidDisp.reset ( new vector_Type ( *solidDisp, M_importerSolid->mapType() ) );
698 
699  M_data->dataFluid()->dataTime()->updateTime(), M_data->dataSolid()->dataTime()->updateTime();
700 }
701 
702 
703 void Problem::initializeStokes ( GetPot const& data_file)
704 {
705 
706  using namespace LifeV;
707 
708  typedef FSIOperator::mesh_Type mesh_Type;
709 
710 
711  //Creating the pointer to the filter
712  filterPtr_Type importer;
713  std::string const importerType = data_file ( "importer/type", "hdf5");
714  std::string const filename = data_file ( "importer/fluid/filename", "fluid");
715 
716  std::string const loadInitSol = data_file ( "importer/initSol", "00000");
717  std::string const loadInitSolFD = data_file ("importer/initSolFD", "-1");
718  std::string iterationString;
719 
720  std::cout << "The filename is : " << filename << std::endl;
721  std::cout << "The importerType is: " << importerType << std::endl;
722 
723 #ifdef HAVE_HDF5
724  if (importerType.compare ("hdf5") == 0)
725  {
726  importer.reset ( new hdf5Filter_Type ( data_file, filename) );
727  }
728  else
729 #endif
730  {
731  if (importerType.compare ("none") == 0)
732  {
733  importer.reset ( new ExporterEmpty<RegionMesh<LinearTetra> > ( data_file, M_fsi->FSIOper()->uFESpace().mesh(), "fluid", M_fsi->FSIOper()->uFESpace().map().comm().MyPID() ) );
734  }
735  else
736  {
737  importer.reset ( new ensightFilter_Type ( data_file, filename) );
738  }
739  }
740 
741  importer->setMeshProcId (M_fsi->FSIOper()->uFESpace().mesh(), M_fsi->FSIOper()->uFESpace().map().comm().MyPID() );
742 
743 
744  std::vector<vectorPtr_Type> fluidStencil;;
745 
746  vectorPtr_Type temporarySol (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique) );
747 
748  //UInt offset=dynamic_cast<LifeV::FSIMonolithic*>(M_fsi->FSIOper().get())->offset();
749 
750  iterationString = loadInitSol;
751  for (UInt iterInit = 0; iterInit < M_fsi->FSIOper()->fluidTimeAdvance()->size(); ++iterInit )
752  {
753 
754  /*!
755  definition of the vector to fill with the initialization.
756  */
757  temporarySol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), Unique) );
758 
759  vectorPtr_Type vel (new vector_Type (M_fsi->FSIOper()->uFESpace().map(), M_importerFluid->mapType() ) );
760  vectorPtr_Type pressure (new vector_Type (M_fsi->FSIOper()->pFESpace().map(), M_importerFluid->mapType() ) );
761 
762  *vel *= 0.0;
763  *pressure *= 0.0;
764 
765  LifeV::ExporterData<mesh_Type> initSolFluidVel (LifeV::ExporterData<mesh_Type>::VectorField, std::string ("f-velocity." + iterationString), M_fsi->FSIOper()->uFESpacePtr(), vel, UInt (0), LifeV::ExporterData<mesh_Type>::UnsteadyRegime );
766  LifeV::ExporterData<mesh_Type> initSolFluidPress (LifeV::ExporterData<mesh_Type>::ScalarField, std::string ("f-pressure." + iterationString), M_fsi->FSIOper()->pFESpacePtr(), pressure, UInt (0), LifeV::ExporterData<mesh_Type>::UnsteadyRegime );
767 
768  /*!load of the solutions*/
769  importer->readVariable (initSolFluidVel);
770  importer->readVariable (initSolFluidPress);
771 
772  int iterations = std::atoi (iterationString.c_str() );
773  iterations--;
774 
775  std::ostringstream iter;
776  iter.fill ( '0' );
777  iter << std::setw (5) << ( iterations );
778  iterationString = iter.str();
779 
780  *temporarySol = *vel + *pressure;
781  fluidStencil.push_back (temporarySol);
782  }
783 
784  // M_fsi->initialize(fluidStencil);
785 
786 }
int main(int argc, char **argv)
GetPot(const int argc_, char **argv_, const char *FieldSeparator=0x0)
Definition: GetPot.hpp:507
FSIChecker(GetPot const &_data_file)
GetPot(const char *FileName, const char *CommentStart=0x0, const char *CommentEnd=0x0, const char *FieldSeparator=0x0)
Definition: GetPot.hpp:541
bool search(unsigned No, const char *P,...)
Definition: GetPot.hpp:1217
GetPot(const STRING_VECTOR &FileNameList)
Definition: GetPot.hpp:645
const std::string follow(const char *Default, unsigned No, const char *Option,...)
Definition: GetPot.hpp:1510