cylinder.cpp 17.2 KB
Newer Older
simone's avatar
simone 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
31
32
33
34
35
36
37
38
39
40
41
/* -*- 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
 */

// #include <lifeconfig.h>

// #include <life/lifecore/life.hpp>
// #include <life/lifecore/GetPot.hpp>
// #include <life/lifecore/debug.hpp>

// #include <life/lifefilters/importer.hpp>

//#include "NavierStokesSolverBlockIP.hpp"

//#include "Epetra_SerialComm.h"
passerini's avatar
passerini committed
42
#include <Epetra_ConfigDefs.h>
simone's avatar
simone committed
43
#ifdef EPETRA_MPI
passerini's avatar
passerini committed
44
45
	#include <Epetra_MpiComm.h>
    #include <mpi.h>
simone's avatar
simone committed
46
#else
passerini's avatar
passerini committed
47
	#include <Epetra_SerialComm.h>
simone's avatar
simone committed
48
49
50
51
52
53
54
55
56
57
58
59
#endif
//#include "life/lifesolver/NavierStokesSolver.hpp"
#include <life/lifearray/EpetraMatrix.hpp>
#include <life/lifealg/EpetraMap.hpp>
#include <life/lifemesh/partitionMesh.hpp>
#include <life/lifesolver/dataNavierStokes.hpp>
#include <life/lifefem/FESpace.hpp>
#include <life/lifefem/bdfNS_template.hpp>
#include <life/lifefilters/ensight.hpp>

#include <life/lifesolver/Oseen.hpp>

passerini's avatar
passerini committed
60
#include "cylinder.hpp"
simone's avatar
simone committed
61
62
#include <iostream>

63
64


simone's avatar
simone committed
65
66
67
68
using namespace LifeV;


//cylinder
simone's avatar
simone committed
69

gfourestey's avatar
update    
gfourestey committed
70
71
//#define CYL2D_MESH_SETTINGS
//#define CYL3D_MESH_SETTINGS
72
#define TUBE20_MESH_SETTINGS
simone's avatar
simone committed
73
74


gfourestey's avatar
update    
gfourestey committed
75
76
77
78
79
#ifdef CYL2D_MESH_SETTINGS // cyl2D.mesh
const int INLET    = 1;
const int SLIPWALL = 20;
const int CYLINDER = 2;
#endif
simone's avatar
simone committed
80
#ifdef CYL3D_2_MESH_SETTINGS // cyl3D-2.mesh
simone's avatar
simone committed
81
82
83
84
85
const int INLET    = 40;
const int WALL     = 60;
const int SLIPWALL = 61;
const int OUTLET   = 50;
const int CYLINDER = 70;
simone's avatar
simone committed
86
87
#endif
#ifdef TUBE20_MESH_SETTINGS
gfourestey's avatar
update    
gfourestey committed
88
89
90
91
92
 const int INLET       = 2;
 const int WALL        = 1;
 const int RINGIN      = 20;
 const int RINGOUT     = 30;
 const int OUTLET      = 3;
simone's avatar
simone committed
93
#endif
simone's avatar
simone committed
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127



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)
{
  switch(i) {
  case 1:
    return 0.0;
    break;
  case 3:
       if ( t <= 0.003 )
          return 1.3332e4;
//      return 0.01;
      return 0.0;
      break;
  case 2:
      return 0.0;
//      return 1.3332e4;
//    else
//      return 0.0;
    break;
  }
  return 0;
}

128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
void
postProcessFluxesPressures( Oseen< RegionMesh3D<LinearTetra> >& nssolver,
    BCHandler& bcHandler,
    const LifeV::Real& t, bool _verbose )
{
  LifeV::Real Q, P;
  UInt flag;

  for( BCHandler::Iterator it = bcHandler.begin();
  it != bcHandler.end(); ++it )
    {
      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("");
      }
    }

}
simone's avatar
simone committed
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184


struct Cylinder::Private
{
    Private() :
        //check(false),
        nu(1),
        //rho(1),
        H(1), D(1)
        //H(20), D(1)
        //H(0.41), D(0.1)
        {}
    typedef boost::function<Real ( Real const&, Real const&, Real const&, Real const&, ID const& )> fct_type;

    double Re;

