DarcySolver.hpp 60.3 KB
Newer Older
1
2
3
//@HEADER
/*
*******************************************************************************
4

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

8
   This file is part of the LifeV library
9

10
11
12
13
   LifeV is free software; you can redistribute it and/or modify
   it under the terms of the GNU Lesser General Public License as published by
   the Free Software Foundation; either version 3 of the License, or
   (at your option) any later version.
14

15
16
17
18
   LifeV is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Lesser General Public License for more details.
19

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


24
*******************************************************************************
25
*/
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
//@HEADER

/*!
 *   @file
     @brief This file contains a Darcy solver with mixed-hybrid finite elements

     @date 05/2010
     @author A. Fumagalli <alessio.fumagalli@mail.polimi.it>

     @contributor M. Kern <michel.kern@inria.fr>
     @maintainer M. Kern <michel.kern@inria.fr>
 */

#ifndef _DARCYSOLVER_H_
#define _DARCYSOLVER_H_ 1
41

42
43
#include <Epetra_LAPACK.h>
#include <Epetra_BLAS.h>
fumagalli's avatar
fumagalli committed
44

45
46
#include <lifev/core/fem/AssemblyElemental.hpp>
#include <lifev/core/fem/BCManage.hpp>
47

48
#include <lifev/core/algorithm/SolverAztecOO.hpp>
fumagalli's avatar
fumagalli committed
49

50
//
51
#include <lifev/darcy/solver/DarcyData.hpp>
52
53
//

54

55
56
57
58
59
60
61
62
63
// Local namespace to store the default source term and the default permeability tensor.
namespace
{

/*! @class Default Darcy source term. Needed in DarcySolver class.
  This class implement the default source term for the Darcy problem, it is the zero function.
*/
struct DarcyDefaultSource
{
64
    LifeV::Real operator()( const LifeV::Real&, const LifeV::Real&,
fumagalli's avatar
fumagalli committed
65
66
                const LifeV::Real&, const LifeV::Real&,
                const LifeV::UInt&) const
67
    {
Tiziano Passerini's avatar
Tiziano Passerini committed
68
        return static_cast<LifeV::Real>( 0 );
69
70
71
72
73
74
75
76
    }
};

}

