FSIMonolithicGI.hpp 9.21 KB
 paolo.crosetto committed Dec 13, 2010 1 2 3 //@HEADER /* *******************************************************************************  crosetto committed Aug 03, 2010 4   paolo.crosetto committed Dec 13, 2010 5 6  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University  crosetto committed Aug 03, 2010 7   paolo.crosetto committed Dec 13, 2010 8  This file is part of LifeV.  crosetto committed Aug 03, 2010 9   paolo.crosetto committed Dec 13, 2010 10 11 12 13  LifeV is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.  crosetto committed Aug 03, 2010 14   paolo.crosetto committed Dec 13, 2010 15 16 17 18  LifeV is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.  crosetto committed Aug 03, 2010 19   paolo.crosetto committed Dec 13, 2010 20 21  You should have received a copy of the GNU Lesser General Public License along with LifeV. If not, see .  crosetto committed Aug 03, 2010 22   paolo.crosetto committed Dec 13, 2010 23 *******************************************************************************  crosetto committed Aug 03, 2010 24 */  paolo.crosetto committed Dec 13, 2010 25 //@HEADER  crosetto committed Aug 03, 2010 26 /**  paolo.crosetto committed Dec 13, 2010 27   Cristiano Malossi committed Oct 25, 2012 28 29 30 31 32 33 34 35 36 /*! * @file * @brief File containing the Monolithic Geometry--Implicit FSI Solver * * @date 18-09-2008 * @author Paolo Crosetto * * @maintainer Paolo Crosetto */  paolo.crosetto committed Dec 13, 2010 37   crosetto committed Aug 03, 2010 38 39 40 #ifndef _MONOLITHICGI_HPP #define _MONOLITHICGI_HPP  Radu Popescu committed Mar 06, 2012 41 #include  crosetto committed Aug 03, 2010 42   Alessio Fumagalli committed Feb 12, 2013 43 44 namespace LifeV {  crosetto committed Aug 03, 2010 45 46 47 48 #ifdef OBSOLETE class Epetra_FullMonolithic; #endif  Cristiano Malossi committed Oct 25, 2012 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 typedef FactorySingleton< Factory< FSIOperator, std::string > > FSIFactory_Type; //! FSIMonolithicGI Geometry-Implicit solver /*! * @author Paolo Crosetto * Class handling the nonlinear monolithic solver for FSI problems. The (exact or inexact) * Newton algorithm is used to solve the nonlinearity. * The block structure of the jacobian matrix is * \f$\left(\begin{array}{ccc} C&B&S\\ D&N&0\\ 0&E&H \end{array}\right)\f$ * where \f$N\f$ represents the solid block, \f$C\f$ the fluid block, \f$H\f$ is the harmonic extension block, * while the extra * diagonal blocks represent the coupling. The implementation of the stress continuity coupling condition * is obtained by means of an augmented formulation. * Different possible preconditioners are implemented. * * Important parameters to set properly in the data file: * - useShapeDerivatives: if true the shape derivatives block is added to the Jacobian matrix; * - domainVelImplicit: if true the domain velocity w in the convective term is considered an unknown (at the time n+1); * - convectiveTermDer: false if the convective term is linearized (\f$u^{n+1}\nabla(u^n-w^n)\f$), * otherwise it can be either true (if we use the Newton method to solve the convective term nonlinearity) or false * (fixed-point method); * - semiImplicit: if true only one iteration of the nonlinear solver is performed. Otherwise * the nonlinear iterations continue up to the specified tolerance. Set it to true for the GCE; * - method: can be either monolithicGE, monolithicGI if the geometry is treated respectively explicitly or implicitly, * or exactJacobians, fixedPoint for partitioned strategies; * - blockOper: specifies the matrix type to be used for the linear system: if AdditiveSchwarz, the matrix is the standard * ine for GE; if AdditiveSchwarzRN the coupling blocks are of Robin type instead of Dirichlet and Neumann. The parameters * for the Robin coupling are alphaf and alphas in the data file. NOTE: this method has currently been tested only for * alphas=0. * - DDBlockPrec: specifies the possible preconditioners to use. Can be: AdditiveSchwarz, MonolithicBlockComposedDN, MonolithicBlockComposedDN2, * MonolithicBlockComposedNN, MonolithicBlockComposedDNND.  crosetto committed Aug 03, 2010 80  */  paolo.crosetto committed Dec 13, 2010 81   Cristiano Malossi committed Oct 25, 2012 82 83 84 class FSIMonolithicGI : public FSIMonolithic { public:  85   Cristiano Malossi committed Oct 25, 2012 86 87 88  typedef FSIMonolithic super_Type; typedef Preconditioner prec_Type; typedef boost::shared_ptr< prec_Type > prec_type;  89   Cristiano Malossi committed Oct 25, 2012 90 91  //!@name Constructor and Destructor //@{  crosetto committed Aug 03, 2010 92   Cristiano Malossi committed Oct 25, 2012 93 94  //! Empty Constructor FSIMonolithicGI();  95   Cristiano Malossi committed Oct 25, 2012 96 97  //! Destructor virtual ~FSIMonolithicGI() {}  98   Cristiano Malossi committed Oct 25, 2012 99  //@}  crosetto committed Aug 03, 2010 100   Cristiano Malossi committed Oct 25, 2012 101 102  //!@name Public Methods //@{  103   Cristiano Malossi committed Oct 25, 2012 104 105 106  /** Sets the parameters read from data file */  Alessio Fumagalli committed Feb 12, 2013 107  void setup ( const GetPot& dataFile );  crosetto committed Aug 03, 2010 108   Cristiano Malossi committed Oct 25, 2012 109  //! initializes the fluid and mesh problems, creates the map of the global matrix  Alessio Fumagalli committed Feb 12, 2013 110  void setupFluidSolid ( UInt const fluxes );  crosetto committed Aug 03, 2010 111   Cristiano Malossi committed Oct 25, 2012 112 113  //! builds the constant part of the fluid-structure-mesh motion matrix void buildSystem();  crosetto committed Aug 03, 2010 114   Cristiano Malossi committed Oct 25, 2012 115 116 117 118 119 120  /** evaluates the residual b-Ax \param res: output \param _sol: fluid domain displacement solution \param iter: current NonLinearRichardson (block Gauss Seidel for the tangent system) iteration */  Alessio Fumagalli committed Feb 12, 2013 121  void evalResidual ( vector_Type& res, const vector_Type& sol, const UInt iter );  paolo.crosetto committed Oct 19, 2010 122   Cristiano Malossi committed Oct 25, 2012 123 124 125 126 127 128 129  //!Apply the boundary conditions to each block composing the monolithic problem /** Sets the vectors of: boundary conditions, FESpaces, couplings, offsets, and sets the blocks in the composed operator which constitutes the monolithic problem. Then calls the applyBoundaryConditions of the MonolithicBlockMatrix operator, passing also the right hand side. */ void applyBoundaryConditions();  Paolo Tricerri committed Nov 09, 2012 130   Alessio Fumagalli committed Feb 12, 2013 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145  void updateSolution ( const vector_Type& solution ) { super_Type::updateSolution ( solution ); //The size of the vectors for the ALE is = dimension of the ALE problem //To do the shift right we first need to extract the fluid displacement //And then push it into the ALE timeAdvance class. vectorPtr_Type displacementToSave ( new vector_Type (M_mmFESpace->map() ) ); UInt offset ( M_solidAndFluidDim + nDimensions * M_interface ); displacementToSave->subset (solution, offset); //This updateRHSFirstDerivative has to be done before the shiftRight //In fact it updates the right hand side of the velocity using the //previous times. The method velocity() uses it and then, the compuation //of the velocity is done using the current time and the previous times.  Cristiano Malossi committed Feb 26, 2013 146  //M_ALETimeAdvance->updateRHSFirstDerivative ( M_data->dataFluid()->dataTime()->timeStep() );  Alessio Fumagalli committed Feb 12, 2013 147 148  M_ALETimeAdvance->shiftRight ( *displacementToSave ); }  Paolo Tricerri committed Oct 14, 2012 149   Paolo Tricerri committed Nov 09, 2012 150 151 152 153  //! Set vectors for restart /*! * Set vectors for restart */  Alessio Fumagalli committed Feb 12, 2013 154 155 156  void setALEVectorInStencil (const vectorPtr_Type& fluidDisp, const UInt iter, const bool lastVector);  Paolo Tricerri committed Nov 09, 2012 157   158   Cristiano Malossi committed Oct 25, 2012 159 160  //!@name Get Methods //@{  161   Cristiano Malossi committed Oct 25, 2012 162  //! get the current solution vector.  Alessio Fumagalli committed Feb 12, 2013 163  LIFEV_DEPRECATED ( const vector_Type& solution() const )  Cristiano Malossi committed Oct 25, 2012 164 165 166 167  { if ( M_epetraWorldComm->MyPID() == 0 ) { std::cerr << std::endl << "Warning: FSIMonolithic::solution() is deprecated!" << std::endl  Alessio Fumagalli committed Feb 12, 2013 168  << " You should not access the solution inside FSIOperator or FSIMonolithic!" << std::endl;  Cristiano Malossi committed Oct 25, 2012 169  }  crosetto committed Aug 03, 2010 170   Alessio Fumagalli committed Feb 12, 2013 171  return M_fluidTimeAdvance->singleElement (0);  Cristiano Malossi committed Oct 25, 2012 172  }  crosetto committed Aug 03, 2010 173   Cristiano Malossi committed Oct 25, 2012 174  //! getter for the map of fluid-structure-interface (without the mesh motion)  Alessio Fumagalli committed Feb 12, 2013 175 176 177 178  const MapEpetra& mapWithoutMesh() const { return *M_mapWithoutMesh; }  crosetto committed Aug 03, 2010 179   Cristiano Malossi committed Oct 25, 2012 180  //! getter for the global matrix of the system  Alessio Fumagalli committed Feb 12, 2013 181 182 183 184  const matrixPtr_Type matrixPtr() const { return M_monolithicMatrix->matrix(); }  185   Cristiano Malossi committed Oct 25, 2012 186  static bool S_register;  187   Cristiano Malossi committed Oct 25, 2012 188  //@}  paolo.crosetto committed Sep 23, 2011 189   Cristiano Malossi committed Oct 25, 2012 190 protected:  crosetto committed Aug 03, 2010 191   Cristiano Malossi committed Oct 25, 2012 192 193  //!@name Protected Methods //@{  paolo.crosetto committed Dec 13, 2010 194   Cristiano Malossi committed Oct 25, 2012 195 196  //! set the block preconditioner void setupBlockPrec();  197   Cristiano Malossi committed Oct 25, 2012 198  //@}  paolo.crosetto committed Mar 07, 2011 199   Cristiano Malossi committed Oct 25, 2012 200 private:  crosetto committed Aug 03, 2010 201   Cristiano Malossi committed Oct 25, 2012 202 203  //! @name Private Methods //@{  crosetto committed Aug 03, 2010 204   Cristiano Malossi committed Oct 25, 2012 205  //! Factory method for the system matrix, of type MonolithicBlockBase  Alessio Fumagalli committed Feb 12, 2013 206  void createOperator ( std::string& operType )  Cristiano Malossi committed Oct 25, 2012 207  {  Alessio Fumagalli committed Feb 12, 2013 208 209  M_monolithicMatrix.reset ( MonolithicBlockMatrix::Factory_Type::instance().createObject ( operType ) ); M_monolithicMatrix.reset ( MonolithicBlockMatrix::Factory_Type::instance().createObject ( operType ) );  Cristiano Malossi committed Oct 25, 2012 210  }  paolo.crosetto committed Dec 13, 2010 211   Cristiano Malossi committed Oct 25, 2012 212 213 214 215 216  /** calculates the terms due to the shape derivatives given the mesh increment deltaDisp. The shape derivative block is assembled in a matrix (not in a right hand side representing the matrix-vector multiplication) \param sdMatrix: output. Shape derivatives block to be summed to the Jacobian matrix. */  Antonio Cervone committed Sep 05, 2013 217  void shapeDerivatives ( FSIOperator::fluid_Type::matrixPtr_Type sdMatrix );  crosetto committed Aug 03, 2010 218   Cristiano Malossi committed Oct 25, 2012 219 220 221 222 223  //! assembles the mesh motion matrix. /*!In Particular it diagonalize the part of the matrix corresponding to the Dirichlet condition expressing the coupling \param iter: current iteration: used as flag to distinguish the first nonlinear iteration from the others */  Alessio Fumagalli committed Feb 12, 2013 224  void assembleMeshBlock ( UInt iter );  crosetto committed Aug 03, 2010 225   Cristiano Malossi committed Oct 25, 2012 226  //@}  paolo.crosetto committed Dec 13, 2010 227   228   Cristiano Malossi committed Oct 25, 2012 229 230  //!@name Private Members //@{  231   Paolo Tricerri committed Oct 28, 2012 232 233 234 235 236 237  boost::shared_ptr M_mapWithoutMesh; //This vector is used in the shapeDerivatives method since a //copy of the solution at the current iteration k is necessary vectorPtr_Type M_uk; UInt M_interface; matrixPtr_Type M_meshBlock;  Antonio Cervone committed Sep 05, 2013 238  FSIOperator::fluid_Type::matrixPtr_Type M_shapeDerivativesBlock;  Paolo Tricerri committed Oct 28, 2012 239 240  matrixPtr_Type M_solidDerBlock; //std::vector M_BChsLin;  Cristiano Malossi committed Oct 25, 2012 241  //@}  paolo.crosetto committed Aug 03, 2011 242   Cristiano Malossi committed Oct 25, 2012 243 244 245 246 247  //! Factory method static FSIOperator* instantiate() { return new FSIMonolithicGI(); }  paolo.crosetto committed Aug 03, 2011 248   Cristiano Malossi committed Oct 25, 2012 249 };  crosetto committed Aug 03, 2010 250   Cristiano Malossi committed Dec 05, 2012 251 252 253 254 255 256 //! Factory create function inline FSIMonolithic* createFSIMonolithicGI() { return new FSIMonolithicGI(); }  crosetto committed Aug 03, 2010 257 258 } #endif