FSIMonolithic.cpp 26.9 KB
Newer Older
1
2
//@HEADER
/*
3
*******************************************************************************
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
22
23
    You should have received a copy of the GNU Lesser General Public License
    along with LifeV.  If not, see <http://www.gnu.org/licenses/>.

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


paolo.crosetto's avatar
paolo.crosetto committed
28
29
#include <EpetraExt_Reindex_MultiVector.h>

30

31
#include <lifev/core/LifeV.hpp>
Radu Popescu's avatar
Radu Popescu committed
32
#include <lifev/fsi/solver/FSIMonolithic.hpp>
33

34
35
namespace LifeV
{
36
37

// ===================================================
38
// Constructors, Destructor
39
// ===================================================
Alessio Fumagalli's avatar
Alessio Fumagalli committed
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
FSIMonolithic::FSIMonolithic() :
    super_Type(),
    M_monolithicMap(),
    M_interfaceMap(),
    M_beta(),
    M_monolithicMatrix(),
    M_precPtr(),
    M_rhsFull(),
    M_BCh_flux(),
    M_BChWS(),
    M_offset (0),
    M_solidAndFluidDim (0),
    M_fluidBlock(),
    M_solidBlockPrec(),
    M_linearSolver(),
    M_numerationInterface(),
    M_BChs(),
    M_FESpaces(),
    M_diagonalScale (false),
    M_reusePrec (true),
    M_resetPrec (true),
    M_maxIterSolver (-1),
    M_restarts (false),
    //end of protected attributes
    M_preconditionedSymmetrizedMatrix(),
    M_stress(),
    M_fluxes (0)
67
68
69
70
71
72
73
{}

FSIMonolithic::~FSIMonolithic()
{
}

// ===================================================
74
// Setup Methods
75
76
77
78
79
80
81
// ===================================================
void
FSIMonolithic::setupFEspace()
{
    super_Type::setupFEspace();

    // Monolitic: In the beginning I need a non-partitioned mesh. later we will do the partitioning
Alessio Fumagalli's avatar
Alessio Fumagalli committed
82
83
84
85
    M_dFESpace.reset ( new FESpace<mesh_Type, MapEpetra> ( M_solidMesh,
                                                           M_data->dataSolid()->order(),
                                                           nDimensions,
                                                           M_epetraComm) );
paolo.crosetto's avatar
paolo.crosetto committed
86

87
88
89
}

void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
90
FSIMonolithic::setupDOF ( void )
91
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
92
93
    M_dofStructureToHarmonicExtension    .reset ( new DOFInterface3Dto3D );
    M_dofStructureToFluid    .reset ( new DOFInterface3Dto3D );
94

Alessio Fumagalli's avatar
Alessio Fumagalli committed
95
96
97
98
99
100
    M_dofStructureToHarmonicExtension->setup (   M_mmFESpace->refFE(), M_mmFESpace->dof(),
                                                 M_dFESpace->refFE(), M_dFESpace->dof() );
    M_dofStructureToHarmonicExtension->update ( *M_mmFESpace->mesh(),  M_data->fluidInterfaceFlag(),
                                                *M_dFESpace->mesh(),  M_data->structureInterfaceFlag(),
                                                M_data->interfaceTolerance(),
                                                M_data->fluidInterfaceVertexFlag() );
101

Alessio Fumagalli's avatar
Alessio Fumagalli committed
102
103
104
105
106
107
    M_dofStructureToFluid->setup (   M_uFESpace->refFE(), M_uFESpace->dof(),
                                     M_dFESpace->refFE(), M_dFESpace->dof() );
    M_dofStructureToFluid->update ( *M_uFESpace->mesh(),  M_data->fluidInterfaceFlag(),
                                    *M_dFESpace->mesh(),  M_data->structureInterfaceFlag(),
                                    M_data->interfaceTolerance(),
                                    M_data->fluidInterfaceVertexFlag() );
108

Alessio Fumagalli's avatar
Alessio Fumagalli committed
109
    createInterfaceMaps (M_dofStructureToFluid/*HarmonicExtension*/->localDofMap() );
110
111
112

    M_fluidMesh.reset();
    M_solidMesh.reset();
113
114
}

115
#ifdef HAVE_HDF5
116
void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
117
FSIMonolithic::setupDOF ( meshFilter_Type& filterMesh )
118
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
119
    createInterfaceMaps (*filterMesh.getStoredInterface (0) );
120
}
121
#endif
122
123
124
125

void
FSIMonolithic::setupSystem( )
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
126
127
    M_fluid->setUp ( M_dataFile );
    setup ( M_dataFile );
128
129
130
}