    std::string data_file_name;

gfourestey's avatar
update    
gfourestey committed
185
    double      nu;  /**< viscosity (in m^2/s) */
simone's avatar
simone committed
186
    //const double rho; /**< density is constant (in kg/m^3) */
gfourestey's avatar
update    
gfourestey committed
187
188
189
190
191
    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;
simone's avatar
simone committed
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245

    Epetra_Comm*   comm;
    /**
     * 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
        {
            if ( id == 1 ) {
                if ( centered ) {
                    return Um_3d() * (H+y)*(H-y) * (H+z)*(H-z) / pow(H,4);
                } else {
                    return 16 * Um_3d() * y * z * (H-y) * (H-z) / pow(H,4);
                }
            } else {
                return 0;
            }
        }

    fct_type getU_3d()
        {
            fct_type f;
            f = boost::bind(&Cylinder::Private::u3d, this, _1, _2, _3, _4, _5);
            return f;
        }

    /**
     * u2d flat 2D velocity profile.
     *
     * Define the velocity profile at the inlet for the 2D cylinder
     */
gfourestey's avatar
update    
gfourestey committed
246
247
    Real u2d( const Real& t,
              const Real& x,
simone's avatar
simone committed
248
              const Real& y,
gfourestey's avatar
update    
gfourestey committed
249
              const Real& z,
simone's avatar
simone committed
250
251
              const ID&   id ) const
        {
gfourestey's avatar
update    
gfourestey committed
252
253
254
255
256

#ifdef CYL2D_MESH_SETTINGS // cyl3D-2.mesh
            if ( id == 1 )
                {
                  return 1.;
simone's avatar
simone committed
257
                  return 1./(20.*20.)*(y + 20.)*(20. - y);
gfourestey's avatar
update    
gfourestey committed
258
                  // 	         if ( centered ) {
simone's avatar
simone committed
259
260
261
262
//                      return Um_2d() * (y+H)*(H-y) / (H*H);
//                  } else {
//                      return 4 * Um_2d() * y * (H-y) / (H*H);
//                  }
gfourestey's avatar
update    
gfourestey committed
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
                }
            else
                {
                    return 0.;
                }
#endif
#ifdef TUBE20_MESH_SETTINGS
            switch(id) {
            case 1: // x component
                return 0.0;
                break;
            case 3: // z component
                if ( t <= 0.003 )
                    return 1.3332e4;
                //      return 0.01;
                return 0.0;
                break;
            case 2: // y component
                return 0.0;
                //      return 1.3332e4;
                //    else
                //      return 0.0;
                break;
            }
            return 0;
#endif


simone's avatar
simone committed
291
292
293
294
295
296
297
298
299
        }

    fct_type getU_2d()
        {
            fct_type f;
            f = boost::bind(&Cylinder::Private::u2d, this, _1, _2, _3, _4, _5);
            return f;
        }

simone's avatar
simone committed
300
301
302
303
304
    /**
     * one flat (1,1,1)
     *
     * Define the velocity profile at the inlet for the 2D cylinder
     */
gfourestey's avatar
update    
gfourestey committed
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
    Real poiseuille( const Real& t,
                      const Real& x,
                      const Real& y,
                      const Real& z,
                      const ID&   id ) const
    {
        double r = std::sqrt(x*x + y*y);

        if (id == 3)
            return Um_2d()*2*((D/2.)*(D/2.) - r*r);

        return 0.;
    }

    fct_type getU_pois()
        {
            fct_type f;
            f = boost::bind(&Cylinder::Private::poiseuille, this, _1, _2, _3, _4, _5);
            return f;
        }


simone's avatar
simone committed
327
328
    Real oneU( const Real& /*t*/,
               const Real& /*x*/,
329
               const Real& /*y*/,
simone's avatar
simone committed
330
               const Real& /*z*/,
gfourestey's avatar
update    
gfourestey committed
331
               const ID&   id ) const
simone's avatar
simone committed
332
        {
gfourestey's avatar
gfourestey committed
333
334
            //            if (id == 3)
                return 10.;
gfourestey's avatar
update    
gfourestey committed
335
336

            return 0.;
simone's avatar
simone committed
337
338
339
340
341
342
343
344
345
        }

    fct_type getU_one()
        {
            fct_type f;
            f = boost::bind(&Cylinder::Private::oneU, this, _1, _2, _3, _4, _5);
            return f;
        }

simone's avatar
simone committed
346
347
348
349
350
351
352
353
354
355
356

};

