FSIMonolithicGI.hpp 9.21 KB
Newer Older
1
2
3
//@HEADER
/*
*******************************************************************************
4

5
6
    Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
    Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7

8
    This file is part of LifeV.
9

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.
14

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.
19

20
21
    You should have received a copy of the GNU Lesser General Public License
    along with LifeV.  If not, see <http://www.gnu.org/licenses/>.
22

23
*******************************************************************************
24
*/
25
//@HEADER
26
/**
27

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 <paolo.crosetto@epfl.ch>
 *
 *  @maintainer Paolo Crosetto <paolo.crosetto@epfl.ch>
 */
37

38
39
40
#ifndef _MONOLITHICGI_HPP
#define _MONOLITHICGI_HPP

Radu Popescu's avatar
Radu Popescu committed
41
#include <lifev/fsi/solver/FSIMonolithic.hpp>
42

Alessio Fumagalli's avatar
Alessio Fumagalli committed
43
44
namespace LifeV
{
45
46
47
48
#ifdef OBSOLETE
class Epetra_FullMonolithic;
#endif

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 <paolo.crosetto@epfl.ch>
 *  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.
80
 */
81

82
83
84
class FSIMonolithicGI : public FSIMonolithic
{
public:
85

86
87
88
    typedef FSIMonolithic super_Type;
    typedef Preconditioner prec_Type;
    typedef boost::shared_ptr< prec_Type > prec_type;
89

90
91
    //!@name Constructor and Destructor
    //@{
92

93
94
    //! Empty Constructor
    FSIMonolithicGI();
95

96
97
    //! Destructor
    virtual ~FSIMonolithicGI() {}
98

99
    //@}
100

101
102
    //!@name Public Methods
    //@{
103

104
105
106
    /**
     Sets the parameters read from data file
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
107
    void setup ( const GetPot& dataFile );
108

109
    //! initializes the fluid and mesh problems, creates the map of the global matrix
Alessio Fumagalli's avatar
Alessio Fumagalli committed
110
    void setupFluidSolid ( UInt const fluxes );
111

112
113
    //! builds the constant part of the fluid-structure-mesh motion matrix
    void buildSystem();
114

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's avatar
Alessio Fumagalli committed
121
    void evalResidual ( vector_Type& res, const vector_Type& sol, const UInt iter );
122

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's avatar
Paolo Tricerri committed
130

Alessio Fumagalli's avatar
Alessio Fumagalli committed
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.
146
        //M_ALETimeAdvance->updateRHSFirstDerivative ( M_data->dataFluid()->dataTime()->timeStep() );
Alessio Fumagalli's avatar
Alessio Fumagalli committed
147
148
        M_ALETimeAdvance->shiftRight ( *displacementToSave );
    }
149

Paolo Tricerri's avatar
Paolo Tricerri committed
150
151
152
153
    //! Set vectors for restart
    /*!
     *  Set vectors for restart
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
154
155
156
    void setALEVectorInStencil (const vectorPtr_Type& fluidDisp,
                                const UInt iter,
                                const bool lastVector);
Paolo Tricerri's avatar
Paolo Tricerri committed
157

158

159
160
    //!@name Get Methods
    //@{
161

162
    //! get the current solution vector.
Alessio Fumagalli's avatar
Alessio Fumagalli committed
163
    LIFEV_DEPRECATED ( const vector_Type& solution() const )
164
165
166
167
    {
        if ( M_epetraWorldComm->MyPID() == 0 )
        {
            std::cerr << std::endl << "Warning: FSIMonolithic::solution() is deprecated!" << std::endl
Alessio Fumagalli's avatar
Alessio Fumagalli committed
168
                      << "         You should not access the solution inside FSIOperator or FSIMonolithic!" << std::endl;
169
        }
170

Alessio Fumagalli's avatar
Alessio Fumagalli committed
171
        return  M_fluidTimeAdvance->singleElement (0);
172
    }
173

174
    //! getter for the map of fluid-structure-interface (without the mesh motion)
Alessio Fumagalli's avatar
Alessio Fumagalli committed
175
176
177
178
    const MapEpetra& mapWithoutMesh() const
    {
        return *M_mapWithoutMesh;
    }
179

180
    //! getter for the global matrix of the system
Alessio Fumagalli's avatar
Alessio Fumagalli committed
181
182
183
184
    const matrixPtr_Type matrixPtr() const
    {
        return M_monolithicMatrix->matrix();
    }
185

186
    static bool S_register;
187

188
    //@}
189

190
protected:
191

192
193
    //!@name Protected Methods
    //@{
194

195
196
    //! set the block preconditioner
    void setupBlockPrec();
197

198
    //@}
paolo.crosetto's avatar
paolo.crosetto committed
199

200
private:
201

202
203
    //! @name Private Methods
    //@{
204

205
    //! Factory method for the system matrix, of type MonolithicBlockBase
Alessio Fumagalli's avatar
Alessio Fumagalli committed
206
    void createOperator ( std::string& operType )
207
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
208
209
        M_monolithicMatrix.reset ( MonolithicBlockMatrix::Factory_Type::instance().createObject ( operType ) );
        M_monolithicMatrix.reset ( MonolithicBlockMatrix::Factory_Type::instance().createObject ( operType ) );
210
    }
211

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's avatar
Antonio Cervone committed
217
    void shapeDerivatives ( FSIOperator::fluid_Type::matrixPtr_Type sdMatrix );
218

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's avatar
Alessio Fumagalli committed
224
    void assembleMeshBlock ( UInt iter );
225

226
    //@}
227

228

229
230
    //!@name Private Members
    //@{
231

232
233
234
235
236
237
    boost::shared_ptr<MapEpetra>         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's avatar
Antonio Cervone committed
238
    FSIOperator::fluid_Type::matrixPtr_Type M_shapeDerivativesBlock;
239
240
    matrixPtr_Type                       M_solidDerBlock;
    //std::vector<fluidBchandlerPtr_Type>    M_BChsLin;
241
    //@}
242

243
244
245
246
247
    //! Factory method
    static FSIOperator* instantiate()
    {
        return new FSIMonolithicGI();
    }
248

249
};
250

251
252
253
254
255
256
//! Factory create function
inline FSIMonolithic* createFSIMonolithicGI()
{
    return new FSIMonolithicGI();
}

257
258
}
#endif