void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
131
FSIMonolithic::setup ( const GetPot& dataFile )
132
133
{

Alessio Fumagalli's avatar
Alessio Fumagalli committed
134
    M_linearSolver.reset (new solver_Type (M_epetraComm) );
135

Alessio Fumagalli's avatar
Alessio Fumagalli committed
136
137
138
    M_linearSolver->setDataFromGetPot ( dataFile, "linear_system/solver" );
    std::string prectype = dataFile ("problem/DDBlockPrec", "PRECTYPE_UNDEFINED");
    std::string opertype = dataFile ("problem/blockOper", "PRECTYPE_UNDEFINED");
139

Alessio Fumagalli's avatar
Alessio Fumagalli committed
140
    M_precPtr.reset (BlockPrecFactory::instance().createObject ( prectype ) );
141

Alessio Fumagalli's avatar
Alessio Fumagalli committed
142
143
144
145
146
147
148
149
150
151
    M_precPtr->setupSolver (*M_linearSolver, dataFile);
    std::string section ("linear_system/prec");
    M_precPtr->setComm (M_epetraComm);
    M_precPtr->setDataFromGetPot (dataFile, section); //to avoid if we build the prec from a matrix.
    M_monolithicMatrix->setComm (M_epetraComm);
    M_monolithicMatrix->setDataFromGetPot (dataFile, section); //to avoid if we build the prec from a matrix.
    M_reusePrec     = dataFile ( "linear_system/prec/reuse", true);
    M_maxIterSolver = dataFile ( "linear_system/solver/max_iter", -1);
    M_diagonalScale    = dataFile ( "linear_system/prec/diagonalScaling",  false );
    M_restarts         = dataFile ( "exporter/start"  ,  0   );
152
153
154
155
156
157
}

void
FSIMonolithic::setupFluidSolid( )
{
    M_BCh_flux = M_BCh_u; // For the moment M_BCh_u contains only the fluxes.
Alessio Fumagalli's avatar
Alessio Fumagalli committed
158
    M_fluxes = M_BCh_u->size( );
159

Alessio Fumagalli's avatar
Alessio Fumagalli committed
160
    setupFluidSolid ( M_fluxes );
161

Alessio Fumagalli's avatar
Alessio Fumagalli committed
162
    M_BCh_flux->setOffset (M_offset - M_fluxes);
163
164
    std::vector<BCBase>::iterator fluxIt = M_BCh_flux->begin( );
    for ( UInt i = 0; i < M_fluxes; ++i, ++fluxIt )
Alessio Fumagalli's avatar
Alessio Fumagalli committed
165
166
167
    {
        fluxIt->setOffset ( i );
    }
168

169
}
170
171

void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
172
FSIMonolithic::setupFluidSolid ( UInt const fluxes )
173
174
175
{

    // Note: up to now it works only with matching grids (and poly order) on the interface
Alessio Fumagalli's avatar
Alessio Fumagalli committed
176
    assert (M_fluidInterfaceMap->map (Unique)->NumGlobalElements() == M_solidInterfaceMap->map (Unique)->NumGlobalElements() );
177
178
179
180
181
182

    M_interfaceMap = M_solidInterfaceMap;

    //std::map<ID, ID> const& localDofMap = M_dofStructureToHarmonicExtension->localDofMap();
    std::map<ID, ID>::const_iterator ITrow;

Alessio Fumagalli's avatar
Alessio Fumagalli committed
183
    M_monolithicMap.reset (new MapEpetra (M_uFESpace->map() ) );
184

Alessio Fumagalli's avatar
Alessio Fumagalli committed
185
    std::string opertype = M_dataFile ("problem/blockOper", "AdditiveSchwarz");
186

Alessio Fumagalli's avatar
Alessio Fumagalli committed
187
    createOperator ( opertype );
188

Alessio Fumagalli's avatar
Alessio Fumagalli committed
189
190
191
    *M_monolithicMap += M_pFESpace->map();
    *M_monolithicMap += fluxes;
    *M_monolithicMap += M_dFESpace->map();
192

Alessio Fumagalli's avatar
Alessio Fumagalli committed
193
    M_monolithicMatrix->createInterfaceMap ( *M_interfaceMap, M_dofStructureToFluid->localDofMap(), M_dFESpace->map().map (Unique)->NumGlobalElements() / nDimensions, M_epetraWorldComm );
194
195
196
    *M_monolithicMap += *M_monolithicMatrix->interfaceMap();

    //the map for the interface coupling matrices should be done with respect to the coarser mesh.
Alessio Fumagalli's avatar
Alessio Fumagalli committed
197
198
    M_monolithicMatrix->numerationInterface (M_numerationInterface);
    M_beta.reset  (new vector_Type (M_uFESpace->map() ) );
199

Alessio Fumagalli's avatar
Alessio Fumagalli committed
200
201
202
    M_offset = M_uFESpace->dof().numTotalDof() * nDimensions + fluxes +  M_pFESpace->dof().numTotalDof();
    M_solidAndFluidDim = M_offset + M_dFESpace->dof().numTotalDof() * nDimensions;
    M_BCh_d->setOffset (M_offset);
203
204
205
}

// ===================================================
206
// Public Methods
207
208
// ===================================================
void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
209
FSIMonolithic::monolithicToInterface (vector_Type& lambdaSolid, const vector_Type& disp)
210
211
212
{
    if (disp.mapType() == Repeated)
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
213
214
        vector_Type const  dispUnique (disp, Unique);
        monolithicToInterface (lambdaSolid, dispUnique);
215
216
217
218
        return;
    }
    if (lambdaSolid.mapType() == Repeated)
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
219
220
        vector_Type  lambdaSolidUn (lambdaSolid.map(), Unique);
        monolithicToInterface ( lambdaSolidUn, disp);
221
222
223
224
        lambdaSolid = lambdaSolidUn;
        return;
    }

Alessio Fumagalli's avatar
Alessio Fumagalli committed
225
226
227
228
    MapEpetra subMap (*disp.map().map (Unique), M_offset, disp.map().map (Unique)->NumGlobalElements() );
    vector_Type subDisp (subMap, Unique);
    subDisp.subset (disp, M_offset);
    lambdaSolid = subDisp;
229
230
231
}