Cylinder::Cylinder( int argc,
                    char** argv,
                    LifeV::AboutData const& /*ad*/,
                    LifeV::po::options_description const& /*od*/ )
    :
    d( new Private )
{
    GetPot command_line(argc, argv);
357
    string data_file_name = command_line.follow("data", 2, "-f", "--file");
simone's avatar
simone committed
358
359
360
    GetPot dataFile( data_file_name );
    d->data_file_name = data_file_name;

gfourestey's avatar
update    
gfourestey committed
361
362
    d->Re          = dataFile( "fluid/problem/Re", 1. );
    d->nu          = dataFile( "fluid/physics/viscosity", 1. ) /
simone's avatar
simone committed
363
        dataFile( "fluid/physics/density", 1. );
gfourestey's avatar
update    
gfourestey committed
364
365
366
367
368
    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;
simone's avatar
simone committed
369

370

simone's avatar
simone committed
371
372
373
374
375
#ifdef EPETRA_MPI
    std::cout << "mpi initialization ... " << std::endl;

    //    MPI_Init(&argc,&argv);

376
    int ntasks = 0;
simone's avatar
simone committed
377
378
379
380
381
382
383
384
    d->comm = new Epetra_MpiComm( MPI_COMM_WORLD );
    if (!d->comm->MyPID()) {
        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;
    }
385
386
387
388
389
//    int err = MPI_Comm_size(MPI_COMM_WORLD, &ntasks);
#else
    d->comm = new Epetra_SerialComm();
#endif

simone's avatar
simone committed
390
391
392
393
394
395
396
397
398
399
400
}

void
Cylinder::run()

