cylinder.cpp 15.5 KB
Newer Older
gfourestey's avatar
gfourestey committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
/* -*- mode: c++ -*-

  This file is part of the LifeV Applications.

  Author(s): Christophe Prud'homme <christophe.prudhomme@epfl.ch>
       Date: 2005-04-19

  Copyright (C) 2005 EPFL

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2.1 of the License, or
  (at your option) any later version.

  This program 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
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
  USA
*/
/**
   \file cylinder.cpp
   \author Christophe Prud'homme <christophe.prudhomme@epfl.ch>
   \date 2005-04-19
 */

31
32
33
// Tell the compiler to ignore specific kind of warnings:
#pragma GCC diagnostic ignored "-Wunused-variable"
#pragma GCC diagnostic ignored "-Wunused-parameter"
gfourestey's avatar
gfourestey committed
34
35
36

#include <Epetra_ConfigDefs.h>
#ifdef EPETRA_MPI
Tiziano Passerini's avatar
Tiziano Passerini committed
37
#include <mpi.h>
38
#include <Epetra_MpiComm.h>
gfourestey's avatar
gfourestey committed
39
#else
Tiziano Passerini's avatar
Tiziano Passerini committed
40
#include <Epetra_SerialComm.h>
gfourestey's avatar
gfourestey committed
41
#endif
42
43
44
45
46

//Tell the compiler to restore the warning previously silented
#pragma GCC diagnostic warning "-Wunused-variable"
#pragma GCC diagnostic warning "-Wunused-parameter"

gfourestey's avatar
gfourestey committed
47
//#include "life/lifesolver/NavierStokesSolver.hpp"
48
49
50
51
52
53
54
#include <lifev/core/array/MatrixEpetra.hpp>
#include <lifev/core/array/MapEpetra.hpp>
#include <lifev/core/mesh/MeshData.hpp>
#include <lifev/core/mesh/MeshPartitioner.hpp>
#include <lifev/navier_stokes/solver/OseenData.hpp>
#include <lifev/core/fem/FESpace.hpp>
#include <lifev/core/fem/TimeAdvanceBDFNavierStokes.hpp>
gfourestey's avatar
gfourestey committed
55
#ifdef HAVE_HDF5
56
#include <lifev/core/filter/ExporterHDF5.hpp>
gfourestey's avatar
gfourestey committed
57
#endif
58
#include <lifev/core/filter/ExporterEnsight.hpp>
gfourestey's avatar
gfourestey committed
59

60
#include <lifev/navier_stokes/solver/OseenSolver.hpp>
gfourestey's avatar
gfourestey committed
61
62
63
64
65
66
67
68

#include "cylinder.hpp"
#include <iostream>



using namespace LifeV;

69
typedef RegionMesh<LinearTetra>                             mesh_Type;
gfourestey's avatar
gfourestey committed
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

const int INLET       = 2;
const int WALL        = 1;
const int OUTLET      = 3;
const int RINGIN      = 20;
const int RINGOUT     = 30;


Real zero_scalar( const Real& /* t */,
                  const Real& /* x */,
                  const Real& /* y */,
                  const Real& /* z */,
                  const ID& /* i */ )
{
    return 0.;
}

Real u2(const Real& t, const Real& /*x*/, const Real& /*y*/, const Real& /*z*/, const ID& i)
{
Tiziano Passerini's avatar
Tiziano Passerini committed
89
90
    switch (i)
    {
uvilla's avatar
uvilla committed
91
    case 0:
Tiziano Passerini's avatar
Tiziano Passerini committed
92
93
        return 0.0;
        break;
uvilla's avatar
uvilla committed
94
    case 2:
Tiziano Passerini's avatar
Tiziano Passerini committed
95
96
        if ( t <= 0.003 )
            return 1.3332e4;
gfourestey's avatar
gfourestey committed
97
//      return 0.01;
Tiziano Passerini's avatar
Tiziano Passerini committed
98
99
        return 0.0;
        break;
uvilla's avatar
uvilla committed
100
    case 1:
Tiziano Passerini's avatar
Tiziano Passerini committed
101
        return 0.0;
gfourestey's avatar
gfourestey committed
102
103
104
//      return 1.3332e4;
//    else
//      return 0.0;
Tiziano Passerini's avatar
Tiziano Passerini committed
105
106
107
        break;
    }
    return 0;
gfourestey's avatar
gfourestey committed
108
109
110
}