// LifeV namespace.
namespace LifeV
{
77
78
//! class DarcySolver   This class implements a mixed-hybrid FE Darcy solver.

79
/*!
80
  @author A. Fumagalli <alessio.fumagalli@mail.polimi.it>
81
  @see For applications related to two-phase flow see \cite Fumagalli2011a
82
83

  This class implements a Darcy solver.
fumagalli's avatar
fumagalli committed
84
85
  <br>
  <br>
86
87
88
89
90
91
92
  The classical formulation of this problem is a couple of differential equations of first order with two
  unknowns: \f$ p \in C^1 (\Omega ) \f$, being the pressure or the primal unknown,
  and \f$ \sigma \in (C^1( \Omega ) )^n \f$, being the Darcy velocity or the total flux or the dual unknown,
  such that
  \f[
  \left\{
  \begin{array}{l l l}
93
94
95
96
  \Lambda^{-1} \sigma + \nabla p = f_v & \mathrm{in} & \Omega\,,  \vspace{0.2cm} \\
  \nabla \cdot \sigma - f = 0          & \mathrm{in} & \Omega\,,  \vspace{0.2cm} \\
  p = g_D                              & \mathrm{on} & \Gamma_D\,,\vspace{0.2cm} \\
  \sigma \cdot n + h p = g_R           & \mathrm{on} & \Gamma_R\,.
97
98
99
  \end{array}
  \right.
  \f]
100
101
102
103
104
105
  The first equation is the Darcy equation and the second equation is the conservation law. 
  <br>
  The data in the system are: \f$ \Lambda \f$ the permeability tensor, \f$ f \f$ the source term, \f$ f_v \f$ the
  vector source term, \f$ \Gamma_D \f$ the subset of the boundary of \f$ \Omega \f$ with Dirichlet boundary conditions,
  with datum \f$ g_D \f$, and \f$ \Gamma_R \f$ that is the part of the boundary of \f$ \Omega \f$ with Robin, or Neumann, 
  boundary conditions with data \f$ h \f$ and \f$ g_R \f$. We suppose that \f$ \partial \Omega = \Gamma_D \cup \Gamma_R \f$ and
106
  \f$ \Gamma_D \cap \Gamma_R = \emptyset \f$.
fumagalli's avatar
fumagalli committed
107
  <br>
108
109
110
111
112
113
114
115
116
  Using the hybridization procedure, and introducing a new variable, we may split the problem from
  the entire domain to problems in each elements of the triangulation \f$ \mathcal{T}_h \f$, then we may write
  the weak formulation for the dual problem.
  The new variable is the hybrid unknown \f$ \lambda \f$, e.g. "the trace" of pressure of the primal unknown,
  it is the Lagrange multipliers forcing the continuity of the flux across each faces or edges in the triangulation.
  We introduce the following functional spaces, the first is the space of the primal variable, the third for
  the dual variable, the fourth for the hybrid variable.
  \f[
  \begin{array}{l}
fumagalli's avatar
fumagalli committed
117
  V = L^2 (\Omega ) \,,\vspace{0.2cm} \                 \
118
119
120
  H(div, K) = \displaystyle \left\{ \tau \in ( L^2(K))^n : \nabla \cdot \tau \in L^2(K)\right\}\,, \vspace{0.2cm}\\
  Z = \displaystyle \left\{ \tau \in L^2(\Omega) : \tau\vert_K \in H(div, K) \, \forall K \in  \mathcal{T}_h \right\}\,, \vspace{0.2cm}\\
  \Lambda = \displaystyle \left\{ \lambda \in \prod_{K \in \mathcal{T}_h} H^{1/2} (\partial K): \lambda_K = \lambda_{K'} \
121
  ,\, \mathrm{on} \,\, e_{K\cap K'} \, \forall K \in \mathcal{T}_h,\, \lambda = g_D \,\, \mathrm{on} \,\, \Gamma_D \right\}\,.
122
123
124
125
126
127
128
129
130
131
  \end{array}
  \f]
  Introducing the following bilinear forms and functionals
  \f[
  \begin{array}{l l}
  a(\sigma, \tau) = \displaystyle \sum_{K \in \mathcal{T}_h}\int_K \Lambda^{-1} \sigma \cdot \tau \,,  &
  b(p, \tau) = \displaystyle -\sum_{K \in \mathcal{T}_h} \int_K p \nabla \cdot \tau\,, \vspace{0.2cm}\\
  c(\lambda, \tau) = \displaystyle \sum_{K \in \mathcal{T}_h} \int_{\partial K} \lambda \tau \cdot n\,,&
  h(\lambda, \mu) = \displaystyle \sum_{K \in \mathcal{T}_h} \int_{\partial K \cap \Gamma_R} h \mu \lambda \,,\vspace{0.2cm}\\
  F(v) = \displaystyle \sum_{K \in \mathcal{T}_h} \int_K f v\,,&
132
  F_v(\tau) = \displaystyle \sum_{K \in \mathcal{T}_h} \int_K f_v \tau \,, \vspace{0.2cm}\\
133
134
135
136
137
138
139
  G(\mu) = \displaystyle \sum_{K \in \mathcal{T}_h} \int_{\partial K \cap \Gamma_R} g \mu\,,
  \end{array}
  \f]
  we obtain the Darcy problem in the weak form: find \f$ (\sigma, \, p, \, \lambda) \in Z \times V \times \Lambda \f$ such that
  \f[
  \left\{
  \begin{array}{l l}
140
141
142
  a(\sigma, \tau) + b(p, \tau) + c(\lambda, \tau) = F_v(\tau)\,,  & \forall \tau \in Z \,,\vspace{0.2cm}\\
  b(v, \sigma) = - F(v)\,,                                        & \forall v \in V \,,\vspace{0.2cm}\\
  c(\mu, \sigma) + h(\lambda, \mu) = G(\mu) \,,                   & \forall \mu \in \Lambda\,.
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
  \end{array}
  \right.
  \f]
  At discrete level we introduce the polynomial space, of degree \f$ r \f$, that approximate the finite dimensional
  spaces introduced above \f$ V_h \subset V \f$, \f$ Z_h \subset Z \f$ and \f$\Lambda_h \subset \Lambda \f$
  \f[
  \begin{array}{l}
  V_h = \displaystyle \left\{ v_h \in V: v_h|_K \in P_r (K)\, \forall K \in \mathcal{T}_h \right\}\,, \vspace{0.2cm}\\
  Z_h = \displaystyle \left\{ \tau_h \in Z: \tau_h|_K \in RT_r(K) \, \forall K \in \mathcal{T}_h \right\} \,, \vspace{0.2cm} \\
  \Lambda_h = \displaystyle \left\{ \lambda_h \in \Lambda: \lambda_h|_{\partial K} \in R_r( \partial K ) \, \forall K \in \mathcal{T}_h\right\}\,,
  \end{array}
  \f]
  where \f$ P_r(K) \f$ is the space of polynomial of degree \f$ r \f$ in the element \f$ K \f$, \f$ RT_r(K) \f$ is the space of
  polynomial of Raviart-Thomas of degrees \f$ r \f$ in the element \f$ K \f$ and \f$ R_r(\partial K) \f$ is the space of polynomial of
  degree \f$ r \f$ definite on each of the boundary of the element \f$ K \f$ and discontinuous from one edge to the other.
  <BR>
  The finite dimensional problem is: find \f$ (\sigma_h,\, p_h, \, \lambda_h) \in Z_h \times V_h \times \Lambda_h \f$ such that
  \f[
  \left\{
  \begin{array}{l l}
163
164
165
  a(\sigma_h, \tau_h) + b(p_h, \tau_h) + c(\lambda_h, \tau_h) = F_v(\tau_h) \,,  & \forall \tau_h \in Z_h \,,\vspace{0.2cm}\\
  b(v_h, \sigma_h) = -F(v_h)\,,                                                  & \forall v_h \in V_h \,,\vspace{0.2cm}\\
  c(\mu_h, \sigma_h) + h(\lambda_h, \mu_h) = G(\mu_h) \,,                        & \forall \mu_h \in \Lambda_h\,.
166
167
168
169
170
171
172
173
  \end{array}
  \right.
  \f]
  To solve the problem we use the static condensation procedure, i.e. the unknowns in the discrete
  weak system are not independent and \f$ p_K \f$, \f$\sigma_K \f$ may be written in function of
  \f$ \lambda_K \f$ alone. We introduce the following local matrices
  \f[
  \begin{array}{l l l}
174
175
176
177
178
179
180
  \left[ A \right]_{ij}  = \displaystyle   \int_K \Lambda^{-1} \psi_j \cdot \psi_i \,, &
  \left[ B \right]_{ij}  = \displaystyle - \int_K \phi_j \nabla \cdot \psi_i \,, &
  \left[ C \right]_{ij}  = \displaystyle   \int_{\partial K} \xi_i \psi_j \cdot n \,, \vspace{0.2cm} \\
  \left[ H \right]_{ij}  = \displaystyle   \int_{\partial K \cap \Gamma_R} h \xi_i \xi_j\,, &
  \left[ F \right]_{j}   = \displaystyle   \int_K f \phi_j\,, &
  \left[ F_v \right]_{j} = \displaystyle   \int_K f_v \psi_j\,, \vspace{0.2cm} \\
  \left[ G \right]_{j}   = \displaystyle   \int_{\partial K \cap \Gamma_R } g \xi_j\,,
181
182
183
184
185
186
187
  \end{array}
  \f]
  where we avoid to write the dependence on the triangle \f$ K \f$ in all the matrices and vectors. <BR>
  The local matrix formulation of the finite dimensional problem is
  \f[
  \left\{
  \begin{array}{l}
188
189
  A \sigma_K + B p_K + C \lambda_K = F_v\,, \vspace{0.2cm} \\
  B^T \sigma_K = -F \,,                     \vspace{0.2cm}\\
fumagalli's avatar
fumagalli committed
190
  C^T \sigma_K + H \lambda_K = G\,.
191
192
193
194
195
196
197
198
  \end{array}
  \right.
  \f]
  Or alternatively
  \f[
  \begin{array}{l l l}
  \left[
  \begin{array}{c c c}
fumagalli's avatar
fumagalli committed
199
200
201
  A   & B &  C \vspace{0.2cm} \\
  B^T & 0 &  0 \vspace{0.2cm} \\
  C^T & 0 & H
202
203
204
205
206
207
208
209
210
211
212
213
  \end{array}
  \right] \, \cdot &
  \left[
  \begin{array}{c}
  \sigma_K  \vspace{0.2cm}\\
  p_K       \vspace{0.2cm}\\
  \lambda_K
  \end{array}
  \right] \, &
  =
  \left[
  \begin{array}{c}
214
  F_v \vspace{0.2cm}\\
fumagalli's avatar
fumagalli committed
215
216
  -F \vspace{0.2cm}\\
  G
217
218
219
220
221
222
223
  \end{array}
  \right]\,.
  \end{array}
  \f]
  Introducing the local hybrid matrix and local hybrid right hand side
  \f[
  \begin{array}{l}
fumagalli's avatar
fumagalli committed
224
  L_K = - C^T A^{-1} C + C^T A^{-1} B \left( B^T A^{-1} B \right)^{-1} B^T A^{-1} C + H \,, \vspace{0.2cm} \\
225
  r_K = G + C^T A^{-1} B \left( B^T A^{-1} B \right)^{-1} F + C^T A^{-1} B \left( B^T A^{-1} B \right)^{-1} B^T A^{-1} F_v - C^T A^{-1} F_v \,,
226
227
228
229
230
231
232
233
234
  \end{array}
  \f]
  Imposing that at each edge or face the hybrid unknown is single value we obtain a linear system for the hybrid unknown
  \f[
  L \lambda = r \,.
  \f]
  We recover the primal and dual variable as a post-process from the hybrid variable at element level, so we have
  \f[
  \begin{array}{l}
235
236
  p_K = \left( B^T A^{-1} B \right)^{-1} \left( F + B^T A^{-1} F_v - B^T A^{-1} C \lambda_K \right)\,, \vspace{0.2cm} \\
  \sigma_K = -A^{-1} \left( B p_K - F_v + C \lambda_K \right) \,.
237
238
239
240
241
  \end{array}
  \f]
  @note In the code we do not use the matrix \f$ H \f$ and the vector \f$ G \f$, because all the boundary
  conditions are imposed via BCHandler class.
  @todo Insert any scientific publications that use this solver.
fumagalli's avatar
fumagalli committed
242
  @bug If the save flag for the exporter is setted to 0 the program fails.
243

244
*/
245

grandper's avatar
grandper committed
246
template< typename Mesh, typename SolverType = LifeV::SolverAztecOO >
247
248
249
250
251
class DarcySolver
{

public:

252
    //! @name Public Types
253
254
    //@{

255
256
    typedef SolverType                             solver_Type;

257
    typedef boost::function<Real ( const Real&, const Real&, const Real&,
258
259
260
261
                                   const Real&, const UInt& )> function_Type;

