nsip.hpp 14.5 KB
Newer Older
fernandez's avatar
fernandez committed
1
2
3
/*!
  \file NavierStokesSolverIP.h
  \author M.A. Fernandez
prudhomm's avatar
prudhomm committed
4
  \date 6/2003
fernandez's avatar
fernandez committed
5
6
  \version 1.0

prudhomm's avatar
prudhomm committed
7
8
9
  \brief This file contains a Stokes solver class which implements a implicit
         scheme with an exact factorization. Preconditioning of the Schur Complement
         is done by an algebraic Chorin-Temam pressure-corrected preconditioner
fernandez's avatar
fernandez committed
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
*/
#ifndef _NAVIERSTOKESSOLVERIP_H_
#define _NAVIERSTOKESSOLVERIP_H_

#include "NavierStokesHandler.hpp"
#include "elemMat.hpp"
#include "elemVec.hpp"
#include "elemOper.hpp"
#include "values.hpp"
#include "pattern.hpp"
#include "assemb.hpp"
#include "bc_manage.hpp"
#include "az_aztec.h"
#include "algebraic_facto.hpp"
#include "bcCond.hpp"
#include "chrono.hpp"
#include "dataAztec.hpp"
#include "sobolevNorms.hpp"
#include "geoMap.hpp"
prudhomm's avatar
prudhomm committed
29
30
31
32