void
Tiziano Passerini's avatar
Tiziano Passerini committed
111
postProcessFluxesPressures( OseenSolver< mesh_Type >& nssolver,
Tiziano Passerini's avatar
Tiziano Passerini committed
112
113
                            BCHandler& bcHandler,
                            const LifeV::Real& t, bool _verbose )
gfourestey's avatar
gfourestey committed
114
{
Tiziano Passerini's avatar
Tiziano Passerini committed
115
116
    LifeV::Real Q, P;
    UInt flag;
gfourestey's avatar
gfourestey committed
117

perego's avatar
perego committed
118
    for ( BCHandler::bcBaseIterator_Type it = bcHandler.begin();
Tiziano Passerini's avatar
Tiziano Passerini committed
119
            it != bcHandler.end(); ++it )
gfourestey's avatar
gfourestey committed
120
    {
Tiziano Passerini's avatar
Tiziano Passerini committed
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
        flag = it->flag();

        Q = nssolver.flux(flag);
        P = nssolver.pressure(flag);

        if ( _verbose )
        {
            std::ofstream outfile;
            std::stringstream filenamess;
            std::string filename;

            // file name contains the label
            filenamess << flag;
            // writing down fluxes
            filename = "flux_label" + filenamess.str() + ".m";
            outfile.open(filename.c_str(),std::ios::app);
            outfile << Q << " " << t << "\n";
            outfile.close();
            // writing down pressures
            filename = "pressure_label" + filenamess.str() + ".m";
            outfile.open(filename.c_str(),std::ios::app);
            outfile << P << " " << t << "\n";
            outfile.close();
            // reset ostringstream
            filenamess.str("");
        }
gfourestey's avatar
gfourestey committed
147
148
149
150
151
152
153
154
    }

}


struct Cylinder::Private
{
    Private() :
Tiziano Passerini's avatar
Tiziano Passerini committed
155
156
157
158
159
160
161
            //check(false),
            nu(1),
            //rho(1),
            H(1), D(1)
            //H(20), D(1)
            //H(0.41), D(0.1)
    {}
162
    typedef boost::function<Real ( Real const&, Real const&, Real const&, Real const&, ID const& )> fct_Type;
gfourestey's avatar
gfourestey committed
163
164
165
166
167
168
169
170
171
172
173
174
175

    double Re;

    std::string data_file_name;

    double      nu;  /**< viscosity (in m^2/s) */
    //const double rho; /**< density is constant (in kg/m^3) */
    double      H;   /**< height and width of the domain (in m) */
    double      D;   /**< diameter of the cylinder (in m) */
    bool        centered; /**< true if the cylinder is at the origin */

    std::string initial_sol;

176
    boost::shared_ptr<Epetra_Comm>   comm;
gfourestey's avatar
gfourestey committed
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
    /**
     * get the characteristic velocity
     *
     * @return the characteristic velocity
     */
    double Ubar() const { return nu*Re/D; }

    /**
     * get the magnitude of the profile velocity
     *
     *
     * @return the magnitude of the profile velocity
     */
    double Um_3d() const { return 9*Ubar()/4; }

    double Um_2d() const { return 3*Ubar()/2; }


    /**
     * u3d 3D velocity profile.
     *
     * Define the velocity profile at the inlet for the 3D cylinder
     */
    Real u3d( const Real& /* t */,
              const Real& /* x */,
              const Real& y,
              const Real& z,
              const ID&   id ) const
Tiziano Passerini's avatar
Tiziano Passerini committed
205
    {
uvilla's avatar
uvilla committed
206
        if ( id == 0 )
gfourestey's avatar
gfourestey committed
207
        {
Tiziano Passerini's avatar
Tiziano Passerini committed
208
209
            if ( centered )
            {
210
                return Um_3d() * (H+y)*(H-y) * (H+z)*(H-z) / std::pow(H,4);
Tiziano Passerini's avatar
Tiziano Passerini committed
211
212
213
            }
            else
            {
214
                return 16 * Um_3d() * y * z * (H-y) * (H-z) / std::pow(H,4);
gfourestey's avatar
gfourestey committed
215
216
            }
        }
Tiziano Passerini's avatar
Tiziano Passerini committed
217
        else
gfourestey's avatar
gfourestey committed
218
        {
Tiziano Passerini's avatar
Tiziano Passerini committed
219
            return 0;
gfourestey's avatar
gfourestey committed
220
        }
Tiziano Passerini's avatar
Tiziano Passerini committed
221
222
    }

223
    fct_Type getU_3d()
Tiziano Passerini's avatar
Tiziano Passerini committed
224
    {
225
        fct_Type f;
Tiziano Passerini's avatar
Tiziano Passerini committed
226
227
228
        f = boost::bind(&Cylinder::Private::u3d, this, _1, _2, _3, _4, _5);
        return f;
    }
gfourestey's avatar
gfourestey committed
229
230
231
232
233
234
235