{

    typedef Oseen< RegionMesh3D<LinearTetra> >::vector_type  vector_type;
    typedef boost::shared_ptr<vector_type> vector_ptrtype;
    // Reading from data file
    //
401
    GetPot dataFile( d->data_file_name );
simone's avatar
simone committed
402

403
//    int save = dataFile("fluid/miscellaneous/save", 1);
simone's avatar
simone committed
404
405
406
407
408
409
410
411
412

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

    // Boundary conditions
    BCHandler bcH( 5, BCHandler::HINT_BC_NONE );
    BCFunctionBase uZero( zero_scalar );
    std::vector<ID> zComp(1);
    zComp[0] = 3;

simone's avatar
simone committed
413
    BCFunctionBase uIn  (  d->getU_2d() );
gfourestey's avatar
update    
gfourestey committed
414
415
416
417
418
    BCFunctionBase uOne (  d->getU_one() );
    BCFunctionBase uPois(  d->getU_pois() );


#ifdef CYL2D_MESH_SETTINGS // cyl3D-2.mesh
simone's avatar
simone committed
419
420

    //cylinder
gfourestey's avatar
update    
gfourestey committed
421
    bcH.addBC( "Inlet",    INLET,    Essential, Full,      uIn,   3 );
simone's avatar
simone committed
422
423
424
//     if ( WALL != INLET )
    bcH.addBC( "Slipwall", SLIPWALL, Essential, Component, uZero, zComp );
    bcH.addBC( "Cylinder", CYLINDER, Essential, Full,      uZero, 3 );
gfourestey's avatar
update    
gfourestey committed
425
    //    bcH.addBC( "Slipwall", SLIPWALL, Essential, Full, uZero , 3 );
simone's avatar
simone committed
426
#endif
gfourestey's avatar
gfourestey committed
427

simone's avatar
simone committed
428
429
430
431
432
#ifdef TUBE20_MESH_SETTINGS
    //BCFunctionBase unormal(  d->get_normal() );

    //cylinder

gfourestey's avatar
gfourestey committed
433
434
    //bcH.addBC( "Inlet",    INLET,    Natural,   Full,      uIn, 3 );
    bcH.addBC( "Inlet",    INLET,    Flux,        Full,     uOne, 3);
gfourestey's avatar
update    
gfourestey committed
435
    bcH.addBC( "Outlet",   OUTLET,   Natural,     Full,     uZero, 3 );
simone's avatar
simone committed
436
    //bcH.addBC( "Wall",     WALL,     Natural,     Full,      uNormal, 3 );
gfourestey's avatar
update    
gfourestey committed
437
438
439
440
    //bcH.addBC( "Wall",     WALL,     Natural,     Full,      uNormal, 3 );
    bcH.addBC( "Wall",     WALL,     Essential,   Full,     uZero, 3 );
    bcH.addBC( "RingIn",   RINGIN,  Essential,   Full,      uZero, 3 );
    bcH.addBC( "RingOut",  RINGOUT, Essential,   Full,      uZero, 3 );
simone's avatar
simone committed
441
442
#endif

gfourestey's avatar
gfourestey committed
443
444
    int numLM = 1;

simone's avatar
simone committed
445
446
447
448
449
    DataNavierStokes<RegionMesh3D<LinearTetra> > dataNavierStokes( dataFile );

    partitionMesh< RegionMesh3D<LinearTetra> >   meshPart(*dataNavierStokes.mesh(), *d->comm);

    if (verbose) std::cout << std::endl;
malossi's avatar
malossi committed
450
    if (verbose) std::cout << "Time discretization order " << dataNavierStokes.getBDF_order() << std::endl;
simone's avatar
simone committed
451
452
453

    dataNavierStokes.setMesh(meshPart.mesh());

454
    std::string uOrder =  dataFile( "fluid/space_discretization/vel_order", "P1");
simone's avatar
simone committed
455
456
    if (verbose)
        std::cout << "Building the velocity FE space ... " << std::flush;
457
458

    FESpace< RegionMesh3D<LinearTetra>, EpetraMap > uFESpace(meshPart,uOrder,3,*d->comm);
simone's avatar
simone committed
459
460
461
462

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

463
464
465

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

simone's avatar
simone committed
466
467
468
    if (verbose)
        std::cout << "Building the pressure FE space ... " << std::flush;

469
    FESpace< RegionMesh3D<LinearTetra>, EpetraMap > pFESpace(meshPart,pOrder,1,*d->comm);
simone's avatar
simone committed
470
471
472
473

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

simone's avatar
simone committed
474
475
    UInt totalVelDof   = uFESpace.map().getMap(Unique)->NumGlobalElements();
    UInt totalPressDof = pFESpace.map().getMap(Unique)->NumGlobalElements();
simone's avatar
simone committed
476
477


gfourestey's avatar
gfourestey committed
478

simone's avatar
simone committed
479
480
481
482
483
    if (verbose) std::cout << "Total Velocity Dof = " << totalVelDof << std::endl;
    if (verbose) std::cout << "Total Pressure Dof = " << totalPressDof << std::endl;

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

gfourestey's avatar
gfourestey committed
484
485
    bcH.setOffset("Inlet", totalVelDof + totalPressDof);

simone's avatar
simone committed
486
487
488
    Oseen< RegionMesh3D<LinearTetra> > fluid (dataNavierStokes,
                                              uFESpace,
                                              pFESpace,
gfourestey's avatar
gfourestey committed
489
490
                                              *d->comm,
                                              numLM);
simone's avatar
simone committed
491
492
    EpetraMap fullMap(fluid.getMap());

simone's avatar
simone committed
493

gfourestey's avatar
update    
gfourestey committed
494
495
496
497
498
499
// #ifdef TUBE20_MESH_SETTINGS
//     vector_type vec_lambda_aux(fullMap),	mixtevec_aux(fullMap);
//     vec_lambda_aux.getEpetraVector().PutScalar(1);
//     vec_lambda_aux.GlobalAssemble();
//     mixtevec_aux.getEpetraVector().PutScalar(1);
//     mixtevec_aux.GlobalAssemble();
simone's avatar
simone committed
500

gfourestey's avatar
update    
gfourestey committed
501
//     vector_type vec_lambda(fluid.getMap(), Repeated), mixtevec(fluid.getMap(), Repeated);
simone's avatar
simone committed
502

gfourestey's avatar
update    
gfourestey committed
503
504
//     vec_lambda = vec_lambda_aux;
//     mixtevec = mixtevec_aux;
simone's avatar
simone committed
505

gfourestey's avatar
update    
gfourestey committed
506
507
//     // Robin BC
//     BCVector robin_wall(vec_lambda, uFESpace.dof().numTotalDof());
simone's avatar
simone committed
508

gfourestey's avatar
update    
gfourestey committed
509
510
//     // Neumann BC, normal component
//     //BCVector robin_wall(vec_lambda, uFESpace.dof().numTotalDof(),1);
simone's avatar
simone committed
511

gfourestey's avatar
update    
gfourestey committed
512
513
//     robin_wall.setMixteCoef(0);
//     robin_wall.setMixteVec(mixtevec);
simone's avatar
simone committed
514

gfourestey's avatar
update    
gfourestey committed
515
516
//     vec_lambda_aux.spy("lambda");
//     vec_lambda.spy("lambdaRep");
simone's avatar
simone committed
517
518
519



gfourestey's avatar
update    
gfourestey committed
520
521
522
//     bcH.addBC( "Wall",      1, Mixte, Full,  robin_wall, 3 );
//     //bcH.addBC( "Wall",      1, Natural, Full,  robin_wall, 3 );
// #endif
simone's avatar
simone committed
523
524
525
526




simone's avatar
simone committed
527
528
529
530
531
532
533
534
535
    if (verbose) std::cout << "ok." << std::endl;

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

    MPI_Barrier(MPI_COMM_WORLD);

    // Initialization

malossi's avatar
malossi committed
536
537
538
    Real dt     = dataNavierStokes.getTimeStep();
    Real t0     = dataNavierStokes.getInitialTime();
    Real tFinal = dataNavierStokes.getEndTime();
simone's avatar
simone committed
539
540
541
542


    // bdf object to store the previous solutions

malossi's avatar
malossi committed
543
    BdfTNS<vector_type> bdf(dataNavierStokes.getBDF_order());
simone's avatar
simone committed
544
545
546
547

    vector_type beta( fullMap );
    vector_type rhs ( fullMap );

gfourestey's avatar
update    
gfourestey committed
548
    vector_ptrtype velAndPressure ( new vector_type(fluid.solution(), Repeated ) );
simone's avatar
simone committed
549
550
551
552
553
554
555
556
557
558

    Ensight<RegionMesh3D<LinearTetra> > ensight( dataFile, meshPart.mesh(), "cylinder", d->comm->MyPID());


    ensight.addVariable( ExporterData::Vector, "velocity", velAndPressure,
                         UInt(0), uFESpace.dof().numTotalDof() );

    ensight.addVariable( ExporterData::Scalar, "pressure", velAndPressure,
                         UInt(3*uFESpace.dof().numTotalDof()),
                         UInt(3*uFESpace.dof().numTotalDof()+pFESpace.dof().numTotalDof()) );
gfourestey's avatar
update    
gfourestey committed
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584

    // initialization with stokes solution

    if (d->initial_sol == "stokes")
        {
            if (verbose) std::cout << std::endl;
            if (verbose) std::cout << "Computing the stokes solution ... " << std::endl << std::endl;

            dataNavierStokes.setTime(t0);

            MPI_Barrier(MPI_COMM_WORLD);

            beta *= 0.;
            rhs  *= 0.;

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

//    fluid.postProcess();

            *velAndPressure = fluid.solution();
            ensight.postProcess( 0 );
            fluid.resetPrec();
        }

    bdf.bdf_u().initialize_unk( fluid.solution() );
simone's avatar
simone committed
585
586
587
588
589
590
591
592
593
594
595
596
597
598

    // Temporal loop

    Chrono chrono;
    int iter = 1;

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

        dataNavierStokes.setTime(time);

        if (verbose)
        {
            std::cout << std::endl;
malossi's avatar
malossi committed
599
            std::cout << "We are now at time "<< dataNavierStokes.getTime() << " s. " << std::endl;
simone's avatar
simone committed
600
601
602
603
604
            std::cout << std::endl;
        }

        chrono.start();

malossi's avatar
malossi committed
605
        double alpha = bdf.bdf_u().coeff_der( 0 ) / dataNavierStokes.getTimeStep();
simone's avatar
simone committed
606
607
608

        beta = bdf.bdf_u().extrap();

malossi's avatar
malossi committed
609
        rhs  = fluid.matrMass()*bdf.bdf_u().time_der( dataNavierStokes.getTimeStep() );
simone's avatar
simone committed
610
611
612
613
614
615
616
617
618
619
620

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

        bdf.bdf_u().shift_right( fluid.solution() );

//         if (((iter % save == 0) || (iter == 1 )))
//         {
        *velAndPressure = fluid.solution();
        ensight.postProcess( time );
//         }
gfourestey's avatar
update    
gfourestey committed
621
//         postProcessFluxesPressures(fluid, bcH, time, verbose);
simone's avatar
simone committed
622
623
624
625
626
627
628
629
630
631
632
633
634
635


        MPI_Barrier(MPI_COMM_WORLD);

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

}


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