void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
232
FSIMonolithic::monolithicToX (const vector_Type& disp, vector_Type& dispFluid, MapEpetra& map, UInt offset)
233
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
234
    if (disp.mapType() == Repeated)
235
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
236
237
        vector_Type dispUnique (disp, Unique);
        monolithicToX (dispUnique, dispFluid, map, offset);
238
239
240
        dispFluid = dispUnique;
        return;
    }
Alessio Fumagalli's avatar
Alessio Fumagalli committed
241
    dispFluid.subset (disp, map, offset, offset);
242
243
244
245
246
247
}


void
FSIMonolithic::buildSystem()
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
248
    M_solid->buildSystem ( M_solidTimeAdvance->coefficientSecondDerivative ( 0 ) / (M_data->dataSolid()->dataTime()->timeStep() *M_data->dataSolid()->dataTime()->timeStep() ) );
249
250
251
252
253
254
255
256
}

#ifdef HAVE_TRILINOS_ANASAZI
Real&
FSIMonolithic::computeMaxSingularValue( )
{
    typedef Epetra_Operator                                                operatorPtr_Type;

Alessio Fumagalli's avatar
Alessio Fumagalli committed
257
    M_preconditionedSymmetrizedMatrix.reset (new ComposedOperator<Epetra_Operator> (M_epetraComm) );
258

Alessio Fumagalli's avatar
Alessio Fumagalli committed
259
    boost::shared_ptr<operatorPtr_Type>  ComposedPrecPtr (M_linearSolver->preconditioner()->preconditioner() );
260
261

    //M_monolithicMatrix->getMatrixPtr()->OptimizeStorage();
Alessio Fumagalli's avatar
Alessio Fumagalli committed
262
263
264
265
266
    boost::shared_ptr<Epetra_FECrsMatrix> matrCrsPtr (new Epetra_FECrsMatrix (*M_monolithicMatrix->matrix()->matrixPtr() ) );
    M_preconditionedSymmetrizedMatrix->push_back (boost::static_pointer_cast<operatorPtr_Type> (/*ComposedPrecPtr*/matrCrsPtr) );
    M_preconditionedSymmetrizedMatrix->push_back ( (ComposedPrecPtr/*matrCrsPtr*/), true);
    M_preconditionedSymmetrizedMatrix->push_back ( (ComposedPrecPtr/*matrCrsPtr*/), true, true);
    M_preconditionedSymmetrizedMatrix->push_back (boost::static_pointer_cast<operatorPtr_Type> (/*ComposedPrecPtr*/matrCrsPtr), false, true);
267
268
269
270
271
272

    std::vector<LifeV::Real> real;
    std::vector<LifeV::Real> imaginary;

    boost::shared_ptr<EigenSolver> eig;

Alessio Fumagalli's avatar
Alessio Fumagalli committed
273
274
    UInt nev = M_dataFile ("eigensolver/nevec", 10); //number of eigenvectors
    if (nev)
275
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
276
277
        eig.reset (new EigenSolver (M_preconditionedSymmetrizedMatrix, M_preconditionedSymmetrizedMatrix->OperatorDomainMap(), nev) );
        eig->setDataFromGetPot (M_dataFile, "eigensolver/");
278
        eig->solve();
Alessio Fumagalli's avatar
Alessio Fumagalli committed
279
        eig->eigenvalues (real, imaginary);
280
281
282
283
284
    }
    else
    {
        throw UNDEF_EIGENSOLVER_EXCEPTION();
    }
Alessio Fumagalli's avatar
Alessio Fumagalli committed
285
    for (UInt i = 0; i < real.size(); ++i)
286
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
287
288
        displayer().leaderPrint ("\n real part ", real[i]);
        displayer().leaderPrint ("\n imaginary part ", imaginary[i]);
289
290
291
292
293
294
    }
    return real[0];
}
#endif

void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
295
FSIMonolithic::computeFluidNormals ( vector_Type& normals)
296
297
{
    BCManageNormal<matrix_Type> normalManager;
Alessio Fumagalli's avatar
Alessio Fumagalli committed
298
299
300
301
302
303
    if ( !M_BChWS->bcUpdateDone() ) //possibly to avoid
    {
        M_BChWS->bcUpdate (*M_uFESpace->mesh(), M_uFESpace->feBd(), M_uFESpace->dof() );
    }
    normalManager.init ( (*M_BChWS) [0], 0.);
    normalManager.computeIntegratedNormals (M_uFESpace->dof(), M_uFESpace->feBd(), normals, *M_uFESpace->mesh() );
304
305
306
}

void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
307
FSIMonolithic::solveJac ( vector_Type& step, const vector_Type& res, const Real /*linearRelTol*/ )
308
309
310
{
    setupBlockPrec( );

Alessio Fumagalli's avatar
Alessio Fumagalli committed
311
    checkIfChangedFluxBC ( M_precPtr );
paolo.crosetto's avatar
paolo.crosetto committed
312

313
    M_precPtr->blockAssembling();
Alessio Fumagalli's avatar
Alessio Fumagalli committed
314
    M_precPtr->applyBoundaryConditions ( dataFluid()->dataTime()->time() );
315
316
    M_precPtr->GlobalAssemble();

317
#ifdef HAVE_LIFEV_DEBUG
Alessio Fumagalli's avatar
Alessio Fumagalli committed
318
    M_solid->displayer().leaderPrint ("  M-  Residual NormInf:                        ", res.normInf(), "\n");
319
#endif
Alessio Fumagalli's avatar
Alessio Fumagalli committed
320
    iterateMonolithic (res, step);
321
#ifdef HAVE_LIFEV_DEBUG
Alessio Fumagalli's avatar
Alessio Fumagalli committed
322
    M_solid->displayer().leaderPrint ("  M-  Solution NormInf:                        ", step.normInf(), "\n");
323
#endif
324
325
326
327
328
}