    /**
     * u2d flat 2D velocity profile.
     *
     * Define the velocity profile at the inlet for the 2D cylinder
     */
    Real u2d( const Real& t,
236
237
238
              const Real& /*x*/,
              const Real& /*y*/,
              const Real& /*z*/,
gfourestey's avatar
gfourestey committed
239
              const ID&   id ) const
Tiziano Passerini's avatar
Tiziano Passerini committed
240
    {
gfourestey's avatar
gfourestey committed
241

Tiziano Passerini's avatar
Tiziano Passerini committed
242
243
        switch (id)
        {
uvilla's avatar
uvilla committed
244
        case 0: // x component
Tiziano Passerini's avatar
Tiziano Passerini committed
245
246
            return 0.0;
            break;
uvilla's avatar
uvilla committed
247
        case 2: // z component
Tiziano Passerini's avatar
Tiziano Passerini committed
248
249
250
251
252
            if ( t <= 0.003 )
                return 1.3332e4;
            //      return 0.01;
            return 0.0;
            break;
uvilla's avatar
uvilla committed
253
        case 1: // y component
Tiziano Passerini's avatar
Tiziano Passerini committed
254
255
256
257
258
            return 0.0;
            //      return 1.3332e4;
            //    else
            //      return 0.0;
            break;
gfourestey's avatar
gfourestey committed
259
        }
Tiziano Passerini's avatar
Tiziano Passerini committed
260
261
        return 0;
    }
gfourestey's avatar
gfourestey committed
262

263
    fct_Type getU_2d()
Tiziano Passerini's avatar
Tiziano Passerini committed
264
    {
265
        fct_Type f;
Tiziano Passerini's avatar
Tiziano Passerini committed
266
267
268
        f = boost::bind(&Cylinder::Private::u2d, this, _1, _2, _3, _4, _5);
        return f;
    }
gfourestey's avatar
gfourestey committed
269
270
271
272
273
274

    /**
     * one flat (1,1,1)
     *
     * Define the velocity profile at the inlet for the 2D cylinder
     */
275
    Real poiseuille( const Real& /*t*/,
Tiziano Passerini's avatar
Tiziano Passerini committed
276
277
                     const Real& x,
                     const Real& y,
278
                     const Real& /*z*/,
Tiziano Passerini's avatar
Tiziano Passerini committed
279
                     const ID&   id ) const
gfourestey's avatar
gfourestey committed
280
281
282
    {
        double r = std::sqrt(x*x + y*y);

uvilla's avatar
uvilla committed
283
        if (id == 2)
gfourestey's avatar
gfourestey committed
284
285
286
287
288
            return Um_2d()*2*((D/2.)*(D/2.) - r*r);

        return 0.;
    }

289
    fct_Type getU_pois()
Tiziano Passerini's avatar
Tiziano Passerini committed
290
    {
291
        fct_Type f;
Tiziano Passerini's avatar
Tiziano Passerini committed
292
293
294
        f = boost::bind(&Cylinder::Private::poiseuille, this, _1, _2, _3, _4, _5);
        return f;
    }
gfourestey's avatar
gfourestey committed
295
296
297
298
299
300


    Real oneU( const Real& /*t*/,
               const Real& /*x*/,
               const Real& /*y*/,
               const Real& /*z*/,
301
               const ID&   /*id*/ ) const
Tiziano Passerini's avatar
Tiziano Passerini committed
302
303
304
    {
        //            if (id == 3)
        return 10.;
gfourestey's avatar
gfourestey committed
305

Tiziano Passerini's avatar
Tiziano Passerini committed
306
307
        return -1.;
    }
gfourestey's avatar
gfourestey committed
308

309
    fct_Type getU_one()
Tiziano Passerini's avatar
Tiziano Passerini committed
310
    {
311
        fct_Type f;
Tiziano Passerini's avatar
Tiziano Passerini committed
312
313
314
        f = boost::bind(&Cylinder::Private::oneU, this, _1, _2, _3, _4, _5);
        return f;
    }
gfourestey's avatar
gfourestey committed
315
316
317
318
319


};