namespace LifeV
{
/*!
fernandez's avatar
fernandez committed
33
34
  \class NavierStokesSolverIP

prudhomm's avatar
prudhomm committed
35
36
   This class implements an NavierStokes solver via exact factorization. Preconditioning of the
  Schur Complement is done by an algebraic Chorin-Temam pressure-corrected preconditioner
fernandez's avatar
fernandez committed
37
38
39
40
41

*/
template<typename Mesh>
class NavierStokesSolverIP:
public NavierStokesHandler<Mesh> {
prudhomm's avatar
prudhomm committed
42

fernandez's avatar
fernandez committed
43
 public:
prudhomm's avatar
prudhomm committed
44
   typedef  typename  NavierStokesHandler<Mesh>::Function Function;
fernandez's avatar
fernandez committed
45

prudhomm's avatar
prudhomm committed
46
  //! Constructor
fernandez's avatar
fernandez committed
47
48
  /*!
    \param data_file GetPot data file
prudhomm's avatar
prudhomm committed
49
50
51
    \param refFE reference FE
    \param Qr volumic quadrature rule
    \param bdQr surface quadrature rule
fernandez's avatar
fernandez committed
52
53
54
55
56
    \param BCh boundary conditions for the velocity
  */
  NavierStokesSolverIP(const GetPot& data_file, const RefFE& refFE, const QuadRule& Qr,
	    const QuadRule& bdQr, BC_Handler& BCh);

prudhomm's avatar
prudhomm committed
57
58
59
  //! Update the right  hand side  for time advancing
  /*!
    \param source volumic source
fernandez's avatar
fernandez committed
60
61
62
63
64
65
66
67
68
69
70
71
72
    \param time present time
  */
  void timeAdvance(const Function source, const Real& time);

  //! Update convective term, bc treatment and solve the linearized ns system
  void iterate(const Real& time);

  void eval(Vector& fx0,Vector& gx0,Vector x0, int status);

 private:

  //! Block pattern of M_u
  MSRPatt _pattM_u_block;
prudhomm's avatar
prudhomm committed
73

fernandez's avatar
fernandez committed
74
75
76
77
  //! Pattern for M
  MixedPattern<3,3,MSRPatt> _pattM_u;

  //! Block pattern of C
prudhomm's avatar
prudhomm committed
78
  MSRPatt _pattC;
fernandez's avatar
fernandez committed
79
80
81

  //! Matrix M_u: Vmass
  MixedMatr<3,3,MSRPatt,double> _M_u;
prudhomm's avatar
prudhomm committed
82

fernandez's avatar
fernandez committed
83
84
  //! Matrix CStokes
  MSRMatr<double> _CStokes;
prudhomm's avatar
prudhomm committed
85

fernandez's avatar
fernandez committed
86
87
88
89
  //! Matrix C: 1/dt*Vmass + mu*Vstiff operator + Convective_term
  MSRMatr<double> _C;

  //! Elementary matrices and vectors
prudhomm's avatar
prudhomm committed
90
  ElemMat _elmatC; //velocity stiffnes
fernandez's avatar
fernandez committed
91
  ElemMat _elmatM_u; //velocity mass
prudhomm's avatar
prudhomm committed
92
93
  ElemMat _elmatD;
  ElemMat _elmatDtr;
fernandez's avatar
fernandez committed
94
95
96
97
98
  ElemMat _elmatP;
  ElemVec _elvec; // Elementary right hand side

  //! Right  hand  side for the velocity
  PhysVectUnknown<Vector> _f_u;
prudhomm's avatar
prudhomm committed
99

fernandez's avatar
fernandez committed
100
101
102
  //! Right  hand  side global
  Vector _b;

prudhomm's avatar
prudhomm committed
103
104
105
106
  //! Global solution _u and _p
  Vector _x;

  //! Global solution _u and _p
fernandez's avatar
fernandez committed
107
108
109
110
111
112
113
114
  Vector _diff;

  DataAztec _dataAztec;

  Real _time;
};


prudhomm's avatar
prudhomm committed
115
template<typename Mesh> void NavierStokesSolverIP<Mesh>::
fernandez's avatar
fernandez committed
116
117
118
119
eval(Vector& fx0,Vector& gx0,Vector x0, int status) {

  iterate(0.0);
  for (UInt i=0; i < 3*_dim_u ; ++i)
prudhomm's avatar
prudhomm committed
120
121
    fx0[i] = _u[i];

fernandez's avatar
fernandez committed
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
}

//
//                                         IMPLEMENTATION
//
template<typename Mesh> NavierStokesSolverIP<Mesh>::
NavierStokesSolverIP(const GetPot& data_file, const RefFE& refFE, const QuadRule& Qr,
		     const QuadRule& bdQr, BC_Handler& BCh):
  NavierStokesHandler<Mesh>(data_file,refFE,refFE,Qr,bdQr,Qr,bdQr,BCh),
  _pattM_u_block(_dof_u),
  _pattM_u(_pattM_u_block,"diag"),
  _pattC(_dof_u,_mesh,4),
  _M_u(_pattM_u),
  _CStokes(_pattC),
  _C(_pattC),
prudhomm's avatar
prudhomm committed
137
  _elmatC(_fe_u.nbNode,nDimensions,nDimensions),
fernandez's avatar
fernandez committed
138
  _elmatM_u(_fe_u.nbNode,nDimensions,nDimensions),
prudhomm's avatar
prudhomm committed
139
140
141
142
143
  _elmatD(_fe_u.nbNode,4,nDimensions),
  _elmatDtr(_fe_u.nbNode,nDimensions,4),
  _elmatP(_fe_u.nbNode,4,4),
  _elvec(_fe_u.nbNode,nDimensions),
  _f_u(_dim_u),
fernandez's avatar
fernandez committed
144
145
146
147
  _b(4*_dim_u),
  _x(4*_dim_u),
  _diff(4*_dim_u),
  _dataAztec(data_file,"fluid/aztec") {
prudhomm's avatar
prudhomm committed
148
149
150



fernandez's avatar
fernandez committed
151
    cout << endl;
prudhomm's avatar
prudhomm committed
152
    cout << "O-  Pressure unknowns: " << _dim_p     << endl;
fernandez's avatar
fernandez committed
153
    cout << "O-  Velocity unknowns: " << _dim_u     << endl<<endl;
prudhomm's avatar
prudhomm committed
154
155
    cout << "O-  Computing mass and Stokes matrices... ";

fernandez's avatar
fernandez committed
156
157
    Chrono chrono;
    chrono.start();
prudhomm's avatar
prudhomm committed
158
159

    // Matrices initialization
fernandez's avatar
fernandez committed
160
161
    _M_u.zeros();
    _CStokes.zeros();
prudhomm's avatar
prudhomm committed
162
163

    // Number of velocity components
fernandez's avatar
fernandez committed
164
    UInt nc_u=_u.nbcomp();
prudhomm's avatar
prudhomm committed
165

fernandez's avatar
fernandez committed
166
167
168
    //inverse of dt:
    // Real dti=1./_dt;

prudhomm's avatar
prudhomm committed
169
170


fernandez's avatar
fernandez committed
171
172

    Real gamma;
prudhomm's avatar
prudhomm committed
173
174

    // Elementary computation and matrix assembling
fernandez's avatar
fernandez committed
175
176
    // Loop on elements
    for(UInt i = 1; i <= _mesh.numVolumes(); i++){
prudhomm's avatar
prudhomm committed
177

fernandez's avatar
fernandez committed
178
      _fe_u.updateFirstDeriv(_mesh.volumeList(i));
prudhomm's avatar
prudhomm committed
179

fernandez's avatar
fernandez committed
180
181
182
183
      _elmatC.zero();
      _elmatM_u.zero();
      _elmatD.zero();
      _elmatDtr.zero();
prudhomm's avatar
prudhomm committed
184

fernandez's avatar
fernandez committed
185
186
187
      // stiffness strain
      stiff_strain(2.0*_mu,_elmatC,_fe_u);
      //stiff_div(0.5*_fe_u.diameter(),_elmatC,_fe_u);
prudhomm's avatar
prudhomm committed
188

fernandez's avatar
fernandez committed
189
190
191
      // mass
      //mass(_rho*dti,_elmatM_u,_fe_u,0,0,nDimensions);
      // _elmatC.mat() += _elmatM_u.mat();
prudhomm's avatar
prudhomm committed
192
193
194


      for(UInt ic=0;ic<nc_u;ic++){
fernandez's avatar
fernandez committed
195
196
197
198
199
200
	for(UInt jc=0;jc<nc_u;jc++) {
	  // stiffness
	  assemb_mat(_CStokes,_elmatC,_fe_u,_dof_u,ic,jc);
	}
	// mass
	//assemb_mat(_M_u,_elmatM_u,_fe_u,_dof_u,ic,ic);
prudhomm's avatar
prudhomm committed
201

fernandez's avatar
fernandez committed
202
	// div
prudhomm's avatar
prudhomm committed
203
	grad(ic,1.0,_elmatDtr,_fe_u,_fe_u,ic,3);
fernandez's avatar
fernandez committed
204
	div( ic,-1.0,_elmatD  ,_fe_u,_fe_u,nc_u,ic);
prudhomm's avatar
prudhomm committed
205

fernandez's avatar
fernandez committed
206
207
208
	// assembling
	assemb_mat(_CStokes,_elmatDtr,_fe_u,_dof_u,ic,nc_u);
	assemb_mat(_CStokes,_elmatD,_fe_u,_dof_u,nc_u,ic);
prudhomm's avatar
prudhomm committed
209

fernandez's avatar
fernandez committed
210
211
212
      }
    }

prudhomm's avatar
prudhomm committed
213

fernandez's avatar
fernandez committed
214
    cout << endl;
prudhomm's avatar
prudhomm committed
215

fernandez's avatar
fernandez committed
216
217
218
    UInt iElAd1, iElAd2;
    CurrentFE fe1(_refFE_u,getGeoMap(_mesh),_Qr_u);
    CurrentFE fe2(_refFE_u,getGeoMap(_mesh),_Qr_u);
prudhomm's avatar
prudhomm committed
219

fernandez's avatar
fernandez committed
220
    _elmatP.zero();
prudhomm's avatar
prudhomm committed
221
222

    // Elementary computation and matrix assembling
fernandez's avatar
fernandez committed
223
    // Loop on interior faces
prudhomm's avatar
prudhomm committed
224
225
    for (UInt i=_mesh.numBFaces()+1; i<=_mesh.numFaces(); ++i){

fernandez's avatar
fernandez committed
226
227
228
229

    // Updating face staff
    _feBd_u.updateMeas( _mesh.face(i) );

fernandez's avatar
fernandez committed
230
231
232
    gamma = _feBd_u.measure()/32.0;   // P1
    //  gamma = _feBd_u.measure()/128.0; // P2
    // gamma = _feBd_u.measure()*sqrt(_feBd_u.measure())/8.0; // P1 p non smooth
prudhomm's avatar
prudhomm committed
233
234
235
236

    iElAd1 = _mesh.face(i).ad_first();
    iElAd2 = _mesh.face(i).ad_second();

fernandez's avatar
fernandez committed
237
238
239
    fe1.updateFirstDeriv(_mesh.volumeList( iElAd1 ) );
    fe2.updateFirstDeriv(_mesh.volumeList( iElAd2 ) );

prudhomm's avatar
prudhomm committed
240
241
    ipstab_grad(gamma,_elmatP, fe1, fe1, _feBd_u,3,3);
    assemb_mat(_CStokes,_elmatP,fe1,_dof_u,3,3);
fernandez's avatar
fernandez committed
242
243
244

    ipstab_grad(gamma,_elmatP, fe2, fe2, _feBd_u,3,3);
    assemb_mat(_CStokes,_elmatP,fe2,_dof_u,3,3);
prudhomm's avatar
prudhomm committed
245
246

    ipstab_grad(-gamma,_elmatP, fe1, fe2, _feBd_u,3,3);
fernandez's avatar
fernandez committed
247
    assemb_mat(_CStokes,_elmatP,fe1,fe2,_dof_u,3,3);
prudhomm's avatar
prudhomm committed
248

fernandez's avatar
fernandez committed
249
250
    ipstab_grad(-gamma,_elmatP, fe2, fe1, _feBd_u,3,3);
    assemb_mat(_CStokes,_elmatP,fe2,fe1,_dof_u,3,3);
prudhomm's avatar
prudhomm committed
251

fernandez's avatar
fernandez committed
252
253
254
  }

  _x = 0.0;
prudhomm's avatar
prudhomm committed
255

fernandez's avatar
fernandez committed
256
257
258
  chrono.stop();
  cout << "done in " << chrono.diff() << " s." << endl;
  // _CStokes.spy("CS.m");
prudhomm's avatar
prudhomm committed
259
260


fernandez's avatar
fernandez committed
261
262
}

prudhomm's avatar
prudhomm committed
263
template<typename Mesh>
fernandez's avatar
fernandez committed
264
265
266
267
void NavierStokesSolverIP<Mesh>::
timeAdvance(const Function source, const Real& time) {


prudhomm's avatar
prudhomm committed
268
  _time = time;
fernandez's avatar
fernandez committed
269
270
271
272

  cout << endl;
  cout << "O== Now we are at time "<< _time << " s." << endl;

prudhomm's avatar
prudhomm committed
273
  // Number of velocity components
fernandez's avatar
fernandez committed
274
275
276
277
278
279
280
281
  UInt nc_u=_u.nbcomp();

  cout << "  o-  Updating mass term on right hand side... ";

  Chrono chrono;
  chrono.start();

  // Right hand side for the velocity at time
fernandez's avatar
fernandez committed
282
  _f_u=0.0;
fernandez's avatar
fernandez committed
283
284
285
286
287
288

  // loop on volumes: assembling source term
  for(UInt i=1; i<=_mesh.numVolumes(); ++i){
    _elvec.zero();
    _fe_u.updateJacQuadPt(_mesh.volumeList(i));

prudhomm's avatar
prudhomm committed
289
    for (UInt ic=0; ic<nc_u; ++ic){
fernandez's avatar
fernandez committed
290
      compute_vec(source,_elvec,_fe_u,_time,ic); // compute local vector
prudhomm's avatar
prudhomm committed
291
      assemb_vec(_f_u,_elvec,_fe_u,_dof_u,ic); // assemble local vector into global one
fernandez's avatar
fernandez committed
292
293
294
295
    }
  }

  //  For the moment steady: without mass term
prudhomm's avatar
prudhomm committed
296
  //
fernandez's avatar
fernandez committed
297
298
299
300
301
302
303
304
305
  //  _f_u += _M_u * _u;

  chrono.stop();
  cout << "done in " << chrono.diff() << " s." << endl;
}




prudhomm's avatar
prudhomm committed
306
template<typename Mesh>
fernandez's avatar
fernandez committed
307
308
309
310
311
312
313
314
315
316
317
318
void NavierStokesSolverIP<Mesh>::iterate(const Real& time) {

  Chrono  chrono;


  _time = time;
  chrono.start();
  _C=_CStokes;

  chrono.stop();

  cout << "  o-  Stokes matrix was copied in " << chrono.diff() << "s." << endl;
prudhomm's avatar
prudhomm committed
319

fernandez's avatar
fernandez committed
320
  cout << "  o-  Updating convective term... ";
prudhomm's avatar
prudhomm committed
321
322

  // Number of velocity components
fernandez's avatar
fernandez committed
323
324
  UInt nc_u=_u.nbcomp();

prudhomm's avatar
prudhomm committed
325

fernandez's avatar
fernandez committed
326
327
328
329
  chrono.start();

  // loop on volumes
  for(UInt i=1; i<=_mesh.numVolumes(); ++i){
prudhomm's avatar
prudhomm committed
330

fernandez's avatar
fernandez committed
331
    _fe_u.updateFirstDeriv(_mesh.volumeList(i)); // as updateFirstDer
prudhomm's avatar
prudhomm committed
332

fernandez's avatar
fernandez committed
333
334
335
    _elmatC.zero();

    UInt eleID = _fe_u.currentId();
prudhomm's avatar
prudhomm committed
336
    // Non linear term, Semi-implicit approach
fernandez's avatar
fernandez committed
337
338
339
    // ULoc contains the velocity values in the nodes
    for (UInt k=0 ; k<(UInt)_fe_u.nbNode ; k++){
      UInt  iloc = _fe_u.patternFirst(k);
prudhomm's avatar
prudhomm committed
340
341
342
      for (UInt ic=0; ic<nc_u; ++ic){
	UInt ig=_dof_u.localToGlobal(eleID,iloc+1)-1+ic*_dim_u;
	_elvec.vec()[iloc+ic*_fe_u.nbNode] = _rho * _u(ig);
fernandez's avatar
fernandez committed
343
344
345
346
347
348
      }
    }

    // Stabilising term: div u^n u v
    //    mass_divw(0.5*_rho,_elvec,_elmatC,_fe_u,0,0,nc_u);

prudhomm's avatar
prudhomm committed
349
    // loop on components
fernandez's avatar
fernandez committed
350
    for (UInt ic=0; ic<nc_u; ++ic){
prudhomm's avatar
prudhomm committed
351
      // compute local convective term and assembling
fernandez's avatar
fernandez committed
352
353
354
355
356
357
      grad(0,_elvec,_elmatC,_fe_u,_fe_u,ic,ic);
      grad(1,_elvec,_elmatC,_fe_u,_fe_u,ic,ic);
      grad(2,_elvec,_elmatC,_fe_u,_fe_u,ic,ic);
      assemb_mat(_C,_elmatC,_fe_u,_dof_u,ic,ic);
    }
  }
prudhomm's avatar
prudhomm committed
358

fernandez's avatar
fernandez committed
359
360
361
  UInt iElAd1, iElAd2, ig, iFaEl,iloc;
  CurrentFE fe1(_refFE_u,getGeoMap(_mesh),_Qr_u);
  CurrentFE fe2(_refFE_u,getGeoMap(_mesh),_Qr_u);
prudhomm's avatar
prudhomm committed
362
  Real gamma_b,gamma_u ,bmax,bmax_u;
fernandez's avatar
fernandez committed
363
364
365
366
367

  ElemVec beta(_feBd_u.nbNode,nDimensions); // local trace of the velocity

  typedef ID (*FTOP)(ID const _localFace, ID const _point);

fernandez's avatar
fernandez committed
368
  FTOP ftop=0;
fernandez's avatar
fernandez committed
369
370
371
372
373
374
375
376
377
378
379
380
381

  switch( _fe_u.nbNode ) {
  case 4:
    ftop = LinearTetra::fToP;
    break;
  case 10:
    ftop = QuadraticTetra::fToP;
    break;
  case 8:
    ftop = LinearHexa::fToP;
    break;
  case 20:
    ftop = QuadraticHexa::fToP;
prudhomm's avatar
prudhomm committed
382
    break;
fernandez's avatar
fernandez committed
383
384
385
  default:
    ERROR_MSG("This refFE is not allowed with IP stabilisation");
    break;
prudhomm's avatar
prudhomm committed
386
  }
fernandez's avatar
fernandez committed
387
388
389

  _elmatC.zero();

prudhomm's avatar
prudhomm committed
390
  // Elementary computation and matrix assembling
fernandez's avatar
fernandez committed
391
  // Loop on interior faces
prudhomm's avatar
prudhomm committed
392
  for (UInt i=_mesh.numBFaces()+1; i<=_mesh.numFaces(); ++i){
fernandez's avatar
fernandez committed
393
394
395
396
397


    // Updating face staff
    _feBd_u.updateMeas( _mesh.face(i) );

prudhomm's avatar
prudhomm committed
398
399
400
401

    iElAd1 = _mesh.face(i).ad_first();
    iElAd2 = _mesh.face(i).ad_second();

fernandez's avatar
fernandez committed
402
403
404
405
    fe1.updateFirstDeriv(_mesh.volumeList( iElAd1 ));
    fe2.updateFirstDeriv(_mesh.volumeList( iElAd2 ));
    iFaEl = _mesh.face(i).pos_first(); // local id of the face in its adjacent element

prudhomm's avatar
prudhomm committed
406
    for(int in=0; in < _feBd_u.nbNode; ++in) {
fernandez's avatar
fernandez committed
407
      iloc = ftop(iFaEl,in+1);
prudhomm's avatar
prudhomm committed
408
409
410
      for(int ic=0; ic < fe1.nbCoor; ++ic) {
	ig=_dof_u.localToGlobal(iElAd1,iloc+1)-1+ic*_dim_u;
	beta.vec()[ic*_feBd_u.nbNode + in]= _u(ig);
fernandez's avatar
fernandez committed
411
412
      }
    }
prudhomm's avatar
prudhomm committed
413

fernandez's avatar
fernandez committed
414
415
416
417
418
    bmax = abs(beta.vec()[0]);
    for (int l=1; l < int(fe1.nbCoor*_feBd_u.nbNode); ++l) {
      if ( bmax < abs(beta.vec()[l]) )
	bmax = abs(beta.vec()[l]);
    }
prudhomm's avatar
prudhomm committed
419

fernandez's avatar
fernandez committed
420
421
422
    bmax_u = bmax;
    if ( bmax_u < _feBd_u.measure() )
      bmax_u = _feBd_u.measure();
prudhomm's avatar
prudhomm committed
423

fernandez's avatar
fernandez committed
424
425
    gamma_b = 0.125*_feBd_u.measure()/bmax_u;
    gamma_u = 0.125*_feBd_u.measure()*bmax;
prudhomm's avatar
prudhomm committed
426
427
428
    //gamma_u = 0.125*sqrt(_feBd_u.measure())*bmax;


fernandez's avatar
fernandez committed
429
    _elmatC.zero();
prudhomm's avatar
prudhomm committed
430
    //ipstab_grad(gamma_u,_elmatC, fe1, fe1, _feBd_u,0,0,3);
fernandez's avatar
fernandez committed
431
432
433
    ipstab_bgrad(gamma_b,_elmatC, fe1, fe1,beta,_feBd_u,0,0,3);
    ipstab_div(gamma_u,_elmatC, fe1, fe1, _feBd_u);
    for (UInt ic=0; ic<3; ++ic)
prudhomm's avatar
prudhomm committed
434
      for (UInt jc=0; jc<3; ++jc)
fernandez's avatar
fernandez committed
435
	assemb_mat(_C,_elmatC,fe1,_dof_u,ic,jc);
prudhomm's avatar
prudhomm committed
436

fernandez's avatar
fernandez committed
437
    _elmatC.zero();
fernandez's avatar
fernandez committed
438
    //ipstab_grad(gamma_u,_elmatC, fe2, fe2, _feBd_u,0,0,3);
fernandez's avatar
fernandez committed
439
440
441
442
443
444
445
    ipstab_bgrad(gamma_b,_elmatC, fe2, fe2,beta,_feBd_u,0,0,3);
    ipstab_div(gamma_u,_elmatC, fe2, fe2, _feBd_u);
    for (UInt ic=0; ic<3; ++ic)
      for (UInt jc=0; jc<3; ++jc)
	assemb_mat(_C,_elmatC,fe2,_dof_u,ic,jc);

    _elmatC.zero();
prudhomm's avatar
prudhomm committed
446
447
448
449
    //ipstab_grad(-gamma_u,_elmatC, fe1, fe2, _feBd_u,0,0,3);
    ipstab_bgrad(-gamma_b,_elmatC,fe1,fe2,beta,_feBd_u,0,0,3);
    ipstab_div(  -gamma_u,_elmatC,fe1,fe2,_feBd_u);
    for (UInt ic=0; ic<3; ++ic)
fernandez's avatar
fernandez committed
450
451
452
453
      for (UInt jc=0; jc<3; ++jc)
	assemb_mat(_C,_elmatC,fe1,fe2,_dof_u,ic,jc);

    _elmatC.zero();
fernandez's avatar
fernandez committed
454
    //ipstab_grad(-gamma_u,_elmatC, fe2, fe1, _feBd_u,0,0,3);
fernandez's avatar
fernandez committed
455
456
    ipstab_bgrad(-gamma_b,_elmatC, fe2, fe1,beta,_feBd_u,0,0,3);
    ipstab_div(-gamma_u,_elmatC, fe2, fe1, _feBd_u);
prudhomm's avatar
prudhomm committed
457
    for (UInt ic=0; ic<3; ++ic)
fernandez's avatar
fernandez committed
458
      for (UInt jc=0; jc<3; ++jc)
prudhomm's avatar
prudhomm committed
459
460
461
	assemb_mat(_C,_elmatC,fe2,fe1,_dof_u,ic,jc);

  }
fernandez's avatar
fernandez committed
462
463
464
465


  chrono.stop();
  cout << "done in " << chrono.diff() << "s." << endl;
prudhomm's avatar
prudhomm committed
466

fernandez's avatar
fernandez committed
467
468
469
470


  // for BC treatment (done at each time-step)
  cout << "  o-  Applying boundary conditions... ";
prudhomm's avatar
prudhomm committed
471
472
  chrono.start();

fernandez's avatar
fernandez committed
473
474
475
476
477
478
479


  _b = 0.0;
  for (UInt i=0; i<3*_dim_u; ++i) {
    _b[i]=_f_u[i];
  }

prudhomm's avatar
prudhomm committed
480

fernandez's avatar
fernandez committed
481
  // BC manage for the velocity
prudhomm's avatar
prudhomm committed
482
483
  if ( !_BCh_u.bdUpdateDone() )
    _BCh_u.bdUpdate(_mesh, _feBd_u, _dof_u);
fernandez's avatar
fernandez committed
484
485
  bc_manage(_C, _b, _mesh, _dof_u, _BCh_u, _feBd_u, 1.0, _time);

prudhomm's avatar
prudhomm committed
486
487

  //if (_BCh_u.fullEssential())
fernandez's avatar
fernandez committed
488
489
490
491
492
493
494
  _C.diagonalize(3*_dim_u, 1.0, _b, pexact(_mesh.point(1).x(),_mesh.point(1).y() ,_mesh.point(1).z()) );
  //  _C.diagonalize_row(3*_dim_u, 1.0);
  //_b[3*_dim_u]= pexact(_mesh.point(1).x(),_mesh.point(1).y() ,_mesh.point(1).z());

  chrono.stop();
  cout << "done in " << chrono.diff() << "s." << endl;

prudhomm's avatar
prudhomm committed
495
496


fernandez's avatar
fernandez committed
497
498
  // AZTEC specifications for the first system
  int    data_org[AZ_COMM_SIZE];   // data organisation for C
prudhomm's avatar
prudhomm committed
499
  int    proc_config[AZ_PROC_SIZE];// Processor information:
fernandez's avatar
fernandez committed
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
  int    options[AZ_OPTIONS_SIZE]; // Array used to select solver options.
  double params[AZ_PARAMS_SIZE];   // User selected solver paramters.
  double status[AZ_STATUS_SIZE];   // Information returned from AZ_solve()
                                     // indicating success or failure.
  AZ_set_proc_config(proc_config, AZ_NOT_MPI);

  //AZTEC matrix and preconditioner
  AZ_MATRIX *C;
  AZ_PRECOND *prec_C;

  int N_eq= 4*_dim_u; // number of DOF for each component
 // data_org assigned "by hands" while no parallel computation is performed
  data_org[AZ_N_internal]= N_eq;
  data_org[AZ_N_border]= 0;
  data_org[AZ_N_external]= 0;
  data_org[AZ_N_neigh]= 0;
  data_org[AZ_name]= DATA_NAME_AZTEC;

prudhomm's avatar
prudhomm committed
518
  // create matrix and preconditionner
fernandez's avatar
fernandez committed
519
520
521
522
523
  C= AZ_matrix_create(N_eq);
  prec_C= AZ_precond_create(C, AZ_precondition, NULL);

  AZ_set_MSR( C, (int*)_pattC.giveRaw_bindx(), (double*)_C.giveRaw_value(), data_org, 0, NULL, AZ_LOCAL);

prudhomm's avatar
prudhomm committed
524
525
  _dataAztec.aztecOptionsFromDataFile(options,params);

fernandez's avatar
fernandez committed
526
527
528
529
530
531
532
533
534
535

  // ---------------
  // C * V = F
  // ---------------

  for (UInt i=0; i<3*_dim_u; ++i)
    _x[i]=_u[i];

  _diff = _x;

prudhomm's avatar
prudhomm committed
536
  cout << "  o-  Solving system...  ";
fernandez's avatar
fernandez committed
537
538
539
540
  chrono.start();
  AZ_iterate(&_x[0], &_b[0], options, params, status, proc_config, C, prec_C, NULL);
  chrono.stop();
  cout << "done in " << chrono.diff() << " s." << endl;
prudhomm's avatar
prudhomm committed
541

fernandez's avatar
fernandez committed
542
  for (UInt i=0; i<3*_dim_u; ++i) {
prudhomm's avatar
prudhomm committed
543
    _u[i]=_x[i];
fernandez's avatar
fernandez committed
544
545
    //   cout << i+1 << ": " << _x[i] << endl;
  }
prudhomm's avatar
prudhomm committed
546

fernandez's avatar
fernandez committed
547
  for (UInt i=0; i<_dim_u; ++i) {
prudhomm's avatar
prudhomm committed
548
     _p[i]=_x[i+3*_dim_u];
fernandez's avatar
fernandez committed
549
550
551
552
553
554
555
     // cout << i+1 << ": " << _x[i+3*_dim_u] << endl;
  }

  Real norm_p=0.0;
  Real norm_u=0.0;

  for(UInt i = 1; i <= _mesh.numVolumes(); i++){
prudhomm's avatar
prudhomm committed
556

fernandez's avatar
fernandez committed
557
    _fe_u.updateFirstDeriv(_mesh.volumeList(i));
prudhomm's avatar
prudhomm committed
558

fernandez's avatar
fernandez committed
559
560
    norm_p += elem_L2_diff_2(_p,pexact,_fe_u,_dof_u);
    norm_u += elem_L2_diff_2(_u,uexact,_fe_u,_dof_u,0.0,int(nc_u));
prudhomm's avatar
prudhomm committed
561

fernandez's avatar
fernandez committed
562
563
564
565
566
567
568
569
570
  }
  cout << endl;
  cout << " - L2 pressure error = " << sqrt(norm_p) << endl;
  cout << " - L2 velocity error = " << sqrt(norm_u) << endl;
  cout << endl;

  AZ_matrix_destroy(&C);
  AZ_precond_destroy(&prec_C);
}
prudhomm's avatar
prudhomm committed
571
}
fernandez's avatar
fernandez committed
572
#endif