void
FSIMonolithic::updateSystem()
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
329
330
    *M_rhs *= 0;
    *M_rhsFull *= 0;
331
332
333
    this->M_fluid->resetStabilization();
}

334

335
// ===================================================
336
// Protected Methods
337
338
// ===================================================
void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
339
FSIMonolithic::iterateMonolithic (const vector_Type& rhs, vector_Type& step)
340
{
Radu Popescu's avatar
Radu Popescu committed
341
    LifeChrono chrono;
342

Alessio Fumagalli's avatar
Alessio Fumagalli committed
343
    displayer().leaderPrint ("  M-  Solving the system ... \n" );
344
345
346

    //M_monolithicMatrix->GlobalAssemble();
    //necessary if we did not imposed Dirichlet b.c.
347

Alessio Fumagalli's avatar
Alessio Fumagalli committed
348
    M_linearSolver->setOperator (*M_monolithicMatrix->matrix()->matrixPtr() );
349

Alessio Fumagalli's avatar
Alessio Fumagalli committed
350
    M_linearSolver->setReusePreconditioner ( (M_reusePrec) && (!M_resetPrec) );
351

Alessio Fumagalli's avatar
Alessio Fumagalli committed
352
    int numIter = M_precPtr->solveSystem ( rhs, step, M_linearSolver );
353
354
355
356
357

    if (numIter < 0)
    {
        chrono.start();

Alessio Fumagalli's avatar
Alessio Fumagalli committed
358
        M_solid->displayer().leaderPrint ("   Iterative solver failed, numiter = ", -numIter );
359
360

        if (numIter <= -M_maxIterSolver)
Alessio Fumagalli's avatar
Alessio Fumagalli committed
361
362
363
        {
            M_solid->displayer().leaderPrint ("   ERROR: Iterative solver failed.\n");
        }
364
365
    }

366
    //M_solid->getDisplayer().leaderPrint("  M-  System solved.\n" );
367
368
369
}

void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
370
FSIMonolithic::couplingRhs (vectorPtr_Type rhs) // not working with non-matching grids
371
{
372
    std::map<ID, ID> const& localDofMap = M_dofStructureToFluid->localDofMap();
373
374
    std::map<ID, ID>::const_iterator ITrow;

Alessio Fumagalli's avatar
Alessio Fumagalli committed
375
376
    vector_Type rhsStructureVelocity (M_solidTimeAdvance->rhsContributionFirstDerivative() *M_solid->rescaleFactor(), Unique);
    vector_Type lambda (*M_interfaceMap, Unique);
paolo.crosetto's avatar
paolo.crosetto committed
377

Alessio Fumagalli's avatar
Alessio Fumagalli committed
378
    this->monolithicToInterface (lambda, rhsStructureVelocity);
paolo.crosetto's avatar
paolo.crosetto committed
379

Alessio Fumagalli's avatar
Alessio Fumagalli committed
380
381
    UInt interface (M_monolithicMatrix->interface() );
    UInt totalDofs (M_dFESpace->dof().numTotalDof() );
382
383


Alessio Fumagalli's avatar
Alessio Fumagalli committed
384
    for (UInt dim = 0; dim < nDimensions; ++dim)
385
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
386
        for ( ITrow = localDofMap.begin(); ITrow != localDofMap.end() ; ++ITrow)
387
        {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
388
            if (M_interfaceMap->map (Unique)->LID (ITrow->second /*+ dim*solidDim*/) >= 0 ) //to avoid repeated stuff
389
            {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
390
391
392
393
                if (rhs.get() )
                {
                    (*rhs) [  (int) (*M_numerationInterface) [ITrow->second ] + dim * interface + M_solidAndFluidDim ] = -lambda ( ITrow->second + dim * totalDofs );
                }
394
395
396
397
398
399
400
            }
        }
    }
}

void
FSIMonolithic::
Alessio Fumagalli's avatar
Alessio Fumagalli committed
401
evalResidual ( const vector_Type& sol, const vectorPtr_Type& rhs, vector_Type& res, bool diagonalScaling)
402
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
403
    if ( diagonalScaling )
Cristiano Malossi's avatar
Cristiano Malossi committed
404
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
405
406
407
408
409
        diagonalScale (*rhs, M_monolithicMatrix->matrix() );
    }
    if (! (M_data->dataSolid()->solidType().compare ("exponential") && M_data->dataSolid()->solidType().compare ("neoHookean") ) )
    {
        M_solid->apply (sol * M_solid->rescaleFactor(), res);
Cristiano Malossi's avatar
Cristiano Malossi committed
410
        M_fluidBlock->globalAssemble();
411

Alessio Fumagalli's avatar
Alessio Fumagalli committed
412
        res += ( (*M_fluidBlock) * sol);
413

Alessio Fumagalli's avatar
Alessio Fumagalli committed
414
        res += *M_monolithicMatrix->coupling() * sol;
Cristiano Malossi's avatar
Cristiano Malossi committed
415
416
417
    }
    else
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
418
        res = * (M_monolithicMatrix->matrix() ) * sol;
Cristiano Malossi's avatar
Cristiano Malossi committed
419
420
        res -= *rhs; // Ax-b
    }