Cylinder::Cylinder( int argc,
320
                    char** argv )
Tiziano Passerini's avatar
Tiziano Passerini committed
321
322
        :
        d( new Private )
gfourestey's avatar
gfourestey committed
323
324
325
326
327
328
329
330
{
    GetPot command_line(argc, argv);
    string data_file_name = command_line.follow("data", 2, "-f", "--file");
    GetPot dataFile( data_file_name );
    d->data_file_name = data_file_name;

    d->Re          = dataFile( "fluid/problem/Re", 1. );
    d->nu          = dataFile( "fluid/physics/viscosity", 1. ) /
Tiziano Passerini's avatar
Tiziano Passerini committed
331
                     dataFile( "fluid/physics/density", 1. );
gfourestey's avatar
gfourestey committed
332
333
334
335
336
337
338
339
340
341
342
343
344
    d->H           = 20.;//dataFile( "fluid/problem/H", 20. );
    d->D           =               dataFile( "fluid/problem/D", 1. );
    d->centered    = (bool)        dataFile( "fluid/problem/centered", 0 );
    d->initial_sol = (std::string) dataFile( "fluid/problem/initial_sol", "stokes");
    std::cout << d->initial_sol << std::endl;


#ifdef EPETRA_MPI
    std::cout << "mpi initialization ... " << std::endl;

    //    MPI_Init(&argc,&argv);

    int ntasks = 0;
345
    d->comm.reset( new Epetra_MpiComm( MPI_COMM_WORLD ) );
Tiziano Passerini's avatar
Tiziano Passerini committed
346
347
    if (!d->comm->MyPID())
    {
gfourestey's avatar
gfourestey committed
348
349
350
351
352
353
354
355
        std::cout << "My PID = " << d->comm->MyPID() << " out of " << ntasks << " running." << std::endl;
        std::cout << "Re = " << d->Re << std::endl
                  << "nu = " << d->nu << std::endl
                  << "H  = " << d->H  << std::endl
                  << "D  = " << d->D  << std::endl;
    }
//    int err = MPI_Comm_size(MPI_COMM_WORLD, &ntasks);
#else
356
    d->comm.reset( new Epetra_SerialComm() );
gfourestey's avatar
gfourestey committed
357
358
359
360
361
362
363
364
#endif

}

void
Cylinder::run()