    typedef boost::function<Vector ( const Real&, const Real&, const Real&,
                                     const Real&, const UInt& )> vectorFunction_Type;
262

263
264
    typedef inversePermeability<Mesh, SolverType>  permeability_Type;
    typedef boost::shared_ptr<permeability_Type>   permeabilityPtr_Type;
265

Michel Kern's avatar
Michel Kern committed
266
    typedef DarcyData<Mesh>                        data_Type;
267

268
    typedef Mesh                                   mesh_Type;
269

270
271
    typedef BCHandler                              bchandler_raw_Type;
    typedef boost::shared_ptr<bchandler_raw_Type>  bchandler_Type;
272

273
274
    typedef typename solver_Type::matrix_type      matrix_Type;
    typedef boost::shared_ptr<matrix_Type>         matrixPtr_Type;
275

276
277
    typedef typename solver_Type::vector_type      vector_Type;
    typedef boost::shared_ptr<vector_Type>         vectorPtr_Type;
278

279
280
    typedef typename solver_Type::prec_raw_type    prec_raw_Type;
    typedef typename solver_Type::prec_type        prec_Type;
281

282
283
    typedef Epetra_Comm                            comm_Type;
    typedef boost::shared_ptr< comm_Type >         commPtr_Type;
fumagalli's avatar
fumagalli committed
284

285
286
287
288
289
    //@}