421
}
422
423
424

void
FSIMonolithic::
Alessio Fumagalli's avatar
Alessio Fumagalli committed
425
updateSolidSystem ( vectorPtr_Type& rhsFluidCoupling )
426
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
427
428
    Real coefficient ( M_data->dataSolid()->dataTime()->timeStep() * M_data->dataSolid()->dataTime()->timeStep() * M_solid->rescaleFactor() /  M_solidTimeAdvance->coefficientSecondDerivative ( 0 ) );
    M_solidTimeAdvance->updateRHSContribution ( M_data->dataSolid()->dataTime()->timeStep() );
429
    *rhsFluidCoupling += *M_solid->massMatrix() * ( M_solidTimeAdvance->rhsContributionSecondDerivative() * coefficient );
430
431
    // TODO NOTE: this mass * vector multiplication in serial may lead to a NaN for unclear reasons
    // (both the matrix and the vector does not contain a NaN before the multiplication..)
432
433
}

Alessio Fumagalli's avatar
Alessio Fumagalli committed
434
435
436
437
438
void FSIMonolithic::setVectorInStencils ( const vectorPtr_Type& vel,
                                          const vectorPtr_Type& pressure,
                                          const vectorPtr_Type& solidDisp,
                                          //const vectorPtr_Type& fluidDisp,
                                          const UInt iter)
439
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
440
441
442
    setFluidVectorInStencil (vel, pressure, iter);
    setSolidVectorInStencil (solidDisp, iter);
    //  setALEVectorInStencil(fluidDisp, iter);
443
444
445

}

grandper's avatar
grandper committed
446
447
448
void FSIMonolithic::setFluidVectorInStencil ( const vectorPtr_Type& vel,
                                              const vectorPtr_Type& pressure,
                                              const UInt /*iter*/ )
449
450
451
452
453
454
{

    //The fluid and solid TimeAdvance classes have a stencil of dimension
    //as big as the coupled problem.

    //Fluid Problem
Alessio Fumagalli's avatar
Alessio Fumagalli committed
455
456
    vectorPtr_Type vectorMonolithicFluidVelocity (new vector_Type (*M_monolithicMap, Unique, Zero) );
    vectorPtr_Type vectorMonolithicFluidPressure (new vector_Type (*M_monolithicMap, Unique, Zero) );
457
458
459
460

    *vectorMonolithicFluidVelocity *= 0.0;
    *vectorMonolithicFluidPressure *= 0.0;

Alessio Fumagalli's avatar
Alessio Fumagalli committed
461
462
    vectorMonolithicFluidVelocity->subset (*vel, vel->map(), UInt (0), UInt (0) ) ;
    vectorMonolithicFluidPressure->subset ( *pressure, pressure->map(), UInt (0), (UInt) 3 * M_uFESpace->dof().numTotalDof() );
463
464
465

    *vectorMonolithicFluidVelocity += *vectorMonolithicFluidPressure;

Alessio Fumagalli's avatar
Alessio Fumagalli committed
466
467
    vector_Type* normalPointerToFluidVector ( new vector_Type (*vectorMonolithicFluidVelocity) );
    (M_fluidTimeAdvance->stencil() ).push_back ( normalPointerToFluidVector );
468
469
470
}


Alessio Fumagalli's avatar
Alessio Fumagalli committed
471
472
void FSIMonolithic::setSolidVectorInStencil ( const vectorPtr_Type& solidDisp,
                                              const UInt iter)
473
474
{
    //Solid problem
Alessio Fumagalli's avatar
Alessio Fumagalli committed
475
476
477
    vectorPtr_Type vectorMonolithicSolidDisplacement (new vector_Type (*M_monolithicMap, Unique, Zero) );
    *vectorMonolithicSolidDisplacement *= 0.0;
    vectorMonolithicSolidDisplacement->subset ( *solidDisp, solidDisp->map(), (UInt) 0, M_offset);
478
479
    *vectorMonolithicSolidDisplacement *= 1.0 / M_solid->rescaleFactor();

480
481
482
    //The fluid timeAdvance has size = orderBDF because it is seen as an equation frist order in time.
    //We need to add the solidVector to the fluidVector in the fluid TimeAdvance because we have the
    //extrapolation on it.
Alessio Fumagalli's avatar
Alessio Fumagalli committed
483
484
485
486
    if ( iter <= M_fluidTimeAdvance->size() - 1 )
    {
        * ( M_fluidTimeAdvance->stencil() [ iter ] ) += *vectorMonolithicSolidDisplacement;
    }
487

Alessio Fumagalli's avatar
Alessio Fumagalli committed
488
489
    vector_Type* normalPointerToSolidVector ( new vector_Type (*vectorMonolithicSolidDisplacement) );
    (M_solidTimeAdvance->stencil() ).push_back ( normalPointerToSolidVector );
490
491
492

}