{
Tiziano Passerini's avatar
Tiziano Passerini committed
365
366
367
368
    typedef FESpace< mesh_Type, MapEpetra >       feSpace_Type;
    typedef boost::shared_ptr<feSpace_Type>       feSpacePtr_Type;
    typedef OseenSolver< mesh_Type >::vector_Type vector_Type;
    typedef boost::shared_ptr<vector_Type>        vectorPtr_Type;
gfourestey's avatar
gfourestey committed
369
370
371
372
373
374
375
376
377
    // Reading from data file
    //
    GetPot dataFile( d->data_file_name );

    //    int save = dataFile("fluid/miscellaneous/save", 1);

    bool verbose = (d->comm->MyPID() == 0);

    // Boundary conditions
perego's avatar
perego committed
378
    BCHandler bcH;
gfourestey's avatar
gfourestey committed
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
    BCFunctionBase uZero( zero_scalar );
    std::vector<ID> zComp(1);
    zComp[0] = 3;

    BCFunctionBase uIn  (  d->getU_2d() );
    BCFunctionBase uOne (  d->getU_one() );
    BCFunctionBase uPois(  d->getU_pois() );


    //BCFunctionBase unormal(  d->get_normal() );

    //cylinder

    bcH.addBC( "Inlet",    INLET,    Essential,     Full,     uPois  , 3 );
    bcH.addBC( "Ringin",   RINGIN,   Essential,     Full,     uZero  , 3 );
    bcH.addBC( "Ringout",  RINGOUT,  Essential,     Full,     uZero  , 3 );
    bcH.addBC( "Outlet",   OUTLET,   Natural,     Full,     uZero, 3 );
    bcH.addBC( "Wall",     WALL,     Essential,   Full,     uZero, 3 );

    int numLM = 0;

400
401
    boost::shared_ptr<OseenData> oseenData(new OseenData());
    oseenData->setup( dataFile );
gfourestey's avatar
gfourestey committed
402

403
404
    MeshData meshData;
    meshData.setup(dataFile, "fluid/space_discretization");
gfourestey's avatar
gfourestey committed
405

406
    boost::shared_ptr<mesh_Type> fullMeshPtr ( new mesh_Type( d->comm ) );
407
    readMesh(*fullMeshPtr, meshData);
gfourestey's avatar
gfourestey committed
408

409
    boost::shared_ptr<mesh_Type> meshPtr;
410
411
    {
        MeshPartitioner< mesh_Type >   meshPart(fullMeshPtr, d->comm);
412
        meshPtr = meshPart.meshPartition();
413
    }
gfourestey's avatar
gfourestey committed
414
    if (verbose) std::cout << std::endl;
415
    if (verbose) std::cout << "Time discretization order " << oseenData->dataTimeAdvance()->orderBDF() << std::endl;
gfourestey's avatar
gfourestey committed
416

417
    //oseenData.meshData()->setMesh(meshPtr);
gfourestey's avatar
gfourestey committed
418
419
420
421
422

    std::string uOrder =  dataFile( "fluid/space_discretization/vel_order", "P1");
    if (verbose)
        std::cout << "Building the velocity FE space ... " << std::flush;

423
    feSpacePtr_Type uFESpacePtr( new feSpace_Type(meshPtr,uOrder,3,d->comm) );
gfourestey's avatar
gfourestey committed
424
425
426
427
428
429
430
431
432
433

    if (verbose)
        std::cout << "ok." << std::endl;


    std::string pOrder =  dataFile( "fluid/space_discretization/press_order", "P1");

    if (verbose)
        std::cout << "Building the pressure FE space ... " << std::flush;

434
    feSpacePtr_Type pFESpacePtr( new feSpace_Type(meshPtr,pOrder,1,d->comm) );
gfourestey's avatar
gfourestey committed
435
436
437
438

    if (verbose)
        std::cout << "ok." << std::endl;

Tiziano Passerini's avatar
Tiziano Passerini committed
439
440
    UInt totalVelDof   = uFESpacePtr->map().map(Unique)->NumGlobalElements();
    UInt totalPressDof = pFESpacePtr->map().map(Unique)->NumGlobalElements();
gfourestey's avatar
gfourestey committed
441
442
443



quinodoz's avatar
quinodoz committed
444
445
    if (verbose) std::cout << "Total Velocity DOF = " << totalVelDof << std::endl;
    if (verbose) std::cout << "Total Pressure DOF = " << totalPressDof << std::endl;
gfourestey's avatar
gfourestey committed
446
447
448

    if (verbose) std::cout << "Calling the fluid constructor ... ";

449
    bcH.setOffset( "Inlet", totalVelDof + totalPressDof - 1 );
gfourestey's avatar
gfourestey committed
450

Tiziano Passerini's avatar
Tiziano Passerini committed
451
452
453
454
    OseenSolver< mesh_Type > fluid (oseenData,
                                    *uFESpacePtr,
                                    *pFESpacePtr,
                                    d->comm, numLM);
455
    MapEpetra fullMap(fluid.getMap());
gfourestey's avatar
gfourestey committed
456
457
458
459
460
461
462
463
464
465

    if (verbose) std::cout << "ok." << std::endl;

    fluid.setUp(dataFile);
    fluid.buildSystem();

    MPI_Barrier(MPI_COMM_WORLD);

    // Initialization

466
467
468
    Real dt     = oseenData->dataTime()->timeStep();
    Real t0     = oseenData->dataTime()->initialTime();
    Real tFinal = oseenData->dataTime()->endTime();
gfourestey's avatar
gfourestey committed
469
470
471
472


    // bdf object to store the previous solutions

473
    TimeAdvanceBDFNavierStokes<vector_Type> bdf;
474
    bdf.setup(oseenData->dataTimeAdvance()->orderBDF());
gfourestey's avatar
gfourestey committed
475

476
477
    vector_Type beta( fullMap );
    vector_Type rhs ( fullMap );
gfourestey's avatar
gfourestey committed
478
479

#ifdef HAVE_HDF5
480
    ExporterHDF5<mesh_Type > ensight( dataFile, meshPtr, "cylinder", d->comm->MyPID());
gfourestey's avatar
gfourestey committed
481
#else
482
    ExporterEnsight<mesh_Type > ensight( dataFile, meshPtr, "cylinder", d->comm->MyPID());
gfourestey's avatar
gfourestey committed
483
484
#endif

485
    vectorPtr_Type velAndPressure ( new vector_Type(*fluid.solution(), ensight.mapType() ) );
gfourestey's avatar
gfourestey committed
486

Tiziano Passerini's avatar
Tiziano Passerini committed
487
488
    ensight.addVariable( ExporterData<mesh_Type>::VectorField, "velocity", uFESpacePtr,
                         velAndPressure, UInt(0) );
gfourestey's avatar
gfourestey committed
489

Tiziano Passerini's avatar
Tiziano Passerini committed
490
491
    ensight.addVariable( ExporterData<mesh_Type>::ScalarField, "pressure", pFESpacePtr,
                         velAndPressure, UInt(3*uFESpacePtr->dof().numTotalDof() ) );
gfourestey's avatar
gfourestey committed
492
493
494
495

    // initialization with stokes solution

    if (d->initial_sol == "stokes")
Tiziano Passerini's avatar
Tiziano Passerini committed
496
497
498
    {
        if (verbose) std::cout << std::endl;
        if (verbose) std::cout << "Computing the stokes solution ... " << std::endl << std::endl;
gfourestey's avatar
gfourestey committed
499

500
        oseenData->dataTime()->setTime(t0);
gfourestey's avatar
gfourestey committed
501

Tiziano Passerini's avatar
Tiziano Passerini committed
502
        MPI_Barrier(MPI_COMM_WORLD);
gfourestey's avatar
gfourestey committed
503

Tiziano Passerini's avatar
Tiziano Passerini committed
504
505
        beta *= 0.;
        rhs  *= 0.;
gfourestey's avatar
gfourestey committed
506

Tiziano Passerini's avatar
Tiziano Passerini committed
507
508
        fluid.updateSystem(0, beta, rhs );
        fluid.iterate( bcH );
gfourestey's avatar
gfourestey committed
509
510
511

//    fluid.postProcess();

Tiziano Passerini's avatar
Tiziano Passerini committed
512
513
        *velAndPressure = *fluid.solution();
        ensight.postProcess( 0 );
Tiziano Passerini's avatar
Tiziano Passerini committed
514
        fluid.resetPreconditioner();
Tiziano Passerini's avatar
Tiziano Passerini committed
515
    }
gfourestey's avatar
gfourestey committed
516

517
    bdf.bdfVelocity().setInitialCondition( *fluid.solution() );
gfourestey's avatar
gfourestey committed
518
519
520

    // Temporal loop

521
    LifeChrono chrono;
gfourestey's avatar
gfourestey committed
522
523
524
525
526
    int iter = 1;

    for ( Real time = t0 + dt ; time <= tFinal + dt/2.; time += dt, iter++)
    {

527
        oseenData->dataTime()->setTime(time);
gfourestey's avatar
gfourestey committed
528
529
530
531

        if (verbose)
        {
            std::cout << std::endl;
532
            std::cout << "We are now at time "<< oseenData->dataTime()->time() << " s. " << std::endl;
gfourestey's avatar
gfourestey committed
533
534
535
536
537
            std::cout << std::endl;
        }

        chrono.start();

538
        double alpha = bdf.bdfVelocity().coefficientFirstDerivative( 0 ) / oseenData->dataTime()->timeStep();
539
    //beta = bdf.bdfVelocity().extrapolation(  beta);
540
        bdf.bdfVelocity().extrapolation(beta);
541
        bdf.bdfVelocity().updateRHSContribution( oseenData->dataTime()->timeStep());
542
        rhs  = fluid.matrixMass()*bdf.bdfVelocity().rhsContributionFirstDerivative();
gfourestey's avatar
gfourestey committed
543
544
545
546

        fluid.updateSystem( alpha, beta, rhs );
        fluid.iterate( bcH );

547
        bdf.bdfVelocity().shiftRight( *fluid.solution() );
gfourestey's avatar
gfourestey committed
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568

//         if (((iter % save == 0) || (iter == 1 )))
//         {
        *velAndPressure = *fluid.solution();
        ensight.postProcess( time );
//         }
//         postProcessFluxesPressures(fluid, bcH, time, verbose);


        MPI_Barrier(MPI_COMM_WORLD);

        chrono.stop();
        if (verbose) std::cout << "Total iteration time " << chrono.diff() << " s." << std::endl;
    }

}


//////////////////////