LifeV
NavierStokesSolverBlocks.cpp
Go to the documentation of this file.
1 #include <lifev/navier_stokes_blocks/solver/NavierStokesSolverBlocks.hpp>
2 
3 
4 namespace LifeV
5 {
6 
7 class SignFunction
8 {
9 public:
10  typedef Real return_Type;
11 
12  inline return_Type operator() (const Real& a)
13  {
14  if ( a < 0 )
15  {
16  return a;
17  }
18  else
19  {
20  return 0.0;
21  }
22  }
23 
24  SignFunction( const std::shared_ptr<Epetra_Comm>& communicator ) { k = 0; z = 0; comm = communicator; }
25  SignFunction (const SignFunction&) {}
26  ~SignFunction() {}
27  int k, z;
28  std::shared_ptr<Epetra_Comm> comm;
29  int interrogate_k () { int k_global; comm->MaxAll(&k,&k_global,1); return k_global; }
30  int interrogate_z () { int z_global; comm->MaxAll(&z,&z_global,1); return z_global; }
31  void clear_k () { k = 0; }
32  void clear_z () { z = 0; }
33 };
34 
35 NavierStokesSolverBlocks::NavierStokesSolverBlocks(const dataFile_Type dataFile, const commPtr_Type& communicator):
36  M_comm(communicator),
37  M_dataFile(dataFile),
38  M_displayer(communicator),
39  M_graphIsBuilt(false),
40  M_oper(new Operators::NavierStokesOperator),
41  M_operLoads(new Operators::NavierStokesOperator),
42  M_invOper(),
43  M_fullyImplicit(false),
44  M_graphPCDisBuilt(false),
45  M_steady ( dataFile("fluid/miscellaneous/steady", false) ),
46  M_density ( dataFile("fluid/physics/density", 1.0 ) ),
47  M_viscosity ( dataFile("fluid/physics/viscosity", 0.035 ) ),
48  M_useStabilization ( dataFile("fluid/stabilization/use", false) ),
49  M_stabilizationType ( dataFile("fluid/stabilization/type", "none") ),
50  M_nonconforming ( dataFile("fluid/nonconforming", false) ),
51  M_imposeWeakBC ( dataFile("fluid/weak_bc/use", false) ),
52  M_flagWeakBC ( dataFile("fluid/weak_bc/flag", 31) ),
53  M_meshSizeWeakBC ( dataFile("fluid/weak_bc/mesh_size", 0.0) ),
54  M_computeAerodynamicLoads ( dataFile("fluid/forces/compute", false) ),
55  M_flagBody ( dataFile("fluid/forces/flag", 31 ) ),
56  M_solve_blocks ( dataFile("fluid/solve_blocks", true ) ),
57  M_useFastAssembly ( dataFile("fluid/use_fast_assembly", false ) ),
58  M_orderBDF ( dataFile("fluid/time_discretization/BDF_order", 2 ) ),
59  M_orderVel ( dataFile("fluid/stabilization/vel_order", 2 ) )
60 {
61  M_prec.reset ( Operators::NSPreconditionerFactory::instance().createObject (dataFile("fluid/preconditionerType","SIMPLE")));
62 }
63 
64 NavierStokesSolverBlocks::~NavierStokesSolverBlocks()
65 {
66 }
67 
68 void NavierStokesSolverBlocks::setParameters( )
69 {
70  std::string optionsPrec = M_dataFile("fluid/options_preconditioner","solverOptionsFast");
71  optionsPrec += ".xml";
72  Teuchos::RCP<Teuchos::ParameterList> solversOptions = Teuchos::getParametersFromXmlFile (optionsPrec);
73  M_prec->setOptions(*solversOptions);
74  setSolversOptions(*solversOptions);
75 }
76 
77 void NavierStokesSolverBlocks::setup(const meshPtr_Type& mesh, const int& id_domain)
78 {
79  std::string uOrder;
80  std::string pOrder;
81 
82  if ( !M_nonconforming )
83  {
84  uOrder = M_dataFile("fluid/space_discretization/vel_order","P1");
85  pOrder = M_dataFile("fluid/space_discretization/pres_order","P1");
86  M_stiffStrain = M_dataFile("fluid/space_discretization/stiff_strain", false );
87 
88  }
89  else
90  {
91  if ( id_domain == 0 )
92  {
93  uOrder = M_dataFile("fluid/master_space_discretization/vel_order","P1");
94  pOrder = M_dataFile("fluid/master_space_discretization/pres_order","P1");
95  M_stiffStrain = M_dataFile("fluid/master_space_discretization/stiff_strain", false );
96 
97  }
98  else if ( id_domain == 1 )
99  {
100  uOrder = M_dataFile("fluid/slave_space_discretization/vel_order","P1");
101  pOrder = M_dataFile("fluid/slave_space_discretization/pres_order","P1");
102  M_stiffStrain = M_dataFile("fluid/slave_space_discretization/stiff_strain", false );
103  }
104  }
105 
106  // Penalization reverse flow
107  M_penalizeReverseFlow = M_dataFile("fluid/penalize_reverse_flow", false);
108  M_flagPenalizeReverseFlow = M_dataFile("fluid/flag_outflow", 2);
109 
110  M_uOrder = uOrder;
111  M_pOrder = pOrder;
112 
113  M_useGraph = M_dataFile("fluid/use_graph", false);
114 
115  M_velocityFESpace.reset (new FESpace<mesh_Type, map_Type> (mesh, uOrder, 3, M_comm) );
116  M_pressureFESpace.reset (new FESpace<mesh_Type, map_Type> (mesh, pOrder, 1, M_comm) );
117  M_velocityFESpaceScalar.reset (new FESpace<mesh_Type, map_Type> (mesh, uOrder, 1, M_comm) );
118 
119  M_fespaceUETA.reset( new ETFESpace_velocity(M_velocityFESpace->mesh(), &(M_velocityFESpace->refFE()), M_comm));
120  M_fespaceUETA_scalar.reset( new ETFESpace_pressure(M_velocityFESpaceScalar->mesh(), &(M_velocityFESpaceScalar->refFE()), M_comm));
121  M_fespacePETA.reset( new ETFESpace_pressure(M_pressureFESpace->mesh(), &(M_pressureFESpace->refFE()), M_comm));
122 
123  M_uExtrapolated.reset( new vector_Type ( M_velocityFESpace->map(), Repeated ) );
124  M_uExtrapolated->zero();
125 
126  M_velocity.reset( new vector_Type(M_velocityFESpace->map()) );
127  M_pressure.reset( new vector_Type(M_pressureFESpace->map()) );
128 
129  M_block00.reset( new matrix_Type(M_velocityFESpace->map() ) );
130  M_block01.reset( new matrix_Type(M_velocityFESpace->map() ) );
131  M_block10.reset( new matrix_Type(M_pressureFESpace->map() ) );
132  M_block11.reset( new matrix_Type(M_pressureFESpace->map() ) );
133 
134  M_block00->zero();
135  M_block01->zero();
136  M_block10->zero();
137  M_block11->zero();
138 
139  M_fullyImplicit = M_dataFile ( "newton/convectiveImplicit", false);
140 
141  M_monolithicMap.reset( new map_Type ( M_velocityFESpace->map() ) );
142  *M_monolithicMap += M_pressureFESpace->map();
143 
144  if ( M_useFastAssembly ) // currently works for fully implicit NS within FSI or for P2-P1
145  {
146  M_fastAssembler.reset( new FastAssemblerNS ( M_velocityFESpace->mesh(), M_comm,
147  &(M_velocityFESpace->refFE()), &(M_pressureFESpace->refFE()),
148  M_velocityFESpace, M_pressureFESpace,
149  &(M_velocityFESpace->qr()) ) );
150  M_fastAssembler->allocateSpace ( &(M_velocityFESpace->fe()), true );
151  }
152 
153 
154  if ( M_fullyImplicit )
155  {
156  M_displayer.leaderPrint ( " F - solving nonlinear Navier-Stokes\n");
157  M_relativeTolerance = M_dataFile ( "newton/abstol", 1.e-4);
158  M_absoluteTolerance = M_dataFile ( "newton/reltol", 1.e-4);
159  M_etaMax = M_dataFile ( "newton/etamax", 1e-4);
160  M_maxiterNonlinear = M_dataFile ( "newton/maxiter", 10);
161 
162  M_nonLinearLineSearch = M_dataFile ( "newton/NonLinearLineSearch", 0);
163 
164  if (M_comm->MyPID() == 0)
165  M_out_res.open ("residualsNewton");
166 
167  M_solution.reset( new vector_Type ( *M_monolithicMap ) );
168  M_solution->zero();
169  }
170 
171  if ( M_useStabilization )
172  {
173  if ( id_domain == 0 )
174  {
175  M_stabilization.reset ( StabilizationFactory::instance().createObject ( "SUPG_SEMI_IMPLICIT_ALE" ) );
176  }
177  else
178  {
179  M_stabilization.reset ( StabilizationFactory::instance().createObject ( M_dataFile("fluid/stabilization/type","none") ) );
180 
181  }
182 
183  if (M_comm->MyPID() == 0)
184  std::cout << "\nUsing the " << M_stabilization->label() << " stabilization\n\n";
185 
186  M_stabilization->setVelocitySpace(M_velocityFESpace);
187  M_stabilization->setPressureSpace(M_pressureFESpace);
188  M_stabilization->setDensity(M_density);
189  M_stabilization->setViscosity(M_viscosity);
190  int vel_order = M_dataFile ( "fluid/stabilization/vel_order", 1 );
191  M_stabilization->setConstant( vel_order );
192  M_stabilization->setBDForder ( M_dataFile ( "fluid/time_discretization/BDF_order", 1 ) );
193  M_stabilization->setCommunicator ( M_velocityFESpace->map().commPtr() );
194  M_stabilization->setTimeStep ( M_dataFile ( "fluid/time_discretization/timestep", 0.001 ) );
195  M_stabilization->setETvelocitySpace ( M_fespaceUETA );
196  M_stabilization->setETpressureSpace ( M_fespacePETA );
197  M_stabilization->setUseGraph( M_useGraph );
198  M_stabilization->setUseODEfineScale( M_dataFile("fluid/stabilization/ode_fine_scale", false ) );
199  }
200 
201  M_displayer.leaderPrintMax ( " Number of DOFs for the velocity = ", M_velocityFESpace->dof().numTotalDof()*3 ) ;
202  M_displayer.leaderPrintMax ( " Number of DOFs for the pressure = ", M_pressureFESpace->dof().numTotalDof() ) ;
203 
204 }
205 
206 void NavierStokesSolverBlocks::setExportFineScaleVelocity( ExporterHDF5<mesh_Type>& exporter, const int& numElementsTotal)
207 {
208  M_stabilization->setExportFineScaleVelocity ( exporter, numElementsTotal );
209 }
210 
211 void NavierStokesSolverBlocks::setSolversOptions(const Teuchos::ParameterList& solversOptions)
212 {
213  std::shared_ptr<Teuchos::ParameterList> monolithicOptions;
214  monolithicOptions.reset(new Teuchos::ParameterList(solversOptions.sublist("MonolithicOperator")) );
215  M_pListLinSolver = monolithicOptions;
216 }
217 
218 void NavierStokesSolverBlocks::buildGraphs()
219 {
220  M_displayer.leaderPrint ( " F - Pre-building the graphs... ");
221  LifeChrono chrono;
222  chrono.start();
223 
224  {
225  using namespace ExpressionAssembly;
226 
227  if ( !M_steady )
228  {
229  // Graph velocity mass -> block (0,0)
230  M_Mu_graph.reset (new Epetra_FECrsGraph (Copy, * (M_velocityFESpace->map().map (Unique) ), 0) );
231  buildGraph ( elements (M_fespaceUETA->mesh() ),
232  quadRuleTetra4pt,
233  M_fespaceUETA,
234  M_fespaceUETA,
235  M_density * dot ( phi_i, phi_j )
236  ) >> M_Mu_graph;
237  M_Mu_graph->GlobalAssemble();
238  M_Mu_graph->OptimizeStorage();
239  }
240 
241  // Graph block (0,1) of NS
242  M_Btranspose_graph.reset (new Epetra_FECrsGraph (Copy, * (M_velocityFESpace->map().map (Unique) ), 0) );
243  buildGraph ( elements (M_fespaceUETA->mesh() ),
244  quadRuleTetra4pt,
245  M_fespaceUETA,
246  M_fespacePETA,
247  value(-1.0) * phi_j * div(phi_i)
248  ) >> M_Btranspose_graph;
249  M_Btranspose_graph->GlobalAssemble( *(M_pressureFESpace->map().map (Unique)), *(M_velocityFESpace->map().map (Unique)) );
250  M_Btranspose_graph->OptimizeStorage();
251 
252  // Graph block (1,0) of NS
253  M_B_graph.reset (new Epetra_FECrsGraph (Copy, *(M_pressureFESpace->map().map (Unique)), 0) );
254  buildGraph ( elements (M_fespaceUETA->mesh() ),
255  quadRuleTetra4pt,
256  M_fespacePETA,
257  M_fespaceUETA,
258  phi_i * div(phi_j)
259  ) >> M_B_graph;
260  M_B_graph->GlobalAssemble( *(M_velocityFESpace->map().map (Unique)), *(M_pressureFESpace->map().map (Unique)) );
261  M_B_graph->OptimizeStorage();
262 
263  // Graph convective term, block (0,0)
264  M_C_graph.reset (new Epetra_FECrsGraph (Copy, * (M_velocityFESpace->map().map (Unique) ), 0) );
265  buildGraph ( elements (M_fespaceUETA->mesh() ),
266  quadRuleTetra4pt,
267  M_fespaceUETA,
268  M_fespaceUETA,
269  dot( M_density*value(M_fespaceUETA, *M_uExtrapolated)*grad(phi_j), phi_i)
270  ) >> M_C_graph;
271  M_C_graph->GlobalAssemble();
272  M_C_graph->OptimizeStorage();
273 
274  // Graph stiffness, block (0,0)
275  M_A_graph.reset (new Epetra_FECrsGraph (Copy, * (M_velocityFESpace->map().map (Unique) ), 0) );
276 
277  buildGraph ( elements (M_fespaceUETA->mesh() ),
278  quadRuleTetra4pt,
279  M_fespaceUETA,
280  M_fespaceUETA,
281  value( 0.5 * M_viscosity ) * dot( grad(phi_i) , grad(phi_j) + transpose(grad(phi_j)) )
282  ) >> M_A_graph;
283  M_A_graph->GlobalAssemble();
284  M_A_graph->OptimizeStorage();
285 
286  // Graph of entire block (0,0)
287  M_F_graph.reset (new Epetra_FECrsGraph (Copy, * (M_velocityFESpace->map().map (Unique) ), 0) );
288  if (M_stiffStrain)
289  {
290  buildGraph ( elements (M_fespaceUETA->mesh() ),
291  quadRuleTetra4pt,
292  M_fespaceUETA,
293  M_fespaceUETA,
294  dot ( phi_i, phi_j ) + // mass
295  dot( M_density*value(M_fespaceUETA, *M_uExtrapolated)*grad(phi_j), phi_i) + // convective term
296  value( 0.5 * M_viscosity ) * dot( grad(phi_i) + transpose(grad(phi_i)) , grad(phi_j) + transpose(grad(phi_j)) ) // stiffness
297  ) >> M_F_graph;
298 
299 
300  M_Jacobian_graph.reset (new Epetra_FECrsGraph (Copy, * (M_velocityFESpace->map().map (Unique) ), 0) );
301  buildGraph ( elements (M_fespaceUETA->mesh() ),
302  quadRuleTetra4pt,
303  M_fespaceUETA,
304  M_fespaceUETA,
305  dot ( phi_i, phi_j ) + // mass
306  dot( M_density*value(M_fespaceUETA, *M_uExtrapolated)*grad(phi_j), phi_i) + // convective term
307  dot( M_density * phi_j * grad(M_fespaceUETA, *M_uExtrapolated), phi_i ) + // part of the Jacobian
308  value( 0.5 * M_viscosity ) * dot( grad(phi_i) + transpose(grad(phi_i)) , grad(phi_j) + transpose(grad(phi_j)) ) // stiffness
309  ) >> M_Jacobian_graph;
310  M_Jacobian_graph->GlobalAssemble();
311  M_Jacobian_graph->OptimizeStorage();
312 
313  }
314  else
315  {
316  buildGraph ( elements (M_fespaceUETA->mesh() ),
317  quadRuleTetra4pt,
318  M_fespaceUETA,
319  M_fespaceUETA,
320  dot ( phi_i, phi_j ) + // mass
321  dot( M_density*value(M_fespaceUETA, *M_uExtrapolated)*grad(phi_j), phi_i) + // convective term
322  M_viscosity * dot( grad(phi_i) , grad(phi_j) + transpose(grad(phi_j)) ) // stiffness
323  ) >> M_F_graph;
324 
325 
326  M_Jacobian_graph.reset (new Epetra_FECrsGraph (Copy, * (M_velocityFESpace->map().map (Unique) ), 0) );
327  buildGraph ( elements (M_fespaceUETA->mesh() ),
328  quadRuleTetra4pt,
329  M_fespaceUETA,
330  M_fespaceUETA,
331  dot ( phi_i, phi_j ) + // mass
332  dot( M_density*value(M_fespaceUETA, *M_uExtrapolated)*grad(phi_j), phi_i) + // convective term
333  dot( M_density * phi_j * grad(M_fespaceUETA, *M_uExtrapolated), phi_i ) + // part of the Jacobian
334  M_viscosity * dot( grad(phi_i) , grad(phi_j) + transpose(grad(phi_j)) ) // stiffness
335  ) >> M_Jacobian_graph;
336  M_Jacobian_graph->GlobalAssemble();
337  M_Jacobian_graph->OptimizeStorage();
338 
339  }
340  M_F_graph->GlobalAssemble();
341  M_F_graph->OptimizeStorage();
342 
343  if ( M_useStabilization )
344  {
345  M_stabilization->buildGraphs();
346  }
347 
348  }
349 
350  M_graphIsBuilt = true;
351 
352  chrono.stop();
353  M_displayer.leaderPrintMax ( " done in ", chrono.diff() ) ;
354 }
355 
356 void NavierStokesSolverBlocks::buildSystem()
357 {
358  if ( M_useFastAssembly )
359  {
360  M_displayer.leaderPrint ( " F - Using fast assembly\n");
361 
362  if ( M_orderVel == 1 )
363  {
364  M_fastAssembler->setConstants_NavierStokes(M_density, M_viscosity, M_timeStep, M_orderBDF, 30.0, M_alpha );
365  }
366  else if ( M_orderVel == 2 )
367  {
368  M_fastAssembler->setConstants_NavierStokes(M_density, M_viscosity, M_timeStep, M_orderBDF, 60.0, M_alpha );
369  }
370 
371  M_fastAssembler->updateGeoQuantities ( &(M_velocityFESpace->fe()) );
372  }
373 
374  if ( M_useGraph && !M_graphIsBuilt )
375  buildGraphs();
376 
377  M_displayer.leaderPrint ( " F - Assembling constant terms... ");
378  LifeChrono chrono;
379  chrono.start();
380 
381  if ( M_useGraph )
382  {
383  M_Mu.reset (new matrix_Type ( M_velocityFESpace->map(), *M_Mu_graph ) );
384  M_Mu->zero();
385 
386  M_Btranspose.reset (new matrix_Type ( M_velocityFESpace->map(), *M_Btranspose_graph ) );
387  M_Btranspose->zero();
388 
389  M_B.reset (new matrix_Type ( M_pressureFESpace->map(), *M_B_graph ) );
390  M_B->zero();
391 
392  M_A.reset (new matrix_Type ( M_velocityFESpace->map(), *M_A_graph ) );
393  M_A->zero();
394 
395  M_C.reset (new matrix_Type ( M_velocityFESpace->map(), *M_C_graph ) );
396  M_C->zero();
397 
398  M_F.reset (new matrix_Type ( M_velocityFESpace->map(), *M_F_graph ) );
399  M_F->zero();
400 
401  M_Jacobian.reset (new matrix_Type ( M_velocityFESpace->map(), *M_Jacobian_graph ) );
402  M_Jacobian->zero();
403  }
404  else
405  {
406  M_Mu.reset (new matrix_Type ( M_velocityFESpace->map() ) ); // mass velocity
407  M_Mu->zero();
408 
409  M_Btranspose.reset (new matrix_Type ( M_velocityFESpace->map() ) ); // grad term
410  M_Btranspose->zero();
411 
412  M_B.reset (new matrix_Type ( M_pressureFESpace->map() ) ); // div term
413  M_B->zero();
414 
415  M_A.reset (new matrix_Type ( M_velocityFESpace->map() ) ); // stiffness
416  M_A->zero();
417 
418  M_C.reset (new matrix_Type ( M_velocityFESpace->map() ) ); // convective
419  M_C->zero();
420 
421  M_F.reset (new matrix_Type ( M_velocityFESpace->map() ) ); // linearization convective
422  M_F->zero();
423 
424  M_Jacobian.reset (new matrix_Type ( M_velocityFESpace->map() ) ); // jacobian
425  M_Jacobian->zero();
426  }
427 
428  if ( M_useFastAssembly )
429  {
430  M_fastAssembler->assemble_constant_terms( M_Mu, M_A, M_Btranspose, M_B );
431  M_Mu->globalAssemble();
432  M_Btranspose->globalAssemble( M_pressureFESpace->mapPtr(), M_velocityFESpace->mapPtr() );
433  M_B->globalAssemble( M_velocityFESpace->mapPtr(), M_pressureFESpace->mapPtr());
434  M_A->globalAssemble();
435  }
436  else
437  {
438  using namespace ExpressionAssembly;
439 
440  if ( !M_steady )
441  {
442  // Velocity mass -> block (0,0)
443  integrate ( elements (M_fespaceUETA->mesh() ),
444  M_velocityFESpace->qr(),
445  M_fespaceUETA,
446  M_fespaceUETA,
447  value(M_density) * dot ( phi_i, phi_j )
448  ) >> M_Mu;
449  M_Mu->globalAssemble();
450  }
451 
452  integrate( elements(M_fespaceUETA->mesh()),
453  M_velocityFESpace->qr(),
454  M_fespaceUETA,
455  M_fespacePETA,
456  value(-1.0) * phi_j * div(phi_i)
457  ) >> M_Btranspose;
458 
459  M_Btranspose->globalAssemble( M_pressureFESpace->mapPtr(), M_velocityFESpace->mapPtr() );
460 
461  integrate( elements(M_fespaceUETA->mesh()),
462  M_pressureFESpace->qr(),
463  M_fespacePETA,
464  M_fespaceUETA,
465  phi_i * div(phi_j)
466  ) >> M_B;
467 
468  M_B->globalAssemble( M_velocityFESpace->mapPtr(), M_pressureFESpace->mapPtr());
469 
470  integrate( elements(M_fespaceUETA->mesh()),
471  M_velocityFESpace->qr(),
472  M_fespaceUETA,
473  M_fespaceUETA,
474  value( M_viscosity ) * dot( grad(phi_i), grad(phi_j) + transpose(grad(phi_j)) )
475  ) >> M_A;
476 
477  M_A->globalAssemble();
478  }
479 
480  M_block01->zero();
481  M_block10->zero();
482  *M_block01 += *M_Btranspose;
483  *M_block10 += *M_B;
484 
485  if ( !M_useStabilization )
486  {
487  M_block01->globalAssemble( M_pressureFESpace->mapPtr(), M_velocityFESpace->mapPtr() );
488  M_block10->globalAssemble( M_velocityFESpace->mapPtr(), M_pressureFESpace->mapPtr());
489  }
490 
491  chrono.stop();
492  M_displayer.leaderPrintMax ( " done in ", chrono.diff() ) ;
493 }
494 
495 void NavierStokesSolverBlocks::updateSystem( const vectorPtr_Type& u_star, const vectorPtr_Type& rhs_velocity )
496 {
497  // Note that u_star HAS to extrapolated from outside. Hence it works also for FSI in this manner.
498  M_velocityExtrapolated.reset( new vector_Type ( *u_star, Unique ) );
499  M_uExtrapolated.reset( new vector_Type ( *u_star, Repeated ) );
500 
501  // Update convective term
502  M_C->zero();
503  {
504  if ( M_useFastAssembly )
505  {
506  M_fastAssembler->assembleConvective( M_C, *M_uExtrapolated );
507  }
508  else
509  {
510  using namespace ExpressionAssembly;
511  integrate( elements(M_fespaceUETA->mesh()),
512  M_velocityFESpace->qr(),
513  M_fespaceUETA,
514  M_fespaceUETA,
515  dot( M_density * value(M_fespaceUETA, *M_uExtrapolated)*grad(phi_j), phi_i) // semi-implicit treatment of the convective term
516  )
517  >> M_C;
518  }
519  if ( M_imposeWeakBC )
520  {
521  // Reference for essential BCs imposed weakly:
522  // - Y. Bazilevs, K. Takizawa and E. Tezduyar. Computational Fluid-Structre Interaction. Page 70.
523  using namespace ExpressionAssembly;
524  QuadratureBoundary myBDQR (buildTetraBDQR (quadRuleTria7pt) );
525 
526  M_block00_weakBC.reset ( new matrix_Type ( M_velocityFESpace->map() ) );
527  M_block00_weakBC->zero();
528 
529  integrate( boundary(M_fespaceUETA->mesh(), M_flagWeakBC),
530  myBDQR,
531  M_fespaceUETA,
532  M_fespaceUETA,
533  value(-1.0)*value(M_viscosity)*dot( phi_i, ( grad(phi_j) + transpose( grad(phi_j) ) ) * Nface )
534  -value(M_viscosity)*dot( ( grad(phi_i) + transpose( grad(phi_i) ) ) * Nface, phi_j )
535  -value(M_density)*dot( phi_i, dot( value(M_fespaceUETA, *M_uExtrapolated), Nface ) * phi_j )
536  +value(4.0)*value(M_viscosity)/value(M_meshSizeWeakBC)*dot( phi_i - dot( phi_i, Nface)*Nface, phi_j - dot( phi_j, Nface)*Nface )
537  +value(4.0)*value(M_viscosity)/value(M_meshSizeWeakBC)*( dot( phi_i, Nface)*dot( phi_j, Nface) )
538  )
539  >> M_block00_weakBC;
540  M_block00_weakBC->globalAssemble();
541  *M_C += *M_block00_weakBC;
542  }
543 
544  if ( M_penalizeReverseFlow )
545  {
546  this->setupPostProc();
547  std::shared_ptr<SignFunction> signEvaluation(new SignFunction(M_comm));
548 
549  using namespace ExpressionAssembly;
550 
551  // Penalizing back flow
552  QuadratureBoundary myBDQR (buildTetraBDQR (quadRuleTria7pt) );
553  integrate (
554  boundary (M_fespaceUETA->mesh(), M_flagPenalizeReverseFlow),
555  myBDQR,
556  M_fespaceUETA,
557  M_fespaceUETA,
558  value(-1.0*M_density)*eval( signEvaluation, dot(value(M_fespaceUETA, *M_uExtrapolated),Nface))*dot(phi_i, phi_j)
559  )
560  >> M_C;
561  }
562  }
563  M_C->globalAssemble();
564 
565  // Get the matrix corresponding to the block (0,0)
566  M_block00->zero();
567  *M_block00 += *M_Mu;
568  *M_block00 *= M_alpha/M_timeStep;
569  *M_block00 += *M_A;
570  *M_block00 += *M_C;
571  if ( !M_useStabilization )
572  M_block00->globalAssemble();
573 
574  // Get the right hand side with inertia contribution
575  M_rhs.reset( new vector_Type ( M_velocityFESpace->map(), Unique ) );
576  M_rhs->zero();
577 
578  if ( !M_steady )
579  *M_rhs = *M_Mu* (*rhs_velocity);
580 
581  M_block01->zero();
582  *M_block01 += *M_Btranspose;
583 
584  M_block10->zero();
585  *M_block10 += *M_B;
586 
587  if ( !M_fullyImplicit && M_useStabilization )
588  {
589  M_displayer.leaderPrint ( "\tF - Assembling semi-implicit stabilization terms... ");
590  LifeChrono chrono;
591  chrono.start();
592 
593  if ( M_useStabilization )
594  {
595  if ( M_stabilizationType.compare("VMSLES_SEMI_IMPLICIT")==0 )
596  M_stabilization->apply_matrix( *u_star, *M_pressure_extrapolated, *rhs_velocity);
597  else
598  {
599  if ( M_useFastAssembly )
600  {
601  M_stabilization->setFastAssembler( M_fastAssembler );
602  M_fastAssembler->setAlpha(M_alpha);
603  M_fastAssembler->setTimeStep(M_timeStep);
604  }
605 
606  M_stabilization->apply_matrix( *u_star );
607  }
608  }
609 
610  *M_block00 += *M_stabilization->block_00();
611 
612  *M_block01 += *M_stabilization->block_01();
613 
614  *M_block10 += *M_stabilization->block_10();
615 
616  M_block11->zero();
617  *M_block11 += *M_stabilization->block_11();
618  M_block11->globalAssemble();
619 
620  M_rhs_pressure.reset( new vector_Type ( M_pressureFESpace->map(), Unique ) );
621  M_rhs_pressure->zero();
622 
623  M_stabilization->apply_vector( M_rhs, M_rhs_pressure, *u_star, *rhs_velocity);
624 
625  M_rhs_pressure->globalAssemble();
626 
627  chrono.stop();
628  M_displayer.leaderPrintMax ( " done in ", chrono.diff() ) ;
629  }
630 
631  if ( M_imposeWeakBC )
632  {
633  // Reference for essential BCs imposed weakly:
634  // - Y. Bazilevs, K. Takizawa and E. Tezduyar. Computational Fluid-Structre Interaction. Page 70.
635  using namespace ExpressionAssembly;
636 
637  QuadratureBoundary myBDQR (buildTetraBDQR (quadRuleTria7pt) );
638 
639  M_block01_weakBC.reset( new matrix_Type(M_velocityFESpace->map() ) );
640 
641  integrate( boundary(M_fespaceUETA->mesh(), M_flagWeakBC),
642  myBDQR,
643  M_fespaceUETA,
644  M_fespacePETA,
645  dot( phi_i, phi_j * Nface)
646  )
647  >> M_block01_weakBC;
648  M_block01_weakBC->globalAssemble( M_pressureFESpace->mapPtr(), M_velocityFESpace->mapPtr() );
649  *M_block01 += *M_block01_weakBC;
650 
651  integrate( boundary(M_fespaceUETA->mesh(), M_flagWeakBC),
652  myBDQR,
653  M_fespacePETA,
654  M_fespaceUETA,
655  value(-1.0)*dot(phi_i*Nface, phi_j)
656  )
657  >> M_block10;
658  }
659 
660  M_rhs->globalAssemble();
661  M_block00->globalAssemble();
662  M_block01->globalAssemble( M_pressureFESpace->mapPtr(), M_velocityFESpace->mapPtr() );
663  M_block10->globalAssemble(M_velocityFESpace->mapPtr(), M_pressureFESpace->mapPtr());
664 
665 }
666 
667 void NavierStokesSolverBlocks::applyBoundaryConditions ( bcPtr_Type & bc, const Real& time )
668 {
669  if ( M_computeAerodynamicLoads )
670  {
671  M_displayer.leaderPrint( "\tNS operator - Compute Loads: TRUE");
672 
673  M_block00_noBC.reset( new matrix_Type(M_velocityFESpace->map() ) );
674  M_block01_noBC.reset( new matrix_Type(M_velocityFESpace->map() ) );
675 
676  M_block00_noBC->zero();
677  M_block01_noBC->zero();
678 
679  *M_block00_noBC += *M_block00;
680  *M_block01_noBC += *M_block01;
681 
682  if ( M_imposeWeakBC )
683  {
684  *M_block00_noBC -= *M_block00_weakBC;
685  *M_block01_noBC -= *M_block01_weakBC;
686  }
687 
688  M_block00_noBC->globalAssemble();
689  M_block01_noBC->globalAssemble( M_pressureFESpace->mapPtr(), M_velocityFESpace->mapPtr() );
690 
691  // Right hand side
692  M_rhs_noBC.reset( new vector_Type ( M_velocityFESpace->map(), Unique ) );
693  M_rhs_noBC->zero();
694  *M_rhs_noBC += *M_rhs;
695  }
696 
697  updateBCHandler(bc);
698  bcManage ( *M_block00, *M_rhs, *M_velocityFESpace->mesh(), M_velocityFESpace->dof(), *bc, M_velocityFESpace->feBd(), 1.0, time );
699  bcManageMatrix( *M_block01, *M_velocityFESpace->mesh(), M_velocityFESpace->dof(), *bc, M_velocityFESpace->feBd(), 0.0, 0.0);
700 }
701 
702 void NavierStokesSolverBlocks::applyBoundaryConditions ( bcPtr_Type & bc, const Real& time, const vectorPtr_Type& velocities )
703 {
704  /* Used only for example aorta semi implicit */
705 
706  updateBCHandler(bc);
707  bcManage ( *M_block00, *M_rhs, *M_velocityFESpace->mesh(), M_velocityFESpace->dof(), *bc, M_velocityFESpace->feBd(), 1.0, time );
708  bcManageMatrix( *M_block01, *M_velocityFESpace->mesh(), M_velocityFESpace->dof(), *bc, M_velocityFESpace->feBd(), 0.0, 0.0);
709 
710  *M_rhs += *velocities;
711 }
712 
713 void NavierStokesSolverBlocks::iterate( bcPtr_Type & bc, const Real& time, const vectorPtr_Type& velocities )
714 {
715  applyBoundaryConditions ( bc, time, velocities);
716  solveTimeStep();
717  if ( M_dataFile("fluid/stabilization/ode_fine_scale", false ) )
718  {
719  if ( M_stabilizationType.compare("SUPG_SEMI_IMPLICIT")==0 )
720  M_stabilization->updateODEfineScale ( M_velocity, M_pressure );
721  else if ( M_stabilizationType.compare("VMSLES_NEW")==0 )
722  M_stabilization->updateODEfineScale ( M_velocity, M_pressure, M_velocityExtrapolated );
723  }
724 }
725 
726 void NavierStokesSolverBlocks::iterate( bcPtr_Type & bc, const Real& time )
727 {
728  applyBoundaryConditions ( bc, time );
729  solveTimeStep();
730  if (M_useStabilization)
731  if ( M_dataFile("fluid/stabilization/ode_fine_scale", false ) )
732  {
733  if ( M_stabilizationType.compare("SUPG_SEMI_IMPLICIT")==0 )
734  M_stabilization->updateODEfineScale ( M_velocity, M_pressure );
735  else if ( M_stabilizationType.compare("VMSLES_NEW")==0 )
736  M_stabilization->updateODEfineScale ( M_velocity, M_pressure, M_velocityExtrapolated );
737  }
738 }
739 
740 void NavierStokesSolverBlocks::solveTimeStep( )
741 {
742  //(1) Set up the OseenOperator
743  M_displayer.leaderPrint( "\tNS operator - set up the block operator...");
744  LifeChrono chrono;
745  chrono.start();
746 
747  Operators::NavierStokesOperator::operatorPtrContainer_Type operData(2,2);
748  operData(0,0) = M_block00->matrixPtr();
749  operData(0,1) = M_block01->matrixPtr();
750  operData(1,0) = M_block10->matrixPtr();
751  if ( M_useStabilization )
752  operData(1,1) = M_block11->matrixPtr();
753 
754  M_oper->setUp(operData, M_displayer.comm());
755  chrono.stop();
756  M_displayer.leaderPrintMax(" done in " , chrono.diff() );
757 
758  //(2) Set the data for the preconditioner
759 
760  M_displayer.leaderPrint( "\tPreconditioner operator - set up the block operator...");
761  chrono.reset();
762  chrono.start();
763 
764  if ( std::strcmp(M_prec->Label(),"aSIMPLEOperator")==0 )
765  {
766  if (M_useStabilization)
767  {
768  M_prec->setUp(M_block00, M_block10, M_block01, M_block11);
769  }
770  else
771  {
772  M_prec->setUp(M_block00, M_block10, M_block01);
773  }
774 
775  M_prec->setDomainMap(M_oper->OperatorDomainBlockMapPtr());
776  M_prec->setRangeMap(M_oper->OperatorRangeBlockMapPtr());
777  M_prec->updateApproximatedMomentumOperator();
778  M_prec->updateApproximatedSchurComplementOperator();
779  }
780 
781  chrono.stop();
782  M_displayer.leaderPrintMax(" done in " , chrono.diff() );
783 
784  //(3) Set the solver for the linear system
785  M_displayer.leaderPrint( "\tset up the Trilinos solver...");
786  chrono.start();
787  std::string solverType(M_pListLinSolver->get<std::string>("Linear Solver Type"));
788  M_invOper.reset(Operators::InvertibleOperatorFactory::instance().createObject(solverType));
789 
790  M_invOper->setParameterList(M_pListLinSolver->sublist(solverType));
791  M_invOper->setOperator(M_oper);
792  M_invOper->setPreconditioner(M_prec);
793 
794  chrono.stop();
795  M_displayer.leaderPrintMax(" done in " , chrono.diff() );
796 
797  // Solving the system
798  BlockEpetra_Map upMap;
799  upMap.setUp ( M_velocityFESpace->map().map(Unique), M_pressureFESpace->map().map(Unique));
800 
801  if ( M_useStabilization )
802  {
803  vector_Type rhs ( *M_monolithicMap, Unique );
804  vector_Type sol ( *M_monolithicMap, Unique );
805 
806  rhs.zero();
807  sol.zero();
808 
809  rhs.subset ( *M_rhs, M_velocityFESpace->map(), 0, 0 );
810  rhs.subset ( *M_rhs_pressure, M_pressureFESpace->map(), 0, M_velocityFESpace->map().mapSize() );
811 
812  M_invOper->ApplyInverse(rhs.epetraVector(), sol.epetraVector());
813 
814  M_velocity->subset ( sol, M_velocityFESpace->map(), 0, 0 );
815  M_pressure->subset ( sol, M_pressureFESpace->map(), M_velocityFESpace->map().mapSize(), 0 );
816  }
817  else
818  {
819  BlockEpetra_MultiVector up(upMap, 1), rhs(upMap, 1);
820  rhs.block(0).Update(1.0, M_rhs->epetraVector(), 0.);
821 
822  // Solving the linear system
823  M_invOper->ApplyInverse(rhs,up);
824 
825  M_velocity->epetraVector().Update(1.0,up.block(0),0.0);
826  M_pressure->epetraVector().Update(1.0,up.block(1),0.0);
827  }
828 
829 
830  if ( M_computeAerodynamicLoads )
831  {
832  // Computation of aerodynamic forces in Residual form, see [Forti and Dede, 2015]
833  M_forces.reset ( new vector_Type ( M_velocityFESpace->map(), Unique ) );
834 
835  M_forces->zero();
836 
837  *M_forces += *M_rhs_noBC;
838 
839  *M_forces -= *M_block00_noBC * ( *M_velocity );
840 
841  *M_forces -= *M_block01_noBC * ( *M_pressure );
842  }
843 }
844 
845 VectorSmall<2> NavierStokesSolverBlocks::computeForces( BCHandler& bcHDrag,
846  BCHandler& bcHLift )
847 {
848  bcHDrag.bcUpdate ( *M_velocityFESpace->mesh(), M_velocityFESpace->feBd(), M_velocityFESpace->dof() );
849  bcHLift.bcUpdate ( *M_velocityFESpace->mesh(), M_velocityFESpace->feBd(), M_velocityFESpace->dof() );
850 
851  vector_Type onesOnBodyDrag(M_velocityFESpace->map(), Unique);
852  onesOnBodyDrag.zero();
853 
854  vector_Type onesOnBodyLift(M_velocityFESpace->map(), Unique);
855  onesOnBodyLift.zero();
856 
857  bcManageRhs ( onesOnBodyDrag, *M_velocityFESpace->mesh(), M_velocityFESpace->dof(), bcHDrag, M_velocityFESpace->feBd(), 1., 0.);
858  bcManageRhs ( onesOnBodyLift, *M_velocityFESpace->mesh(), M_velocityFESpace->dof(), bcHLift, M_velocityFESpace->feBd(), 1., 0.);
859 
860  Real drag (0.0);
861  Real lift (0.0);
862 
863  drag = M_forces->dot(onesOnBodyDrag);
864  lift = M_forces->dot(onesOnBodyLift);
865 
866  VectorSmall<2> Forces;
867  Forces[0] = drag;
868  Forces[1] = lift;
869 
870  return Forces;
871 }
872 
873 void NavierStokesSolverBlocks::iterate_nonlinear( const Real& time )
874 {
875  applyBoundaryConditionsSolution ( time );
876 
877  M_time = time;
878 
879  // Call Newton
880  UInt status = NonLinearRichardson ( *M_solution, *this, M_absoluteTolerance, M_relativeTolerance, M_maxiterNonlinear, M_etaMax,
881  M_nonLinearLineSearch, 0, 2, M_out_res, 0.0);
882 
883  if ( M_computeAerodynamicLoads )
884  {
885  M_forces.reset ( new vector_Type ( *M_monolithicMap, Unique ) );
886  M_forces->zero();
887  computeForcesNonLinear(M_forces, M_solution);
888  }
889 
890  if (status == EXIT_FAILURE)
891  {
892  M_displayer.leaderPrint(" WARNING: Newton failed " );
893  }
894 
895  M_velocity->subset ( *M_solution, M_velocityFESpace->map(), 0, 0 );
896  M_pressure->subset ( *M_solution, M_pressureFESpace->map(), M_velocityFESpace->map().mapSize(), 0 );
897 }
898 
899 void NavierStokesSolverBlocks::iterate_steady( )
900 {
901  // Initialize the solution
902  M_solution->zero();
903 
904  applyBoundaryConditionsSolution ( 0.0 ); // the second argument is zero since the problem is steady
905 
906  // Call Newton
907  UInt status;
908  status = NonLinearRichardson ( *M_solution, *this, M_absoluteTolerance, M_relativeTolerance, M_maxiterNonlinear, M_etaMax,
909  M_nonLinearLineSearch, 0, 2, M_out_res, 0.0);
910 
911  M_velocity->subset ( *M_solution, M_velocityFESpace->map(), 0, 0 );
912  M_pressure->subset ( *M_solution, M_pressureFESpace->map(), M_velocityFESpace->map().mapSize(), 0 );
913 }
914 
915 void NavierStokesSolverBlocks::applyBoundaryConditionsSolution ( const Real& time )
916 {
917  updateBCHandler(M_bc_sol);
918  bcManageRhs ( *M_solution, *M_velocityFESpace->mesh(), M_velocityFESpace->dof(), *M_bc_sol, M_velocityFESpace->feBd(), 1.0, time );
919 }
920 
921 void NavierStokesSolverBlocks::updateBCHandler( bcPtr_Type & bc )
922 {
923  bc->bcUpdate ( *M_velocityFESpace->mesh(), M_velocityFESpace->feBd(), M_velocityFESpace->dof() );
924 }
925 
926 void NavierStokesSolverBlocks::evaluateResidual( const vectorPtr_Type& convective_velocity,
927  const vectorPtr_Type& velocity_km1,
928  const vectorPtr_Type& pressure_km1,
929  const vectorPtr_Type& rhs_velocity,
930  vectorPtr_Type& residualVelocity,
931  vectorPtr_Type& residualPressure)
932 {
933  residualVelocity->zero();
934  residualPressure->zero();
935 
936  // Residual vector for the velocity and pressure components
937  vectorPtr_Type res_velocity ( new vector_Type ( M_velocityFESpace->map(), Repeated ) );
938  vectorPtr_Type res_pressure ( new vector_Type ( M_pressureFESpace->map(), Repeated ) );
939  res_velocity->zero();
940  res_pressure->zero();
941 
942  // Get repeated versions of input vectors for the assembly
943  vectorPtr_Type convective_velocity_repeated ( new vector_Type (*convective_velocity, Repeated) );
944  vectorPtr_Type u_km1_repeated ( new vector_Type (*velocity_km1, Repeated) );
945  vectorPtr_Type p_km1_repeated ( new vector_Type (*pressure_km1, Repeated) );
946  vectorPtr_Type rhs_velocity_repeated ( new vector_Type (*rhs_velocity, Repeated) );
947 
948  {
949  using namespace ExpressionAssembly;
950 
951  if ( M_stiffStrain )
952  {
953  integrate ( elements ( M_fespaceUETA->mesh() ),
954  M_velocityFESpace->qr(),
955  M_fespaceUETA,
956  value(M_density) * value ( M_alpha / M_timeStep ) * dot ( value ( M_fespaceUETA, *u_km1_repeated ), phi_i ) -
957  value(M_density) * dot ( value ( M_fespaceUETA, *rhs_velocity_repeated ), phi_i ) +
958  value(M_viscosity) * dot ( grad ( phi_i ) + transpose ( grad ( phi_i ) ), grad ( M_fespaceUETA, *u_km1_repeated ) + transpose ( grad ( M_fespaceUETA, *u_km1_repeated ) ) ) +
959  value(M_density) * dot ( value ( M_fespaceUETA, *convective_velocity_repeated ) * grad ( M_fespaceUETA, *u_km1_repeated ), phi_i ) +
960  value ( -1.0 ) * value ( M_fespacePETA, *p_km1_repeated ) * div ( phi_i )
961  ) >> res_velocity;
962  }
963  else
964  {
965  integrate ( elements ( M_fespaceUETA->mesh() ),
966  M_velocityFESpace->qr(),
967  M_fespaceUETA,
968  value(M_density) * value ( M_alpha / M_timeStep ) * dot ( value ( M_fespaceUETA, *u_km1_repeated ), phi_i ) -
969  value(M_density) * dot ( value ( M_fespaceUETA, *rhs_velocity_repeated ), phi_i ) +
970  value(M_viscosity) * dot ( grad ( phi_i ), grad ( M_fespaceUETA, *u_km1_repeated ) + transpose ( grad ( M_fespaceUETA, *u_km1_repeated ) ) ) +
971  value(M_density) * dot ( value ( M_fespaceUETA, *convective_velocity_repeated ) * grad ( M_fespaceUETA, *u_km1_repeated ), phi_i ) +
972  value ( -1.0 ) * value ( M_fespacePETA, *p_km1_repeated ) * div ( phi_i )
973  ) >> res_velocity;
974  }
975 
976  integrate ( elements ( M_fespaceUETA->mesh() ),
977  M_pressureFESpace->qr(),
978  M_fespacePETA,
979  trace ( grad ( M_fespaceUETA, *u_km1_repeated ) ) * phi_i
980  ) >> res_pressure;
981  }
982 
983  if ( M_useStabilization )
984  {
985  M_displayer.leaderPrint ( "[F] - Assembly residual of the stabilization \n" ) ;
986  M_stabilization->apply_vector(res_velocity, res_pressure, *convective_velocity, *velocity_km1, *pressure_km1, *rhs_velocity);
987  }
988 
989  res_velocity->globalAssemble();
990  res_pressure->globalAssemble();
991 
992  vector_Type res_velocity_unique ( *res_velocity, Unique );
993  vector_Type res_pressure_unique ( *res_pressure, Unique );
994 
995  *residualVelocity = res_velocity_unique;
996  *residualPressure = res_pressure_unique;
997 }
998 
999 void NavierStokesSolverBlocks::evaluateResidual( const vectorPtr_Type& convective_velocity,
1000  const vectorPtr_Type& velocity_km1,
1001  const vectorPtr_Type& pressure_km1,
1002  const vectorPtr_Type& rhs_velocity,
1003  vectorPtr_Type& residual)
1004 {
1005  // This methos is used in FSI to assemble the fluid residual component
1006  residual->zero();
1007 
1008  // Residual vector for the velocity and pressure components
1009  vectorPtr_Type res_velocity ( new vector_Type ( M_velocityFESpace->map(), Repeated ) );
1010  vectorPtr_Type res_pressure ( new vector_Type ( M_pressureFESpace->map(), Repeated ) );
1011  res_velocity->zero();
1012  res_pressure->zero();
1013 
1014  // Get repeated versions of input vectors for the assembly
1015  vectorPtr_Type convective_velocity_repeated ( new vector_Type (*convective_velocity, Repeated) );
1016  vectorPtr_Type u_km1_repeated ( new vector_Type (*velocity_km1, Repeated) );
1017  vectorPtr_Type p_km1_repeated ( new vector_Type (*pressure_km1, Repeated) );
1018  vectorPtr_Type rhs_velocity_repeated ( new vector_Type (*rhs_velocity, Repeated) );
1019 
1020  {
1021  using namespace ExpressionAssembly;
1022 
1023  if ( M_stiffStrain )
1024  {
1025  integrate ( elements ( M_fespaceUETA->mesh() ),
1026  M_velocityFESpace->qr(),
1027  M_fespaceUETA,
1028  M_density * value ( M_alpha / M_timeStep ) * dot ( value ( M_fespaceUETA, *u_km1_repeated ), phi_i ) -
1029  M_density * dot ( value ( M_fespaceUETA, *rhs_velocity_repeated ), phi_i ) +
1030  M_viscosity * dot ( grad ( phi_i ) + transpose ( grad ( phi_i ) ), grad ( M_fespaceUETA, *u_km1_repeated ) + transpose ( grad ( M_fespaceUETA, *u_km1_repeated ) ) ) +
1031  M_density * dot ( value ( M_fespaceUETA, *convective_velocity_repeated ) * grad ( M_fespaceUETA, *u_km1_repeated ), phi_i ) +
1032  value ( -1.0 ) * value ( M_fespacePETA, *p_km1_repeated ) * div ( phi_i )
1033  ) >> res_velocity;
1034  }
1035  else
1036  {
1037  integrate ( elements ( M_fespaceUETA->mesh() ),
1038  M_velocityFESpace->qr(),
1039  M_fespaceUETA,
1040  M_density * value ( M_alpha / M_timeStep ) * dot ( value ( M_fespaceUETA, *u_km1_repeated ), phi_i ) -
1041  M_density * dot ( value ( M_fespaceUETA, *rhs_velocity_repeated ), phi_i ) +
1042  M_viscosity * dot ( grad ( phi_i ), grad ( M_fespaceUETA, *u_km1_repeated ) + transpose ( grad ( M_fespaceUETA, *u_km1_repeated ) ) ) +
1043  M_density * dot ( value ( M_fespaceUETA, *convective_velocity_repeated ) * grad ( M_fespaceUETA, *u_km1_repeated ), phi_i ) +
1044  value ( -1.0 ) * value ( M_fespacePETA, *p_km1_repeated ) * div ( phi_i )
1045  ) >> res_velocity;
1046  }
1047 
1048  integrate ( elements ( M_fespaceUETA->mesh() ),
1049  M_pressureFESpace->qr(),
1050  M_fespacePETA,
1051  trace ( grad ( M_fespaceUETA, *u_km1_repeated ) ) * phi_i
1052  ) >> res_pressure;
1053  }
1054 
1055  if ( M_useStabilization )
1056  {
1057  M_displayer.leaderPrint ( "[F] - Assembly residual of the stabilization \n" ) ;
1058  M_stabilization->apply_vector(res_velocity, res_pressure, *convective_velocity, *velocity_km1, *pressure_km1, *rhs_velocity);
1059  }
1060 
1061  res_velocity->globalAssemble();
1062  res_pressure->globalAssemble();
1063 
1064  vector_Type res_velocity_unique ( *res_velocity, Unique );
1065  vector_Type res_pressure_unique ( *res_pressure, Unique );
1066 
1067  residual->subset ( res_velocity_unique, M_velocityFESpace->map(), 0, 0 );
1068  residual->subset ( res_pressure_unique, M_pressureFESpace->map(), 0, M_velocityFESpace->map().mapSize() );
1069 
1070 }
1071 
1072 void NavierStokesSolverBlocks::assembleInterfaceMass( matrixPtr_Type& mass_interface, const mapPtr_Type& interface_map,
1073  markerID_Type interfaceFlag, const vectorPtr_Type& numerationInterface,
1074  const UInt& offset)
1075 {
1076  // INITIALIZE MATRIX WITH THE MAP OF THE INTERFACE
1077  matrixPtr_Type fluid_interfaceMass( new matrix_Type ( M_velocityFESpace->map(), 50 ) );
1078  fluid_interfaceMass->zero();
1079 
1080  // ASSEMBLE MASS MATRIX AT THE INTERFACE
1081  QuadratureBoundary myBDQR (buildTetraBDQR (quadRuleTria7pt) );
1082  {
1083  using namespace ExpressionAssembly;
1084  integrate ( boundary (M_fespaceUETA->mesh(), interfaceFlag ),
1085  myBDQR,
1086  M_fespaceUETA,
1087  M_fespaceUETA,
1088  //Boundary Mass
1089  dot(phi_i,phi_j)
1090  )
1091  >> fluid_interfaceMass;
1092  }
1093  fluid_interfaceMass->globalAssemble();
1094 
1095  // RESTRICT MATRIX TO INTERFACE DOFS ONLY
1096  mass_interface.reset(new matrix_Type ( *interface_map, 50 ) );
1097  fluid_interfaceMass->restrict ( interface_map, numerationInterface, offset, mass_interface );
1098 }
1099 
1100 void NavierStokesSolverBlocks::computeForcesNonLinear(vectorPtr_Type& force, const vectorPtr_Type& solution)
1101 {
1102  // Forces to zero
1103  force->zero();
1104 
1105  // Extract the velocity and the pressure from the solution vector
1106  vectorPtr_Type u_km1 ( new vector_Type ( M_velocityFESpace->map(), Repeated ) );
1107  vectorPtr_Type p_km1 ( new vector_Type ( M_pressureFESpace->map(), Repeated ) );
1108  u_km1->zero();
1109  p_km1->zero();
1110  u_km1->subset ( *solution, M_velocityFESpace->map(), 0, 0 );
1111  p_km1->subset ( *solution, M_pressureFESpace->map(), M_velocityFESpace->map().mapSize(), 0 );
1112 
1113  // Force vector for the velocity and pressure components
1114  vectorPtr_Type res_velocity ( new vector_Type ( M_velocityFESpace->map(), Unique ) );
1115  vectorPtr_Type res_pressure ( new vector_Type ( M_pressureFESpace->map(), Unique ) );
1116  res_velocity->zero();
1117  res_pressure->zero();
1118 
1119  vectorPtr_Type rhs_velocity_repeated;
1120  rhs_velocity_repeated.reset( new vector_Type ( *M_velocityRhs, Repeated ) );
1121 
1122  using namespace ExpressionAssembly;
1123 
1124  if ( M_stiffStrain )
1125  {
1126  integrate ( elements ( M_fespaceUETA->mesh() ),
1127  M_velocityFESpace->qr(),
1128  M_fespaceUETA,
1129  M_density * value ( M_alpha / M_timeStep ) * dot ( value ( M_fespaceUETA, *u_km1 ), phi_i ) -
1130  M_density * dot ( value ( M_fespaceUETA, *rhs_velocity_repeated ), phi_i ) +
1131  value ( 0.5 ) * M_viscosity * dot ( grad ( phi_i ) + transpose ( grad ( phi_i ) ), grad ( M_fespaceUETA, *u_km1 ) +
1132  transpose ( grad ( M_fespaceUETA, *u_km1 ) ) ) +
1133  M_density * dot ( value ( M_fespaceUETA, *u_km1 ) * grad ( M_fespaceUETA, *u_km1 ), phi_i ) +
1134  value ( -1.0 ) * value ( M_fespacePETA, *p_km1 ) * div ( phi_i )
1135  ) >> res_velocity;
1136  }
1137  else
1138  {
1139  integrate ( elements ( M_fespaceUETA->mesh() ),
1140  M_velocityFESpace->qr(),
1141  M_fespaceUETA,
1142  M_density * value ( M_alpha / M_timeStep ) * dot ( value ( M_fespaceUETA, *u_km1 ), phi_i ) -
1143  M_density * dot ( value ( M_fespaceUETA, *rhs_velocity_repeated ), phi_i ) +
1144  M_viscosity * dot ( grad ( phi_i ), grad ( M_fespaceUETA, *u_km1 ) + transpose ( grad ( M_fespaceUETA, *u_km1 ) ) ) +
1145  M_density * dot ( value ( M_fespaceUETA, *u_km1 ) * grad ( M_fespaceUETA, *u_km1 ), phi_i ) +
1146  value ( -1.0 ) * value ( M_fespacePETA, *p_km1 ) * div ( phi_i )
1147  ) >> res_velocity;
1148  }
1149 
1150  integrate ( elements ( M_fespaceUETA->mesh() ),
1151  M_pressureFESpace->qr(),
1152  M_fespacePETA,
1153  trace ( grad ( M_fespaceUETA, *u_km1 ) ) * phi_i
1154  ) >> res_pressure;
1155 
1156  if ( M_useStabilization )
1157  {
1158  M_stabilization->apply_vector(res_velocity, res_pressure, *u_km1, *p_km1, *rhs_velocity_repeated);
1159  }
1160 
1161  res_velocity->globalAssemble();
1162  res_pressure->globalAssemble();
1163 
1164  force->subset ( *res_velocity, M_velocityFESpace->map(), 0, 0 );
1165  force->subset ( *res_pressure, M_pressureFESpace->map(), 0, M_velocityFESpace->map().mapSize() );
1166 }
1167 
1168 void NavierStokesSolverBlocks::evalResidual(vector_Type& residual, const vector_Type& solution, const UInt /*iter_newton*/ )
1169 {
1170  // Residual to zero
1171  residual.zero();
1172 
1173  // Extract the velocity and the pressure at the previous Newton step ( index k minus 1, i.e. km1)
1174  vectorPtr_Type u_km1_subs ( new vector_Type ( M_velocityFESpace->map(), Unique ) );
1175  vectorPtr_Type p_km1_subs ( new vector_Type ( M_pressureFESpace->map(), Unique ) );
1176  u_km1_subs->zero();
1177  p_km1_subs->zero();
1178  u_km1_subs->subset ( solution, M_velocityFESpace->map(), 0, 0 );
1179  p_km1_subs->subset ( solution, M_pressureFESpace->map(), M_velocityFESpace->map().mapSize(), 0 );
1180 
1181  // Repeated vectors
1182  vectorPtr_Type u_km1 ( new vector_Type ( *u_km1_subs, Repeated ) );
1183  vectorPtr_Type p_km1 ( new vector_Type ( *p_km1_subs, Repeated ) );
1184 
1185  // Residual vector for the velocity and pressure components
1186  vectorPtr_Type res_velocity ( new vector_Type ( M_velocityFESpace->map(), Repeated ) );
1187  vectorPtr_Type res_pressure ( new vector_Type ( M_pressureFESpace->map(), Repeated ) );
1188  res_velocity->zero();
1189  res_pressure->zero();
1190 
1191  vectorPtr_Type rhs_velocity_repeated;
1192 
1193  if ( !M_steady )
1194  {
1195  rhs_velocity_repeated.reset( new vector_Type ( *M_velocityRhs, Repeated ) );
1196  }
1197 
1198  if ( M_steady )
1199  {
1200  // Assemble the residual vector:
1201  using namespace ExpressionAssembly;
1202 
1203  if ( M_stiffStrain )
1204  {
1205  integrate ( elements ( M_fespaceUETA->mesh() ),
1206  M_velocityFESpace->qr(),
1207  M_fespaceUETA,
1208  value ( 0.5 ) * M_viscosity * dot ( grad ( phi_i ) + transpose ( grad ( phi_i ) ), grad ( M_fespaceUETA, *u_km1 ) + transpose ( grad ( M_fespaceUETA, *u_km1 ) ) ) +
1209  M_density * dot ( value ( M_fespaceUETA, *u_km1 ) * grad ( M_fespaceUETA, *u_km1 ), phi_i ) +
1210  value ( -1.0 ) * value ( M_fespacePETA, *p_km1 ) * div ( phi_i )
1211 
1212  ) >> res_velocity;
1213  }
1214  else
1215  {
1216  integrate ( elements ( M_fespaceUETA->mesh() ),
1217  M_velocityFESpace->qr(),
1218  M_fespaceUETA,
1219  M_viscosity * dot ( grad ( phi_i ), grad ( M_fespaceUETA, *u_km1 ) + transpose ( grad ( M_fespaceUETA, *u_km1 ) ) ) +
1220  M_density * dot ( value ( M_fespaceUETA, *u_km1 ) * grad ( M_fespaceUETA, *u_km1 ), phi_i ) +
1221  value ( -1.0 ) * value ( M_fespacePETA, *p_km1 ) * div ( phi_i )
1222  ) >> res_velocity;
1223  }
1224 
1225  integrate ( elements ( M_fespaceUETA->mesh() ),
1226  M_pressureFESpace->qr(),
1227  M_fespacePETA,
1228  trace ( grad ( M_fespaceUETA, *u_km1 ) ) * phi_i
1229  ) >> res_pressure;
1230  }
1231  else
1232  {
1233  using namespace ExpressionAssembly;
1234 
1235  if ( M_stiffStrain )
1236  {
1237  integrate ( elements ( M_fespaceUETA->mesh() ),
1238  M_velocityFESpace->qr(),
1239  M_fespaceUETA,
1240  M_density * value ( M_alpha / M_timeStep ) * dot ( value ( M_fespaceUETA, *u_km1 ), phi_i ) -
1241  M_density * dot ( value ( M_fespaceUETA, *rhs_velocity_repeated ), phi_i ) +
1242  M_viscosity * dot ( grad ( phi_i ) + transpose ( grad ( phi_i ) ), grad ( M_fespaceUETA, *u_km1 ) + transpose ( grad ( M_fespaceUETA, *u_km1 ) ) ) +
1243  M_density * dot ( value ( M_fespaceUETA, *u_km1 ) * grad ( M_fespaceUETA, *u_km1 ), phi_i ) +
1244  value ( -1.0 ) * value ( M_fespacePETA, *p_km1 ) * div ( phi_i )
1245  ) >> res_velocity;
1246  }
1247  else
1248  {
1249  integrate ( elements ( M_fespaceUETA->mesh() ),
1250  M_velocityFESpace->qr(),
1251  M_fespaceUETA,
1252  M_density * value ( M_alpha / M_timeStep ) * dot ( value ( M_fespaceUETA, *u_km1 ), phi_i ) -
1253  M_density * dot ( value ( M_fespaceUETA, *rhs_velocity_repeated ), phi_i ) +
1254  M_viscosity * dot ( grad ( phi_i ), grad ( M_fespaceUETA, *u_km1 ) + transpose ( grad ( M_fespaceUETA, *u_km1 ) ) ) +
1255  M_density * dot ( value ( M_fespaceUETA, *u_km1 ) * grad ( M_fespaceUETA, *u_km1 ), phi_i ) +
1256  value ( -1.0 ) * value ( M_fespacePETA, *p_km1 ) * div ( phi_i )
1257  ) >> res_velocity;
1258  }
1259 
1260  integrate ( elements ( M_fespaceUETA->mesh() ),
1261  M_pressureFESpace->qr(),
1262  M_fespacePETA,
1263  trace ( grad ( M_fespaceUETA, *u_km1 ) ) * phi_i
1264  ) >> res_pressure;
1265  }
1266 
1267  if ( M_useStabilization )
1268  {
1269  M_displayer.leaderPrint ( "[F] - Assembly residual of the stabilization \n" ) ;
1270  M_stabilization->apply_vector(res_velocity, res_pressure, *u_km1, *p_km1, *rhs_velocity_repeated);
1271  }
1272 
1273  res_velocity->globalAssemble();
1274  res_pressure->globalAssemble();
1275 
1276  // Residual vector for the velocity and pressure components
1277  vector_Type res_velocity_unique ( *res_velocity, Unique );
1278  vector_Type res_pressure_unique ( *res_pressure, Unique );
1279 
1280  if ( M_steady )
1281  {
1282  applyBoundaryConditionsResidual(res_velocity_unique);
1283  }
1284  else
1285  {
1286  applyBoundaryConditionsResidual(res_velocity_unique, M_time );
1287  }
1288 
1289  residual.subset ( res_velocity_unique, M_velocityFESpace->map(), 0, 0 );
1290  residual.subset ( res_pressure_unique, M_pressureFESpace->map(), 0, M_velocityFESpace->map().mapSize() );
1291 
1292  // We need to update the linearized terms in the jacobian coming from the convective part
1293 
1294  updateConvectiveTerm(u_km1);
1295  updateJacobian(*u_km1);
1296 
1297  if ( M_useStabilization )
1298  {
1299  M_stabilization->apply_matrix(*u_km1, *p_km1, *M_velocityRhs);
1300 
1301  *M_block00 += *M_stabilization->block_00();
1302  M_block00->globalAssemble();
1303 
1304  M_block01->zero();
1305  *M_block01 += *M_Btranspose;
1306  *M_block01 += *M_stabilization->block_01();
1307  M_block01->globalAssemble( M_pressureFESpace->mapPtr(), M_velocityFESpace->mapPtr() );
1308 
1309  M_block10->zero();
1310  *M_block10 += *M_B;
1311  *M_block10 += *M_stabilization->block_10();
1312  M_block10->globalAssemble(M_velocityFESpace->mapPtr(), M_pressureFESpace->mapPtr());
1313 
1314  M_block11->zero();
1315  *M_block11 += *M_stabilization->block_11();
1316  M_block11->globalAssemble();
1317  }
1318 
1319 }
1320 
1321 void NavierStokesSolverBlocks::applyBoundaryConditionsResidual ( vector_Type& r_u, const Real& time )
1322 {
1323  //! Extract each component of the input vector
1324  VectorEpetra ru_copy(r_u, Unique);
1325 
1326  //! Apply BC on each component
1327  if ( !M_bc_res_essential->bcUpdateDone() )
1328  M_bc_res_essential->bcUpdate ( *M_velocityFESpace->mesh(),
1329  M_velocityFESpace->feBd(),
1330  M_velocityFESpace->dof() );
1331 
1332  if ( !M_bc_res_natural->bcUpdateDone() )
1333  M_bc_res_natural->bcUpdate ( *M_velocityFESpace->mesh(),
1334  M_velocityFESpace->feBd(),
1335  M_velocityFESpace->dof() );
1336 
1337  bcManageRhs ( ru_copy,
1338  *M_velocityFESpace->mesh(),
1339  M_velocityFESpace->dof(),
1340  *M_bc_res_natural,
1341  M_velocityFESpace->feBd(),
1342  1.0,
1343  time );
1344 
1345  bcManageRhs ( ru_copy,
1346  *M_velocityFESpace->mesh(),
1347  M_velocityFESpace->dof(),
1348  *M_bc_res_essential,
1349  M_velocityFESpace->feBd(),
1350  0.0,
1351  0.0 );
1352 
1353  r_u.zero();
1354  r_u = ru_copy;
1355 }
1356 
1357 void NavierStokesSolverBlocks::updateConvectiveTerm ( const vectorPtr_Type& velocity)
1358 {
1359  // Note that u_star HAS to extrapolated from outside. Hence it works also for FSI in this manner.
1360  vectorPtr_Type velocity_repeated;
1361  velocity_repeated.reset ( new vector_Type ( *velocity, Repeated ) );
1362 
1363  // Update convective term
1364  M_C->zero();
1365  {
1366  if ( M_useFastAssembly )
1367  {
1368  M_fastAssembler->assembleConvective( M_C, *velocity_repeated );
1369  }
1370  else
1371  {
1372  using namespace ExpressionAssembly;
1373  integrate ( elements (M_fespaceUETA->mesh() ),
1374  M_velocityFESpace->qr(),
1375  M_fespaceUETA,
1376  M_fespaceUETA,
1377  dot ( M_density * value (M_fespaceUETA, *velocity_repeated) *grad (phi_j), phi_i)
1378  )
1379  >> M_C;
1380  }
1381  }
1382  M_C->globalAssemble();
1383 
1384  M_block00->zero();
1385  if ( !M_steady )
1386  {
1387  *M_block00 += *M_Mu;
1388  *M_block00 *= M_alpha / M_timeStep;
1389  }
1390  *M_block00 += *M_A;
1391  *M_block00 += *M_C;
1392  if ( !M_fullyImplicit )
1393  M_block00->globalAssemble( );
1394 }
1395 
1396 void NavierStokesSolverBlocks::updateJacobian( const vector_Type& u_k )
1397 {
1398  vector_Type uk_rep ( u_k, Repeated );
1399  M_Jacobian->zero();
1400 
1401  M_displayer.leaderPrint ( "[F] - Update Jacobian convective term\n" ) ;
1402 
1403  if ( M_fullyImplicit )
1404  {
1405  if ( M_useFastAssembly )
1406  {
1407  M_fastAssembler->jacobianNS( M_Jacobian, uk_rep );
1408  }
1409  else
1410  {
1411  using namespace ExpressionAssembly;
1412  integrate( elements(M_fespaceUETA->mesh()),
1413  M_velocityFESpace->qr(),
1414  M_fespaceUETA,
1415  M_fespaceUETA,
1416  dot( M_density * phi_j * grad(M_fespaceUETA, uk_rep), phi_i )
1417  )
1418  >> M_Jacobian;
1419  }
1420  }
1421  M_Jacobian->globalAssemble();
1422 
1423  *M_block00 += *M_Jacobian;
1424  if ( !M_useStabilization )
1425  M_block00->globalAssemble( );
1426 }
1427 
1428 void NavierStokesSolverBlocks::updateStabilization( const vector_Type& convective_velocity_previous_newton_step,
1429  const vector_Type& velocity_previous_newton_step,
1430  const vector_Type& pressure_previous_newton_step,
1431  const vector_Type& velocity_rhs )
1432 {
1433  M_displayer.leaderPrint ( "[F] - Update Jacobian stabilization terms\n" ) ;
1434 
1435  if ( M_useFastAssembly )
1436  {
1437  vector_Type beta_km1( convective_velocity_previous_newton_step, Repeated);
1438  vector_Type u_km1( velocity_previous_newton_step, Repeated);
1439  vector_Type p_km1( pressure_previous_newton_step, Repeated);
1440  vector_Type u_bdf( velocity_rhs, Repeated);
1441 
1442  M_block01->zero();
1443  M_block10->zero();
1444  M_block11->zero();
1445 
1446  matrixPtr_Type systemMatrix_fast_block_00 (new matrix_Type ( M_velocityFESpace->map() ) );
1447  matrixPtr_Type systemMatrix_fast_block_01 (new matrix_Type ( M_velocityFESpace->map() ) );
1448  matrixPtr_Type systemMatrix_fast_block_10 (new matrix_Type ( M_pressureFESpace->map() ) );
1449  matrixPtr_Type systemMatrix_fast_block_11 (new matrix_Type ( M_pressureFESpace->map() ) );
1450  systemMatrix_fast_block_00->zero();
1451  systemMatrix_fast_block_01->zero();
1452  systemMatrix_fast_block_10->zero();
1453  systemMatrix_fast_block_11->zero();
1454 
1455  M_fastAssembler->supg_FI_FSI_terms( systemMatrix_fast_block_00, systemMatrix_fast_block_01, systemMatrix_fast_block_10, systemMatrix_fast_block_11,
1456  beta_km1, u_km1, p_km1, u_bdf);
1457 
1458  systemMatrix_fast_block_00->globalAssemble();
1459  systemMatrix_fast_block_01->globalAssemble( M_pressureFESpace->mapPtr(), M_velocityFESpace->mapPtr() );
1460  systemMatrix_fast_block_10->globalAssemble( M_velocityFESpace->mapPtr(), M_pressureFESpace->mapPtr() );
1461  systemMatrix_fast_block_11->globalAssemble();
1462 
1463  *M_block00 += *systemMatrix_fast_block_00;
1464  M_block00->globalAssemble();
1465 
1466  M_block01->zero();
1467  *M_block01 += *M_Btranspose;
1468  *M_block01 += *systemMatrix_fast_block_01;
1469  M_block01->globalAssemble( M_pressureFESpace->mapPtr(), M_velocityFESpace->mapPtr() );
1470 
1471  M_block10->zero();
1472  *M_block10 += *M_B;
1473  *M_block10 += *systemMatrix_fast_block_10;
1474  M_block10->globalAssemble(M_velocityFESpace->mapPtr(), M_pressureFESpace->mapPtr());
1475 
1476  M_block11->zero();
1477  *M_block11 += *systemMatrix_fast_block_11;
1478  M_block11->globalAssemble();
1479  }
1480  else
1481  {
1482  M_stabilization->apply_matrix(convective_velocity_previous_newton_step,
1483  velocity_previous_newton_step,
1484  pressure_previous_newton_step,
1485  velocity_rhs);
1486 
1487  *M_block00 += *M_stabilization->block_00();
1488  M_block00->globalAssemble();
1489 
1490  M_block01->zero();
1491  *M_block01 += *M_Btranspose;
1492  *M_block01 += *M_stabilization->block_01();
1493  M_block01->globalAssemble( M_pressureFESpace->mapPtr(), M_velocityFESpace->mapPtr() );
1494 
1495  M_block10->zero();
1496  *M_block10 += *M_B;
1497  *M_block10 += *M_stabilization->block_10();
1498  M_block10->globalAssemble(M_velocityFESpace->mapPtr(), M_pressureFESpace->mapPtr());
1499 
1500  M_block11->zero();
1501  *M_block11 += *M_stabilization->block_11();
1502  M_block11->globalAssemble();
1503  }
1504 
1505 
1506 }
1507 
1508 void NavierStokesSolverBlocks::solveJac( vector_Type& increment, const vector_Type& residual, const Real /*linearRelTol*/ )
1509 {
1510  // Apply BCs on the jacobian matrix
1511  applyBoundaryConditionsJacobian ( M_bc_sol );
1512 
1513  //(1) Set up the OseenOperator
1514  M_displayer.leaderPrint( "\tNS operator - set up the block operator...");
1515  LifeChrono chrono;
1516  chrono.start();
1517 
1518  Operators::NavierStokesOperator::operatorPtrContainer_Type operDataJacobian ( 2, 2 );
1519  operDataJacobian ( 0, 0 ) = M_block00->matrixPtr();
1520  operDataJacobian ( 0, 1 ) = M_block01->matrixPtr();
1521  operDataJacobian ( 1, 0 ) = M_block10->matrixPtr();
1522  if ( M_useStabilization )
1523  operDataJacobian ( 1, 1 ) = M_block11->matrixPtr();
1524 
1525  M_oper->setUp(operDataJacobian, M_displayer.comm());
1526  chrono.stop();
1527  M_displayer.leaderPrintMax(" done in " , chrono.diff() );
1528 
1529  M_displayer.leaderPrint( "\tPreconditioner operator - set up the block operator...");
1530  chrono.reset();
1531  chrono.start();
1532 
1533  //(2) Set the data for the preconditioner
1534  if ( std::strcmp(M_prec->Label(),"aSIMPLEOperator")==0 )
1535  {
1536  if (M_useStabilization)
1537  {
1538  M_prec->setUp(M_block00, M_block10, M_block01, M_block11);
1539  }
1540  else
1541  {
1542  M_prec->setUp(M_block00, M_block10, M_block01);
1543  }
1544 
1545  M_prec->setDomainMap(M_oper->OperatorDomainBlockMapPtr());
1546  M_prec->setRangeMap(M_oper->OperatorRangeBlockMapPtr());
1547  M_prec->updateApproximatedMomentumOperator();
1548  M_prec->updateApproximatedSchurComplementOperator();
1549  }
1550 
1551  chrono.stop();
1552  M_displayer.leaderPrintMax(" done in " , chrono.diff() );
1553 
1554  //(3) Set the solver for the linear system
1555  //(3) Set the solver for the linear system
1556  M_displayer.leaderPrint( "\tset up the Trilinos solver...");
1557  chrono.start();
1558  std::string solverType(M_pListLinSolver->get<std::string>("Linear Solver Type"));
1559  M_invOper.reset(Operators::InvertibleOperatorFactory::instance().createObject(solverType));
1560  M_invOper->setParameterList(M_pListLinSolver->sublist(solverType));
1561 
1562  chrono.stop();
1563  M_displayer.leaderPrintMax(" done in " , chrono.diff() );
1564 
1565  M_invOper->setOperator(M_oper);
1566  M_invOper->setPreconditioner(M_prec);
1567 
1568  increment.zero();
1569 
1570  // Solving the linear system
1571  M_invOper->ApplyInverse(residual.epetraVector(), increment.epetraVector());
1572 
1573  vector_Type increment_velocity ( M_velocityFESpace->map(), Unique ) ;
1574  vector_Type increment_pressure ( M_pressureFESpace->map(), Unique ) ;
1575 
1576  increment_velocity.zero();
1577  increment_pressure.zero();
1578 
1579  increment_velocity.subset(increment);
1580  increment_pressure.subset(increment, M_velocityFESpace->dof().numTotalDof() * 3);
1581 
1582 }
1583 
1584 void NavierStokesSolverBlocks::applyBoundaryConditionsJacobian ( bcPtr_Type & bc )
1585 {
1586  bcManageMatrix( *M_block00, *M_velocityFESpace->mesh(), M_velocityFESpace->dof(), *bc, M_velocityFESpace->feBd(), 1.0, 0.0);
1587  bcManageMatrix( *M_block01, *M_velocityFESpace->mesh(), M_velocityFESpace->dof(), *bc, M_velocityFESpace->feBd(), 0.0, 0.0);
1588  M_block01->globalAssemble(M_pressureFESpace->mapPtr(), M_velocityFESpace->mapPtr());
1589 }
1590 
1591 void NavierStokesSolverBlocks::preprocessBoundary(const Real& nx, const Real& ny, const Real& nz, BCHandler& bc, Real& Q_hat, const vectorPtr_Type& Phi_h, const UInt flag,
1592  vectorPtr_Type& Phi_h_flag, vectorPtr_Type& V_hat_x, vectorPtr_Type& V_hat_y, vectorPtr_Type& V_hat_z)
1593 {
1594  Phi_h_flag.reset ( new vector_Type ( M_velocityFESpaceScalar->map() ) );
1595  Phi_h_flag->zero();
1596 
1597  bcManageRhs ( *Phi_h_flag, *M_velocityFESpaceScalar->mesh(), M_velocityFESpaceScalar->dof(), bc, M_velocityFESpaceScalar->feBd(), 1., 0.);
1598 
1599  *Phi_h_flag *= *Phi_h; // Element by element multiplication
1600 
1601  Q_hat = 0.0;
1602 
1603  // Computing the flowrate associated to Phi_h_inflow
1604 
1605  vectorPtr_Type Phi_h_flag_rep;
1606  Phi_h_flag_rep.reset ( new vector_Type ( *Phi_h_flag, Repeated ) );
1607 
1608  vectorPtr_Type Q_hat_vec;
1609  Q_hat_vec.reset ( new vector_Type ( M_velocityFESpaceScalar->map() ) );
1610  Q_hat_vec->zero();
1611 
1612  vectorPtr_Type ones_vec;
1613  ones_vec.reset ( new vector_Type ( M_velocityFESpaceScalar->map() ) );
1614  ones_vec->zero();
1615  *ones_vec += 1.0;
1616 
1617  {
1618  using namespace ExpressionAssembly;
1619  QuadratureBoundary myBDQR (buildTetraBDQR (quadRuleTria6pt) );
1620  integrate (
1621  boundary (M_fespaceUETA_scalar->mesh(), flag),
1622  myBDQR,
1623  M_fespaceUETA_scalar,
1624  value(M_fespaceUETA_scalar, *Phi_h_flag_rep)*phi_i
1625  )
1626  >> Q_hat_vec;
1627  }
1628 
1629  Q_hat = Q_hat_vec->dot(*ones_vec);
1630 
1631  V_hat_x.reset ( new vector_Type ( *Phi_h_flag ) );
1632  V_hat_y.reset ( new vector_Type ( *Phi_h_flag ) );
1633  V_hat_z.reset ( new vector_Type ( *Phi_h_flag ) );
1634 
1635  *V_hat_x *= nx;
1636  *V_hat_y *= ny;
1637  *V_hat_z *= nz;
1638 }
1639 
1640 void NavierStokesSolverBlocks::solveLaplacian( const UInt& /*flag*/, bcPtr_Type& bc_laplacian, vectorPtr_Type& laplacianSolution)
1641 {
1642  // Update BCs for the laplacian problem
1643  bc_laplacian->bcUpdate ( *M_velocityFESpaceScalar->mesh(), M_velocityFESpaceScalar->feBd(), M_velocityFESpaceScalar->dof() );
1644 
1645  vectorPtr_Type Phi_h;
1646  Phi_h.reset ( new vector_Type ( M_velocityFESpaceScalar->map() ) );
1647  Phi_h->zero();
1648 
1649  vectorPtr_Type rhs_laplacian_repeated( new vector_Type (M_velocityFESpaceScalar->map(), Repeated ) );
1650  rhs_laplacian_repeated->zero();
1651 
1652  std::shared_ptr< MatrixEpetra<Real> > Laplacian;
1653  Laplacian.reset ( new MatrixEpetra<Real> ( M_velocityFESpaceScalar->map() ) );
1654  Laplacian->zero();
1655 
1656  {
1657  using namespace ExpressionAssembly;
1658 
1659  integrate(
1660  elements(M_fespaceUETA_scalar->mesh()),
1661  M_velocityFESpaceScalar->qr(),
1662  M_fespaceUETA_scalar,
1663  M_fespaceUETA_scalar,
1664  dot( grad(phi_i) , grad(phi_j) )
1665  ) >> Laplacian;
1666  }
1667 
1668  {
1669  using namespace ExpressionAssembly;
1670 
1671  integrate(
1672  elements(M_fespaceUETA_scalar->mesh()),
1673  M_velocityFESpaceScalar->qr(),
1674  M_fespaceUETA_scalar,
1675  value(1.0)*phi_i
1676  ) >> rhs_laplacian_repeated;
1677  }
1678 
1679 
1680  rhs_laplacian_repeated->globalAssemble();
1681 
1682  vectorPtr_Type rhs_laplacian( new vector_Type (*rhs_laplacian_repeated, Unique ) );
1683 
1684  bcManage ( *Laplacian, *rhs_laplacian, *M_velocityFESpaceScalar->mesh(), M_velocityFESpaceScalar->dof(), *bc_laplacian, M_velocityFESpaceScalar->feBd(), 1.0, 0.0 );
1685  Laplacian->globalAssemble();
1686 
1687  Teuchos::RCP< Teuchos::ParameterList > belosList = Teuchos::rcp ( new Teuchos::ParameterList );
1688  belosList = Teuchos::getParametersFromXmlFile ( "SolverParamListLaplacian.xml" );
1689 
1690  // Preconditioner
1691  prec_Type* precRawPtr;
1692  basePrecPtr_Type precPtr;
1693  precRawPtr = new prec_Type;
1694  precRawPtr->setDataFromGetPot ( M_dataFile, "prec" );
1695  precPtr.reset ( precRawPtr );
1696 
1697  // Linear Solver
1698  LinearSolver solver;
1699  solver.setCommunicator ( M_comm );
1700  solver.setParameters ( *belosList );
1701  solver.setPreconditioner ( precPtr );
1702 
1703  // Solve system
1704  solver.setOperator ( Laplacian );
1705  solver.setRightHandSide ( rhs_laplacian );
1706  solver.solve ( Phi_h );
1707 
1708  *laplacianSolution = *Phi_h;
1709 }
1710 
1711 void
1712 NavierStokesSolverBlocks::setupPostProc( )
1713 {
1714  M_postProcessing.reset ( new PostProcessingBoundary<mesh_Type> ( M_velocityFESpace->mesh(),
1715  &M_velocityFESpace->feBd(),
1716  &M_velocityFESpace->dof(),
1717  &M_pressureFESpace->feBd(),
1718  &M_pressureFESpace->dof(),
1719  *M_monolithicMap ) );
1720 
1721 }
1722 
1723 Real
1724 NavierStokesSolverBlocks::flux ( const markerID_Type& flag, const vector_Type& velocity )
1725 {
1726  vector_Type velocity_rep ( velocity, Repeated );
1727  return M_postProcessing->flux ( velocity_rep, flag );
1728 }
1729 
1730 Real
1731 NavierStokesSolverBlocks::area ( const markerID_Type& flag )
1732 {
1733  return M_postProcessing->measure ( flag );
1734 }
1735 
1736 Vector
1737 NavierStokesSolverBlocks::geometricCenter ( const markerID_Type& flag )
1738 {
1739  return M_postProcessing->geometricCenter ( flag );
1740 }
1741 
1742 Vector
1743 NavierStokesSolverBlocks::normal ( const markerID_Type& flag )
1744 {
1745  return M_postProcessing->normal ( flag );
1746 }
1747 
1748 Real
1749 NavierStokesSolverBlocks::pres ( const markerID_Type& flag, const vector_Type& pressure )
1750 {
1751  vector_Type pressure_rep ( pressure, Repeated );
1752  return M_postProcessing->average ( pressure_rep, flag, 1 ) [0];
1753 }
1754 
1755 void
1756 NavierStokesSolverBlocks::setBoundaryConditions ( const bcPtr_Type& bc )
1757 {
1758  M_bc_sol.reset ( new BCHandler ( ) );
1759  M_bc_res_essential.reset ( new BCHandler ( ) );
1760  M_bc_res_natural.reset ( new BCHandler ( ) );
1761 
1762  for ( std::vector<BCBase>::iterator it = bc->begin(); it != bc->end(); it++ )
1763  {
1764  if ( it->type() > 3 ) // meaning esssential
1765  {
1766  M_bc_sol->addBC( *it );
1767  M_bc_res_essential->addBC( *it );
1768  }
1769  else if ( it->type() == 0 ) // meaning natural
1770  {
1771  M_bc_res_natural->addBC( *it );
1772  }
1773  }
1774 }
1775 
1776 }
LifeV::SignFunction::SignFunction
SignFunction(const SignFunction &)
Definition: NavierStokesSolverBlocks.cpp:25
LifeV::SignFunction::comm
std::shared_ptr< Epetra_Comm > comm
Definition: NavierStokesSolverBlocks.cpp:28
LifeV::SignFunction::operator()
return_Type operator()(const Real &a)
Definition: NavierStokesSolverBlocks.cpp:12
LifeV::SignFunction::SignFunction
SignFunction(const std::shared_ptr< Epetra_Comm > &communicator)
Definition: NavierStokesSolverBlocks.cpp:24
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175