Paolo Tricerri's avatar
Paolo Tricerri committed
493
void FSIMonolithic::finalizeRestart( )
494
495
{
    //Set the initialRHS for the TimeAdvance classes
Alessio Fumagalli's avatar
Alessio Fumagalli committed
496
497
    vector_Type zeroFluidSolid (*M_monolithicMap, Unique, Zero);
    vector_Type zeroALE (M_mmFESpace->map(), Unique, Zero);
498
499
500
501

    zeroFluidSolid *= 0.0;
    zeroALE *= 0.0;

Alessio Fumagalli's avatar
Alessio Fumagalli committed
502
503
504
    M_fluidTimeAdvance->setInitialRHS (zeroFluidSolid);
    M_solidTimeAdvance->setInitialRHS (zeroFluidSolid);
    M_ALETimeAdvance->setInitialRHS (zeroALE);
505
506
507
508
509

    //This updates at the current value (the one when the simulation was stopped) the RHScontribution
    //of the first derivative which is use to compute the velocity in TimeAdvance::velocity().
    //Please note that, even if it is ugly, at this stage, the fluidTimeAdvance is leading the Time Discretization
    //and this is why there  is the dataFluid class to get the dt.
Alessio Fumagalli's avatar
Alessio Fumagalli committed
510
    M_ALETimeAdvance->updateRHSFirstDerivative ( M_data->dataFluid()->dataTime()->timeStep() );
511
512
513
}


514
515
void
FSIMonolithic::
Alessio Fumagalli's avatar
Alessio Fumagalli committed
516
diagonalScale (vector_Type& rhs, matrixPtr_Type matrFull)
517
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
518
    Epetra_Vector diagonal (*rhs.map().map (Unique) );
519
520
521
    //M_matrFull->matrixPtr()->InvRowSums(diagonal);
    //M_matrFull->matrixPtr()->InvRowMaxs(diagonal);
    //M_matrFull->matrixPtr()->InvColSums(diagonal);
Alessio Fumagalli's avatar
Alessio Fumagalli committed
522
523
524
    matrFull->matrixPtr()->InvColMaxs (diagonal);
    matrFull->matrixPtr()->LeftScale (diagonal);
    rhs.epetraVector().Multiply (1, rhs.epetraVector(), diagonal, 0);
525
526
527
}

void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
528
FSIMonolithic::variablesInit (const std::string& dOrder)
529
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
530
531
532
533
    M_dFESpace.reset (new FESpace<mesh_Type, MapEpetra> (M_solidLocalMesh,
                                                         dOrder,
                                                         3,
                                                         M_epetraComm) );
534
    // INITIALIZATION OF THE VARIABLES
Alessio Fumagalli's avatar
Alessio Fumagalli committed
535
536
    M_lambdaFluid.reset (new vector_Type (*M_fluidInterfaceMap, Unique) );
    M_lambdaFluidRepeated.reset (new vector_Type (*M_fluidInterfaceMap, Repeated) );
537
538
}

539
void FSIMonolithic::setupBlockPrec( )
540
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
541
    if (! (M_precPtr->set() ) )
Cristiano Malossi's avatar
Cristiano Malossi committed
542
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
543
544
545
546
547
548
        M_precPtr->push_back_matrix (M_solidBlockPrec, M_structureNonLinear);
        M_precPtr->push_back_matrix (M_fluidBlock, true);
        M_precPtr->setConditions (M_BChs);
        M_precPtr->setSpaces (M_FESpaces);
        M_precPtr->setOffsets (2, M_offset, 0);
        M_precPtr->coupler (M_monolithicMap, M_dofStructureToFluid->localDofMap(), M_numerationInterface, M_data->dataFluid()->dataTime()->timeStep(), M_solidTimeAdvance->coefficientFirstDerivative ( 0 ), M_solid->rescaleFactor() );
Cristiano Malossi's avatar
Cristiano Malossi committed
549
    }
paolo.crosetto's avatar
paolo.crosetto committed
550
551
    else
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
552
553
        M_precPtr->replace_matrix (M_fluidBlock, 1);
        M_precPtr->replace_matrix (M_solidBlockPrec, 0);
paolo.crosetto's avatar
paolo.crosetto committed
554
    }
555
556

#ifdef HAVE_NS_PREC
Cristiano Malossi's avatar
Cristiano Malossi committed
557
    /* This shall be enabled once the branch about NS_PREC is going to be merged to master*/
558
    boost::shared_ptr<MonolithicBlockComposed> blockPrecPointer ( boost::dynamic_pointer_cast<MonolithicBlockComposed> M_precPtr);
559

560
561
    if (blockPrecPointer.get() != 0)
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
562
563
564
        UInt fluidPosition = blockPrecPointer->whereIsBlock (MonolithicBlockComposed::fluid);
        ASSERT (blockPrecPointer->blockPrecs().size() < fluidPosition, "The preconditioner corresponding to the fluid block is probably not PCD. Check in the data file");
        boost::shared_ptr<PreconditionerPCD> prec_PCD ( boost::dynamic_pointer_cast<PreconditionerPCD> blockPrecPointer->blockPrecs() [fluidPosition] );
565
566


Cristiano Malossi's avatar
Cristiano Malossi committed
567
568
        if (prec_PCD.get() != 0)
        {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
569
570
571
572
573
574
575
576
            prec_PCD->setFESpace (M_uFESpace, M_pFESpace);
            prec_PCD->setBCHandler (M_BCh_u);
            prec_PCD->setTimestep (M_data->dataFluid()->dataTime()->timeStep() );
            prec_PCD->setViscosity (M_data->dataFluid()->viscosity() );
            prec_PCD->setDensity (M_data->dataFluid()->density() );
            prec_PCD->setCouplingMatrixView (M_precPtr->couplingVector() [MonolithicBlockComposed::fluid]);
            prec_PCD->setMapStructure (&M_dFESpace->map() );
            prec_PCD->updateBeta (M_fluidTimeAdvance->singleElement (0) );
Cristiano Malossi's avatar
Cristiano Malossi committed
577
        }