    //! @name Constructors and destructor
    //@{

290
291
    //! Full constructor for the class.

292
293
294
295
296
297
298
    /*!
      @param dataFile Data for the problem.
      @param primal_FESpace Primal finite element space.
      @param dual_FESpace Dual element space.
      @param hybrid_FESpace Hybrid finite element space.
      @param VdotN_FESpace Dual basis function dot outward unit normal at each face (3D) or edge (2D) finite element space.
      @param bcHandler Boundary conditions for the problem.
299
      @param comm Shared pointer of the Epetra communicator.
300
    */
Michel Kern's avatar
Michel Kern committed
301
    DarcySolver ( const data_Type&                dataFile,
302
303
304
305
                        FESpace<Mesh, MapEpetra>& primal_FESpace,
                        FESpace<Mesh, MapEpetra>& dual_FESpace,
                        FESpace<Mesh, MapEpetra>& hybrid_FESpace,
                        FESpace<Mesh, MapEpetra>& VdotN_FESpace,
Michel Kern's avatar
Michel Kern committed
306
307
                  const bchandler_raw_Type&       bcHandler,
                  const commPtr_Type&             comm );
308
309

    //! Constructor for the class without the definition of the boundary handler.
310
311
312
313
314
315
316

    /*!
      @param dataFile Data for the problem.
      @param primal_FESpace Primal finite element space.
      @param dual_FESpace Dual finite element space.
      @param hybrid_FESpace Hybrid finite element space.
      @param VdotN_FESpace Dual basis function dot outward unit normal at each face (3D) or edge (2D) finite element space.
317
      @param comm Shared pointer of the Epetra communicator.
318
    */
Michel Kern's avatar
Michel Kern committed
319
    DarcySolver ( const data_Type&                dataFile,
320
321
322
323
                        FESpace<Mesh, MapEpetra>& primal_FESpace,
                        FESpace<Mesh, MapEpetra>& dual_FESpace,
                        FESpace<Mesh, MapEpetra>& hybrid_FESpace,
                        FESpace<Mesh, MapEpetra>& VdotN_FESpace,
Michel Kern's avatar
Michel Kern committed
324
                  const commPtr_Type&             comm );
325
326

    //! Virtual destructor.
fumagalli's avatar
fumagalli committed
327
    virtual ~DarcySolver ();
328
329
330

    //@}

331
    //! @name Methods
332
333
    //@{

334
335
336
337
    //! Build the global hybrid system, the right hand and apply the boundary conditions.
    void buildSystem ();

    //! Set up the linear solver and the preconditioner for the linear system.
fumagalli's avatar
fumagalli committed
338
    virtual void setup ();
339

340
341
342
343
344
345
346
347
348
349
350
351
    //! Solve the global hybrid system.
    void solve ();

    //! Compute primal and dual variables from the hybrid variable as a post process.
    void computePrimalAndDual ();

    //@}

    //! @name Set Methods
    //@{

    //! Set the boundary conditions.
352
353
354
    /*!
      @param bcHandler Boundary condition handler for the problem.
    */
355
    void setBC ( bchandler_raw_Type& bcHandler )
356
357
358
359
360
    {
        M_BCh   = &bcHandler;
        M_setBC = true;
    }

361
    //! Set source term
362
    /*!
363
      Set the source term, the default source term is the zero function.
364
365
      @param source Source term for the problem.
    */
366
    void setSource ( const function_Type& source )
367
368
369
370
    {
        M_source = source;
    }

371
372
373
374
375
376
377
378
379
    //! Set vector source term
    /*!
      @param source Vector source term for the problem.
    */
    void setVectorSource ( const vectorFunction_Type& vectorSource )
    {
        M_vectorSource = vectorSource;
    }

380
    //! Set inverse diffusion tensor
381
    /*!
382
      Set the inverse of diffusion tensor, the default inverse of permeability is the identity matrix.
383
      @param invPerm Inverse of the permeability tensor for the problem.
384
    */
385
    void setInversePermeability ( const permeability_Type& invPerm )
386
    {
387
        M_inversePermeability.reset ( new permeability_Type( invPerm ) );
388
389
    }

390
    //! Set the hybrid solution vector.
391
    /*!
392
      @param hybrid Constant vectorPtr_Type reference of the hybrid vector.
393
     */
394
    void setHybridSolution ( const vectorPtr_Type& hybrid  )
395
396
397
398
    {
        M_hybrid = hybrid;
    }

399
    //! Set the primal solution vector.
400
    /*!
401
      @param primal Constant vector_Type reference of the primal solution.
402
    */
403
    void setPrimalSolution ( const vectorPtr_Type& primal )
404
405
406
407
    {
        M_primal = primal;
    }

408
    //! Set the dual solution vector.
409
    /*!
410
      @param dual Constant vectorPtr_Type reference of the dual solution.
411
    */
412
    void setDualSolution ( const vectorPtr_Type& dual )
413
414
415
416
    {
        M_dual = dual;
    }

417
418
    //@}

419
    //! @name Get Methods
420
421
    //@{

422
    //! Returns pointer to the hybrid solution vector.
423
    /*!
424
      @return Constant vectorPtr_Type reference of the hybrid vector.
425
     */
426
    const vectorPtr_Type& hybridSolution () const
427
    {
428
        return M_hybrid;
429
430
    }

Michel Kern's avatar
Michel Kern committed
431
432
433
434
435
          vectorPtr_Type& hybridSolution ()
    {
        return M_hybrid;
    }

436
    //! Returns pointer to the primal solution vector.
437
    /*!
438
      @return Constant vector_Type reference of the primal solution.
439
    */
440
    const vectorPtr_Type& primalSolution () const
441
    {
442
        return M_primal;
443
444
    }

Michel Kern's avatar
Michel Kern committed
445
446
447
448
449
          vectorPtr_Type& primalSolution ()
    {
        return M_primal;
    }

450
    //! Returns pointer to the dual solution vector.
451
    /*!
452
      @return Constant vectorPtr_Type reference of the dual solution.
453
    */
454
    const vectorPtr_Type& dualSolution () const
455
    {
456
        return M_dual;
457
458
    }

Michel Kern's avatar
Michel Kern committed
459
460
461
462
463
          vectorPtr_Type& dualSolution ()
    {
        return M_dual;
    }

464
    /*!
465
      Returns the pointer of the residual vector.
466
      @return Constant vectorPtr_Type reference of the residual.
467
    */
468
    const vectorPtr_Type& residual () const
469
    {
470
        return M_residual;
471
472
    }

Michel Kern's avatar
Michel Kern committed
473
474
475
476
477
          vectorPtr_Type& residual ()
    {
        return M_residual;
    }

478
    //! Returns whether bounday conditions is set
479
480
481
482
    /*!
      @return Constant boolean with value true if the boundary condition is setted,
      false otherwise
    */
Michel Kern's avatar
Michel Kern committed
483
    bool isBCset ()
484
485
486
487
    {
        return M_setBC;
    }

488
    //! Returns boundary conditions handler.
489
490
491
    /*!
      @return Reference of boundary conditions handler.
    */
Michel Kern's avatar
Michel Kern committed
492
493
494
495
496
497
    const bchandler_Type& bcHandler () const
    {
        return M_BCh;
    }