578
    }
579
#endif
580
581
582
}

void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
583
FSIMonolithic::assembleSolidBlock ( UInt iter, const vector_Type& solution )
584
585
{
    if (iter == 0)
Alessio Fumagalli's avatar
Alessio Fumagalli committed
586
587
588
    {
        updateSolidSystem (this->M_rhs);
    }
589

590

Alessio Fumagalli's avatar
Alessio Fumagalli committed
591
    if (M_data->dataSolid()->solidType().compare ("exponential") && M_data->dataSolid()->solidType().compare ("neoHookean") )
Cristiano Malossi's avatar
Cristiano Malossi committed
592
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
593
594
        M_solid->material()->computeStiffness (solution * M_solid->rescaleFactor(), 1., M_data->dataSolid(), M_solid->mapMarkersVolumes(), M_solid->displayerPtr() );
        M_solidBlockPrec.reset (new matrix_Type (*M_monolithicMap, 1) );
595
        *M_solidBlockPrec += *M_solid->massMatrix();
Cristiano Malossi's avatar
Cristiano Malossi committed
596
597
598
599
600
601
        *M_solidBlockPrec += *M_solid->material()->stiffMatrix();
        M_solidBlockPrec->globalAssemble();
        *M_solidBlockPrec *= M_solid->rescaleFactor();
    }
    else
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
602
603
        M_solid->material()->updateJacobianMatrix ( solution * M_solid->rescaleFactor(), dataSolid(), M_solid->mapMarkersVolumes(), M_solid->displayerPtr() ); // computing the derivatives if nonlinear (comment this for inexact Newton);
        M_solidBlockPrec.reset (new matrix_Type (*M_monolithicMap, 1) );
604
        *M_solidBlockPrec += *M_solid->massMatrix();
Cristiano Malossi's avatar
Cristiano Malossi committed
605
606
607
608
        *M_solidBlockPrec += *M_solid->material()->jacobian(); //stiffMatrix();
        M_solidBlockPrec->globalAssemble();
        *M_solidBlockPrec *= M_solid->rescaleFactor();
    }
paolo.crosetto's avatar
paolo.crosetto committed
609

610
611
612
}

void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
613
FSIMonolithic::assembleFluidBlock (UInt iter, const vector_Type& solution)
614
{
Antonio Cervone's avatar
Antonio Cervone committed
615
    M_fluidBlock.reset (new  FSIOperator::fluid_Type::matrix_Type (*M_monolithicMap) );
616

Alessio Fumagalli's avatar
Alessio Fumagalli committed
617
    Real alpha = M_fluidTimeAdvance->coefficientFirstDerivative ( 0 ) / M_data->dataFluid()->dataTime()->timeStep(); //mesh velocity w
Paolo Tricerri's avatar
Paolo Tricerri committed
618
619

    //This line is based on the hypothesis that the conservativeFormulation flag is set on FALSE
Alessio Fumagalli's avatar
Alessio Fumagalli committed
620
    M_fluid->updateSystem (alpha, *this->M_beta, *this->M_rhs, M_fluidBlock, solution );
Paolo Tricerri's avatar
Paolo Tricerri committed
621
622

    //This in the case of conservativeFormulation == true
623
624
    // else
    //   if (! M_fluid->matrixMassPtr().get() )
Alessio Fumagalli's avatar
Alessio Fumagalli committed
625
    //  M_fluid->buildSystem( );
626

Alessio Fumagalli's avatar
Alessio Fumagalli committed
627
    if (iter == 0)
Cristiano Malossi's avatar
Cristiano Malossi committed
628
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
629
630
631
632
633
634
        M_resetPrec = true;
        M_fluidTimeAdvance->updateRHSContribution ( M_data->dataFluid()->dataTime()->timeStep() );
        if (!M_data->dataFluid()->conservativeFormulation() )
        {
            *M_rhs += M_fluid->matrixMass() * (M_fluidTimeAdvance->rhsContributionFirstDerivative() );
        }
Cristiano Malossi's avatar
Cristiano Malossi committed
635
        else
Alessio Fumagalli's avatar
Alessio Fumagalli committed
636
637
638
639
        {
            *M_rhs += (M_fluidMassTimeAdvance->rhsContributionFirstDerivative() );
        }
        couplingRhs (M_rhs/*, M_fluidTimeAdvance->stencil()[0]*/);
Cristiano Malossi's avatar
Cristiano Malossi committed
640
    }
641
    //the conservative formulation as it is now is of order 1. To have higher order (TODO) we need to store the mass matrices computed at the previous time steps.
Paolo Tricerri's avatar
Paolo Tricerri committed
642
    //At the moment, the flag conservativeFormulation should be always kept on FALSE
Alessio Fumagalli's avatar
Alessio Fumagalli committed
643
644
645
646
647
    if (M_data->dataFluid()->conservativeFormulation() )
    {
        M_fluid->updateSystem (alpha, *this->M_beta, *this->M_rhs, M_fluidBlock, solution );
    }
    this->M_fluid->updateStabilization (*M_fluidBlock);
648
649
}

650
void FSIMonolithic::updateRHS()
651
{
652
    // Update fluid (iter == 0)
Alessio Fumagalli's avatar
Alessio Fumagalli committed
653
654
655
    M_fluidTimeAdvance->updateRHSContribution ( M_data->dataFluid()->dataTime()->timeStep() );
    *M_rhs += M_fluid->matrixMass() * (M_fluidTimeAdvance->rhsContributionFirstDerivative() );
    couplingRhs (M_rhs);
656

657
    // Update solid (iter == 0)
Alessio Fumagalli's avatar
Alessio Fumagalli committed
658
    updateSolidSystem (M_rhs);
659
660

    // Update RHS
661
662
663
    *M_rhsFull = *M_rhs;
}

paolo.crosetto's avatar
paolo.crosetto committed
664
665
namespace
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
666
667
668
669
670
671
672
673
static Real fZero (const Real& /*t*/, const Real& /*x*/, const Real& /*y*/, const Real& /*z*/, const ID& /*i*/)
{
    return 0.;
}
static Real fOne (const Real& /*t*/, const Real& /*x*/, const Real& /*y*/, const Real& /*z*/, const ID& /*i*/)
{
    return 1.;
}
paolo.crosetto's avatar
paolo.crosetto committed
674
675
}

Alessio Fumagalli's avatar
Alessio Fumagalli committed
676
void FSIMonolithic::enableStressComputation (UInt flag)
paolo.crosetto's avatar
paolo.crosetto committed
677
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
678
679
680
681
    M_BChWS.reset (new solidBchandler_Type() );
    BCFunctionBase bcfZero (fZero);
    BCFunctionBase bcfOne (fOne);
    M_bcfWs.setFunctions_Robin (bcfOne, bcfOne);
paolo.crosetto's avatar
paolo.crosetto committed
682

Alessio Fumagalli's avatar
Alessio Fumagalli committed
683
    M_BChWS->addBC ("WS", (UInt) flag, Robin, Full, M_bcfWs, 3);
paolo.crosetto's avatar
paolo.crosetto committed
684
685
686
687
}

FSIMonolithic::vectorPtr_Type FSIMonolithic::computeStress()
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
688
689
    vector_Type lambda (M_monolithicMatrix->interfaceMap() );
    lambda.subset (M_fluidTimeAdvance->singleElement (0), M_solidAndFluidDim);
paolo.crosetto's avatar
paolo.crosetto committed
690

Alessio Fumagalli's avatar
Alessio Fumagalli committed
691
    M_boundaryMass.reset (new matrix_Type (*M_interfaceMap) );
paolo.crosetto's avatar
paolo.crosetto committed
692
    if ( !M_BChWS->bcUpdateDone() )
Alessio Fumagalli's avatar
Alessio Fumagalli committed
693
694
695
696
    {
        M_BChWS->bcUpdate (*M_dFESpace->mesh(), M_dFESpace->feBd(), M_dFESpace->dof() );
    }
    bcManageMatrix (*M_boundaryMass, *M_dFESpace->mesh(), M_dFESpace->dof(), *M_BChWS, M_dFESpace->feBd(), 1., dataSolid()->dataTime()->time() );
697
    M_boundaryMass->globalAssemble();
paolo.crosetto's avatar
paolo.crosetto committed
698

Alessio Fumagalli's avatar
Alessio Fumagalli committed
699
700
701
    solver_Type solverMass (M_epetraComm);
    solverMass.setDataFromGetPot ( M_dataFile, "solid/solver" );
    solverMass.setupPreconditioner (M_dataFile, "solid/prec"); //to avoid if we have already a prec.
paolo.crosetto's avatar
paolo.crosetto committed
702

Alessio Fumagalli's avatar
Alessio Fumagalli committed
703
    boost::shared_ptr<PreconditionerIfpack> P (new PreconditionerIfpack() );
paolo.crosetto's avatar
paolo.crosetto committed
704

Alessio Fumagalli's avatar
Alessio Fumagalli committed
705
706
707
708
    vectorPtr_Type sol (new vector_Type (M_monolithicMatrix->interfaceMap() ) );
    solverMass.setMatrix (*M_boundaryMass);
    solverMass.setReusePreconditioner (false);
    solverMass.solveSystem ( lambda, *sol, M_boundaryMass);
paolo.crosetto's avatar
paolo.crosetto committed
709

Alessio Fumagalli's avatar
Alessio Fumagalli committed
710
711
712
    EpetraExt::MultiVector_Reindex reindexMV (*M_interfaceMap->map (Unique) );
    boost::shared_ptr<MapEpetra> newMap (new MapEpetra ( *M_interfaceMap ) );
    M_stress.reset (new vector_Type (reindexMV (lambda.epetraVector() ), newMap, Unique) );
paolo.crosetto's avatar
paolo.crosetto committed
713
714
715
716
    return M_stress;
}

void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
717
FSIMonolithic::checkIfChangedFluxBC ( precPtr_Type oper )
paolo.crosetto's avatar
paolo.crosetto committed
718
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
719
720
    UInt nfluxes (M_BChs[1]->numberOfBCWithType (Flux) );
    if (M_fluxes != nfluxes)
paolo.crosetto's avatar
paolo.crosetto committed
721
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
722
        for (UInt i = 0; i < M_fluxes; ++i)
paolo.crosetto's avatar
paolo.crosetto committed
723
        {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
724
725
            const BCBase* bc (M_BChs[1]->findBCWithName (M_BCFluxNames[i]) );
            if (bc->type() != Flux)
paolo.crosetto's avatar
paolo.crosetto committed
726
            {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
727
                oper->addToCoupling (1., M_fluxOffset[i], M_fluxOffset[i], 1 );
paolo.crosetto's avatar
paolo.crosetto committed
728
729
730
731
732
733
            }
        }
    }
}


734
}