          bchandler_Type& bcHandler ()
498
499
500
501
    {
        return M_BCh;
    }

502
    //! Returns Epetra local map.
503
    /*!
504
      @return Constant MapEpetra reference of the problem.
505
    */
506
    const MapEpetra& getMap () const
Michel Kern's avatar
Michel Kern committed
507
508
509
510
    {
        return M_localMap;
    }

511
          MapEpetra& getMap ()
512
513
514
515
    {
        return M_localMap;
    }

516
    //! Returns Epetra communicator.
517
    /*!
518
      @return Constant reference of the shared pointer of the communicator of the problem.
519
    */
Michel Kern's avatar
Michel Kern committed
520
521
522
523
524
525
    const commPtr_Type & getComm () const
    {
        return M_displayer.comm();
    }

          commPtr_Type & getComm ()
526
    {
527
        return M_displayer.comm();
528
529
    }

530
    //! Returns displayer.
531
    /*!
532
      @return Constant reference of the displayer of the problem.
533
    */
Michel Kern's avatar
Michel Kern committed
534
535
536
537
538
539
    const Displayer & getDisplayer() const
    {
        return M_displayer;
    }

          Displayer & getDisplayer()
540
    {
541
        return M_displayer;
542
543
    }

544
545
    //@}

546
547
protected:

548
549
    //! @name Protected Methods
    //@{
550

fumagalli's avatar
Add    
fumagalli committed
551
    //! Pre-computes local (element independant) matrices
552
    /*!
553
554
      Compute all the local matrices that are independant
      from the geometrical element.
555
    */
556
    virtual void computeConstantMatrices ();
557

fumagalli's avatar
Add    
fumagalli committed
558
    //! Computes local (element dependent) matrices
559
560
561
562
563
564
    /*!
      Locally update the current finite element for the primal
      and dual finite element space, then compute the Hdiv mass
      matrix.
      @param iElem Id of the current geometrical element.
    */
fumagalli's avatar
fumagalli committed
565
    virtual void localElementComputation ( const UInt & iElem );
566

567
    //! Performs static condensation
568
    /*!
569
      Locally eliminate pressure and velocity DOFs, create the local
570
571
      hybrid matrix and local hybrid right hand side.
    */
fumagalli's avatar
fumagalli committed
572
    virtual void staticCondensation ();
573
574

    //! Compute locally, as a post process, the primal and dual variable.
fumagalli's avatar
fumagalli committed
575
    virtual void localComputePrimalAndDual ();
576

577
    //! Update all problem variables
578
579
580
581
582
583
    /*!
      Update all the variables of the problem before the construction of
      the global hybrid matrix, e.g. reset the global hybrid matrix.
      It is principally used for a time dependent derived class, in fact we
      want to clear each time step all the variables.
    */
fumagalli's avatar
fumagalli committed
584
    virtual void updateVariables ();
585

fumagalli's avatar
Add    
fumagalli committed
586
    //! Apply the boundary condition
587
588
589
590
591
    /* Apply BC to the hybrid global matrix and to the hybrid global right hand side.
    */
    void applyBoundaryConditions ();

    //! Make matrix symmetric
592
593
594
595
596
597
    /*!
      Transform a symmetric matrix that is stored only in the lower or upper
      triangular part in a full symmetric matrix.
      @param UPLO Character that indicate if A is in the upper triangular part
      using flag 'U' or in the lower triangular part using flag 'L'.
      @param N The size of the matrix A.
598
      @param A The matrix to be reordered.
599
600
    */
    template<typename matrix>
601
    void symmetrizeMatrix ( char UPLO, Int N, matrix& A  );
602

603
604
    //@}

605
606
607
608
    //! @name Parallel stuff
    //@{

    //! MPI process identifier.
609
    UInt      M_me;
610
611

    //! Local map.
612
    MapEpetra M_localMap;
613
614
615

    //! Displayer.
    Displayer M_displayer;
616
617
618
619
620
621
622
623

    //@}

    // Data of the problem.
    //! @name Data of the problem
    //@{

    //! Data for Darcy solvers.
624
    const data_Type&     M_data;
625
626

    //! Source function.
627
628
629
630
    function_Type        M_source;

    //! Vector source function.
    vectorFunction_Type  M_vectorSource;
631
632

    //! Permeability tensor.
633
    permeabilityPtr_Type M_inversePermeability;
634
635

    //! Bondary conditions handler.
636
    bchandler_raw_Type*  M_BCh;
637
638
639
640
641
642
643
644
645
646

    //! Flag if the boundary conditions are setted or not.
    bool                 M_setBC;

    //@}

    // Finite element spaces.
    //! @name Finite element spaces
    //@{

647
    //! Primal finite element space.
648
    FESpace<Mesh, MapEpetra>&  M_primal_FESpace;
649
650

    //! Dual finite element space.
651
    FESpace<Mesh, MapEpetra>&  M_dual_FESpace;
652
653

    //! Hybrid finite element space.
654
    FESpace<Mesh, MapEpetra>&  M_hybrid_FESpace;
655
656

    //! Dual function dot outward unit vector finite element space.
657
    FESpace<Mesh, MapEpetra>&  M_VdotN_FESpace;
658
659
660
661
662
663
664
665

    //@}

    // Algebraic stuff.
    //! @name Algebraic stuff
    //@{

    //! Hybrid matrix.
666
    matrixPtr_Type                          M_matrHybrid;
667
668

    //! Right hand side.
669
    vectorPtr_Type                          M_rhs;
670
671

    //! Primal solution.
672
    vectorPtr_Type                          M_primal;
673
674

    //! Dual solution.
675
    vectorPtr_Type                          M_dual;
676
677

    //! Hybrid solution.
678
    vectorPtr_Type                          M_hybrid;
679
680

    //! Residual.
681
    vectorPtr_Type                          M_residual;
682
683

    //! Linear solver.
Michel Kern's avatar
Michel Kern committed
684
    solver_Type                             M_linearSolver;
685
686

    //! Epetra preconditioner for the linear system.
687
    boost::shared_ptr<Preconditioner> M_prec;
688
689
690
691
692
693
694
695

    //@}

    // Elementary matrices and vectors used for the static condensation.
    //! @name Elementary matrices and vectors used for the static condensation.
    //@{

    //! Element vector for the right hand side, lives in RTkHyb(fe).
696
    VectorElemental M_elvecHyb;
697
698

    //! Element vector for the source term of the problem, lives in Qk(fe).
699
    VectorElemental M_elvecSource;
700

701
702
703
    //! Element vector for the vector source term of the problem, lives in RTk(fe)
    VectorElemental M_elvecSourceVector;

704
    //! Element vector for the dual variabl, lives in RTk(fe).
705
    VectorElemental M_elvecFlux;
706
707
708
709
710
711
712
713
714
715
716
717
718

    /*! Mixed hybrid matrix, maps
      \f[
      [A\,| B\,| C] \qquad \mathrm{in} \qquad
      \left[
      \begin{array}{c c c}
      A   & B & C \\
      B^T & 0 & 0 \\
      C^T & 0 & 0 \\
      \end{array}
      \right]
      \f]
    */
719
    MatrixElemental M_elmatMix;
720

fumagalli's avatar
fumagalli committed
721
    //! Element hybrid matrix.
722
    MatrixElemental M_elmatHyb;
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742

    //! Temporary array of size the square of the number of degrees of freedom of the primal variable.
    KNM<Real> M_BtB;

    /*! Temporary array of size the square of the number of degrees of freedom of the hybrid variable.
        It stores of the final hybrid array. */
    KNM<Real> M_CtC;

    /*! Temporary array of size the prodocut between the number of degrees of freedom of the primal
        and the degrees of freedom of the hybrid variable. */
    KNM<Real> M_BtC;

    //@}

}; // class DarcySolver

//
// IMPLEMENTATION
//

743
744
745
746
// ====================================================================================================
// Constructors and destructors
// ====================================================================================================

747
748
749
// Complete constructor.
template<typename Mesh, typename SolverType>
DarcySolver<Mesh, SolverType>::
Michel Kern's avatar
Michel Kern committed
750
DarcySolver ( const data_Type&                 dataFile,
751
752
753
754
              FESpace<Mesh, MapEpetra>&  primal_FESpace,
              FESpace<Mesh, MapEpetra>&  dual_FESpace,
              FESpace<Mesh, MapEpetra>&  hybrid_FESpace,
              FESpace<Mesh, MapEpetra>&  VdotN_FESpace,
Michel Kern's avatar
Michel Kern committed
755
756
              const bchandler_raw_Type&        bcHandler,
              const commPtr_Type&              comm ):
Tiziano Passerini's avatar
Tiziano Passerini committed
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
        // Parallel stuff.
        M_me                     ( comm->MyPID() ),
        M_localMap               ( hybrid_FESpace.map() ),
        M_displayer              ( comm ),
        // Data of the problem.
        M_data                   ( dataFile ),
        M_source                 ( DarcyDefaultSource() ),
        M_BCh                    ( &bcHandler ),
        M_setBC                  ( true ),
        // Finite element spaces.
        M_primal_FESpace         ( primal_FESpace ),
        M_dual_FESpace           ( dual_FESpace ),
        M_hybrid_FESpace         ( hybrid_FESpace ),
        M_VdotN_FESpace          ( VdotN_FESpace ),
        // Algebraic stuff.
772
773
        M_matrHybrid             ( new matrix_Type ( M_localMap ) ),
        M_rhs                    ( new vector_Type ( M_localMap ) ),
fumagalli's avatar
fumagalli committed
774
775
        M_primal             ( new vector_Type ( M_primal_FESpace.map() ) ),
        M_dual           ( new vector_Type ( M_dual_FESpace.map(), Repeated ) ),
776
777
        M_hybrid                 ( new vector_Type ( M_hybrid_FESpace.map() ) ),
        M_residual               ( new vector_Type ( M_localMap ) ),
Tiziano Passerini's avatar
Tiziano Passerini committed
778
779
780
781
782
        M_linearSolver           ( ),
        M_prec                   ( ),
        // Elementary matrices and vectors used for the static condensation.
        M_elvecHyb               ( M_hybrid_FESpace.refFE().nbDof(), 1 ),
        M_elvecSource            ( M_primal_FESpace.refFE().nbDof(), 1 ),
783
        M_elvecSourceVector      ( M_dual_FESpace.refFE().nbDof(), 1 ),
Tiziano Passerini's avatar
Tiziano Passerini committed
784
785
786
787
788
789
790
        M_elvecFlux              ( M_dual_FESpace.refFE().nbDof(), 1 ),
        M_elmatMix               ( M_dual_FESpace.refFE().nbDof(), 1, 1, M_primal_FESpace.refFE().nbDof(), 0, 1,
                                   M_hybrid_FESpace.refFE().nbDof(), 0, 1 ),
        M_elmatHyb               ( M_hybrid_FESpace.refFE().nbDof(), 1, 1 ),
        M_BtB                    ( M_primal_FESpace.refFE().nbDof(), M_primal_FESpace.refFE().nbDof() ),
        M_CtC                    ( M_hybrid_FESpace.refFE().nbDof(), M_hybrid_FESpace.refFE().nbDof() ),
        M_BtC                    ( M_primal_FESpace.refFE().nbDof(), M_hybrid_FESpace.refFE().nbDof() )
791
792
793
794
795
796
797
798
{

} // Constructor


// Constructor without boundary condition handler.
template<typename Mesh, typename SolverType>
DarcySolver<Mesh, SolverType>::
Michel Kern's avatar
Michel Kern committed
799
DarcySolver ( const data_Type&                 dataFile,
800
801
802
803
              FESpace<Mesh, MapEpetra>&  primal_FESpace,
              FESpace<Mesh, MapEpetra>&  dual_FESpace,
              FESpace<Mesh, MapEpetra>&  hybrid_FESpace,
              FESpace<Mesh, MapEpetra>&  VdotN_FESpace,
Michel Kern's avatar
Michel Kern committed
804
              const commPtr_Type&              comm ):
Tiziano Passerini's avatar
Tiziano Passerini committed
805
806
807
808
809
810
811
812
813
814
815
816
817
818
        // Parallel stuff.
        M_me                     ( comm->MyPID() ),
        M_localMap               ( hybrid_FESpace.map() ),
        M_displayer              ( comm ),
        // Data of the problem.
        M_data                   ( dataFile ),
        M_source                 ( DarcyDefaultSource() ),
        M_setBC                  ( false ),
        // Finite element spaces.
        M_primal_FESpace         ( primal_FESpace ),
        M_dual_FESpace           ( dual_FESpace ),
        M_hybrid_FESpace         ( hybrid_FESpace ),
        M_VdotN_FESpace          ( VdotN_FESpace ),
        // Algebraic stuff.
819
820
        M_matrHybrid             ( new matrix_Type ( M_localMap ) ),
        M_rhs                    ( new vector_Type ( M_localMap ) ),
fumagalli's avatar
fumagalli committed
821
822
        M_primal                 ( new vector_Type ( M_primal_FESpace.map() ) ),
        M_dual                   ( new vector_Type ( M_dual_FESpace.map(), Repeated ) ),
823
824
        M_hybrid                 ( new vector_Type ( M_hybrid_FESpace.map() ) ),
        M_residual               ( new vector_Type ( M_localMap ) ),
Tiziano Passerini's avatar
Tiziano Passerini committed
825
826
827
828
829
        M_linearSolver           ( ),
        M_prec                   ( ),
        // Local matrices and vectors.
        M_elvecHyb               ( M_hybrid_FESpace.refFE().nbDof(), 1 ),
        M_elvecSource            ( M_primal_FESpace.refFE().nbDof(), 1 ),
830
        M_elvecSourceVector      ( M_dual_FESpace.refFE().nbDof(), 1 ),
Tiziano Passerini's avatar
Tiziano Passerini committed
831
832
833
834
835
836
837
        M_elvecFlux              ( M_dual_FESpace.refFE().nbDof(), 1 ),
        M_elmatMix               ( M_dual_FESpace.refFE().nbDof(), 1, 1, M_primal_FESpace.refFE().nbDof(), 0, 1,
                                   M_hybrid_FESpace.refFE().nbDof(), 0, 1 ),
        M_elmatHyb               ( M_hybrid_FESpace.refFE().nbDof(), 1, 1 ),
        M_BtB                    ( M_primal_FESpace.refFE().nbDof(), M_primal_FESpace.refFE().nbDof() ),
        M_CtC                    ( M_hybrid_FESpace.refFE().nbDof(), M_hybrid_FESpace.refFE().nbDof() ),
        M_BtC                    ( M_primal_FESpace.refFE().nbDof(), M_hybrid_FESpace.refFE().nbDof() )
838
839
840
841
842
843
844
{

} // Constructor

// Virtual destructor.
template<typename Mesh, typename SolverType>
DarcySolver<Mesh, SolverType>::
fumagalli's avatar
fumagalli committed
845
~DarcySolver ( void )
846
847
848
849
{

} // Destructor

850
851
852
853
854
855
856
857
858
859
860
// ===========================================================================================
// Public methods
// ==========================================================================================

// Build the global hybrid matrix and global hybrid right hand side, with boundary conditions.
template<typename Mesh, typename SolverType>
void
DarcySolver<Mesh, SolverType>::
buildSystem ()
{

861
862
863
864
865
    // LifeChronos.
    LifeChrono chronoStaticCondensation;
    LifeChrono chronoConstantLocalMatrix;
    LifeChrono chronoGlobalAssembleMatrix;
    LifeChrono chronoGlobalAssembleVector;
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907

    // The total number of elements in the mesh.
    UInt meshNumberOfElements = M_primal_FESpace.mesh()->numElements();

    M_displayer.leaderPrint("  Darcy solver - Computing constant matrices   ...          ");

    // Clear all the local matrices and vectors.
    M_elvecHyb.zero();
    M_elvecSource.zero();
    M_elvecFlux.zero();
    M_elmatMix.zero();
    M_elmatHyb.zero();

    chronoConstantLocalMatrix.start();

    // Compute all the constant matrices, e.g. the matrix B and C
    computeConstantMatrices();

    chronoConstantLocalMatrix.stop();

    M_displayer.leaderPrintMax( "done in " , chronoConstantLocalMatrix.diff() );

    // If setted print the constant matrices computed.
    if ( M_displayer.isLeader() &&  M_data.verbose() > static_cast<UInt>(1) )
    {
        M_displayer.leaderPrint( "elmatHyb :\n\n" );
        M_elmatHyb.showMe();

        M_displayer.leaderPrint( "elmatMix :\n\n" );
        M_elmatMix.showMe();
    }

    // Update all the variables, e.g. reset the hybrid matrix
    updateVariables();

    M_displayer.leaderPrint("                 Perform Static Condensation   ...          ");
    chronoStaticCondensation.start();

    //---------------------------------------
    //    LOOP ON ALL THE VOLUME ELEMENTS
    //---------------------------------------

908
    for ( UInt iElem(0); iElem < meshNumberOfElements; ++iElem )
909
910
911
912
913
914
915
916
    {
        // Compute the Hdiv mass matrix as a local matrix depending on the current element.
        localElementComputation( iElem );

        staticCondensation();

        /* Assemble the global hybrid matrix.
           M_primal_FESpace is used instead of M_hybrid_FESpace for currentLocalId,
917
           because currentFE cannot store a ReferenceFEHybrid. */
918
919
920
921
922
923
924
925
926
        assembleMatrix( *M_matrHybrid,
                        M_primal_FESpace.fe().currentLocalId(),
                        M_elmatHyb,
                        M_hybrid_FESpace.refFE().nbDof(),
                        M_hybrid_FESpace.dof(),
                        0,0,0,0);

        /* Assemble the global hybrid right hand side.
           M_primal_FESpace is used instead of M_hybrid_FESpace for currentLocalId,
927
           because currentFE cannot store a ReferenceFEHybrid. */
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
        assembleVector( *M_rhs,
                        M_primal_FESpace.fe().currentLocalId(),
                        M_elvecHyb,
                        M_hybrid_FESpace.refFE().nbDof(),
                        M_hybrid_FESpace.dof(), 0 );
    }

    //-----------------------------------
    // END OF LOOP ON THE VOLUME ELEMENTS
    //-----------------------------------

    chronoStaticCondensation.stop();
    M_displayer.leaderPrintMax( "done in " , chronoStaticCondensation.diff() );

    // Apply boundary conditions to the hybrid global matrix and to the hybrid global right hand side.
    applyBoundaryConditions();

    M_displayer.leaderPrint("                 Assemble Global Hybrid Matrix ...          ");

    chronoGlobalAssembleMatrix.start();

    // Assemble the global hybrid matrix.
950
    M_matrHybrid->globalAssemble();
951
952
953
954
955
956
957
958
959
960

    chronoGlobalAssembleMatrix.stop();

    M_displayer.leaderPrintMax( "done in " , chronoGlobalAssembleMatrix.diff() );

    M_displayer.leaderPrint("                 Assemble Global Hybrid Vector ...          ");

    chronoGlobalAssembleVector.start();

    // Assemble the global hybrid vector.
961
    M_rhs->globalAssemble();
962
963
964
965
966
967
968

    chronoGlobalAssembleVector.stop();

    M_displayer.leaderPrintMax( "done in " , chronoGlobalAssembleVector.diff() );

} // buildSystem

969
// Set up the linear solver and the preconditioner.
970
971
972
template<typename Mesh, typename SolverType>
void
DarcySolver<Mesh, SolverType>::
fumagalli's avatar
fumagalli committed
973
setup ()
974
{
fumagalli's avatar
fumagalli committed
975

976
977
978
    GetPot dataFile( *(M_data.dataFile()) );

    // Set up data for the linear solver and the preconditioner.
979
980
    M_linearSolver.setDataFromGetPot( dataFile, ( M_data.section() + "/solver" ).data() );
    M_linearSolver.setupPreconditioner( dataFile, ( M_data.section() + "/prec" ).data() );
fumagalli's avatar
Add    
fumagalli committed
981
    M_linearSolver.setCommunicator( M_displayer.comm() );
982
983

    // Choose the preconditioner type.
984
    std::string precType = dataFile( ( M_data.section() + "/prec/prectype" ).data() , "Ifpack");
985

Tiziano Passerini's avatar
Tiziano Passerini committed
986
    // Create a preconditioner object.
987
988
989
    M_prec.reset( PRECFactory::instance().createObject( precType ) );
    ASSERT( M_prec.get() != 0, "DarcySolver : Preconditioner not set" );

perego's avatar
perego committed
990
991
992
993
    // make sure mesh facets are updated
    if(! M_dual_FESpace.mesh()->hasLocalFacets() )
    M_dual_FESpace.mesh()->updateElementFacets();

994
995
} // setup

996
// Solve the linear system.
997
998
999
template<typename Mesh, typename SolverType>
void
DarcySolver<Mesh, SolverType>::
1000
solve ()