TimeIterationPolicyLinear.hpp 6.46 KB
Newer Older
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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
//@HEADER
/*
*******************************************************************************

    Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
    Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University

    This file is part of LifeV.

    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.

    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.

    You should have received a copy of the GNU Lesser General Public License
    along with LifeV.  If not, see <http://www.gnu.org/licenses/>.

*******************************************************************************
*/
//@HEADER

/*!
    @file TimeIterationLinear class
    @brief This class contains all the informations necessary
           to perform a time iteration

    @author Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
    @date 06-12-2012
 */

#ifndef TIMEITERATIONPOLICYLINEAR_HPP
#define TIMEITERATIONPOLICYLINEAR_HPP

#include <iostream>
#include <string>
#include <boost/shared_ptr.hpp>

// Tell the compiler to ignore specific kind of warnings:
#pragma GCC diagnostic ignored "-Wunused-variable"
#pragma GCC diagnostic ignored "-Wunused-parameter"

#ifdef HAVE_MPI
#include <Epetra_MpiComm.h>
#else
#include <Epetra_SerialComm.h>
#endif

#include <Teuchos_ParameterList.hpp>
#include <Teuchos_XMLParameterListHelpers.hpp>
#include <Teuchos_RCP.hpp>

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

#include <lifev/core/LifeV.hpp>
#include <lifev/core/array/MatrixEpetra.hpp>
#include <lifev/core/array/VectorEpetra.hpp>
#include <lifev/core/util/Displayer.hpp>
#include <lifev/core/util/LifeChrono.hpp>
#include <lifev/core/mesh/RegionMesh.hpp>
#include <lifev/core/fem/FESpace.hpp>
#include <lifev/core/fem/TimeAdvanceBDF.hpp>
#include <lifev/core/fem/BCHandler.hpp>
#include <lifev/navier_stokes/solver/NavierStokesSolver/SolverPolicyLinearSolver.hpp>


namespace LifeV
{

grandper's avatar
grandper committed
76
template< class mesh_Type, class AssemblyPolicy, class SolverPolicy = SolverPolicyLinearSolver >
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
struct TimeIterationPolicyLinear : private AssemblyPolicy, public virtual SolverPolicy
{
protected:
    typedef MatrixEpetra<Real>                       matrix_Type;
    typedef boost::shared_ptr<matrix_Type>           matrixPtr_Type;
    typedef VectorEpetra                             vector_Type;
    typedef boost::shared_ptr<VectorEpetra>          vectorPtr_Type;
    typedef MeshPartitioner< mesh_Type >             meshPartitioner_Type;
    typedef MapEpetra                                map_Type;
    typedef boost::shared_ptr<map_Type>              mapPtr_Type;
    typedef FESpace< mesh_Type, map_Type >           fespace_Type;
    typedef boost::shared_ptr< fespace_Type >        fespacePtr_Type;
    typedef TimeAdvanceBDF<vector_Type>              bdf_Type;
    typedef boost::shared_ptr< bdf_Type >            bdfPtr_Type;
    typedef BCHandler                                bcContainer_Type;
    typedef boost::shared_ptr<bcContainer_Type>      bcContainerPtr_Type;

    enum { BDFOrder = AssemblyPolicy::BDFOrder };

    void initTimeIteration ( Teuchos::ParameterList& list );

    void iterate ( vectorPtr_Type solution,
                   bcContainerPtr_Type bchandler,
                   const Real& currentTime );

    // Parameters
    bool                    M_computeResidual;

    matrixPtr_Type          M_systemMatrix;
    vectorPtr_Type          M_rhs;
    mapPtr_Type             M_solutionMap;

    virtual Displayer displayer() = 0;
    virtual fespacePtr_Type uFESpace() const = 0;
    virtual fespacePtr_Type pFESpace() const = 0;
};

grandper's avatar
grandper committed
114
template< class mesh_Type, class AssemblyPolicy, class SolverPolicy >
115
void
grandper's avatar
grandper committed
116
TimeIterationPolicyLinear<mesh_Type, AssemblyPolicy, SolverPolicy>::
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
initTimeIteration ( Teuchos::ParameterList& list )
{
    // Parameters
    M_computeResidual = list.get ( "Compute exact residual", false );

    // Initialization
    M_solutionMap.reset ( new map_Type ( uFESpace()->map() + pFESpace()->map() ) );

    // Init the assembler
    Teuchos::ParameterList assemblyList = list.sublist ( "Assembly: Parameter list" );
    AssemblyPolicy::initAssembly ( assemblyList );

    // Init the solver
    Teuchos::ParameterList solverList = list.sublist ( "Solver: Parameter list" );
    SolverPolicy::initSolver ( solverList );
    M_rhs.reset ( new vector_Type ( *M_solutionMap, Unique ) );
}

grandper's avatar
grandper committed
135
template< class mesh_Type, class AssemblyPolicy, class SolverPolicy >
136
void
grandper's avatar
grandper committed
137
TimeIterationPolicyLinear<mesh_Type, AssemblyPolicy, SolverPolicy>::
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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
iterate ( vectorPtr_Type solution,
          bcContainerPtr_Type bchandler,
          const Real& currentTime )
{
    Real rhsIterNorm ( 0.0 );

    //
    // STEP 1: Updating the system
    //
    displayer().leaderPrint ( "Updating the system... " );
    *M_rhs = 0.0;
    M_systemMatrix.reset ( new matrix_Type ( *M_solutionMap ) );
    AssemblyPolicy::assembleSystem ( M_systemMatrix, M_rhs, solution, SolverPolicy::preconditioner() );
    displayer().leaderPrint ( "done\n" );

    //
    // STEP 2: Applying the boundary conditions
    //
    displayer().leaderPrint ( "Applying BC... " );
    bcManage ( *M_systemMatrix, *M_rhs, *uFESpace()->mesh(), uFESpace()->dof(), *bchandler, uFESpace()->feBd(), 1.0, currentTime );
    M_systemMatrix->globalAssemble();
    displayer().leaderPrint ( "done\n" );

    // Extra information if we want to know the exact residual
    if ( M_computeResidual )
    {
        rhsIterNorm = M_rhs->norm2();
    }

    //
    // STEP 3: Solving the system
    //
    displayer().leaderPrint ( "Solving the system... \n" );
    SolverPolicy::solve ( M_systemMatrix, M_rhs, solution );

    if ( M_computeResidual )
    {
        vector_Type Ax ( solution->map() );
        vector_Type res ( *M_rhs );
        M_systemMatrix->matrixPtr()->Apply ( solution->epetraVector(), Ax.epetraVector() );
        res.epetraVector().Update ( -1, Ax.epetraVector(), 1 );
        Real residual;
        res.norm2 ( &residual );
        residual /= rhsIterNorm;
        displayer().leaderPrint ( "Scaled residual: ", residual, "\n" );
    }
}

} // namespace LifeV

#endif /* TIMEITERATIONPOLICYLINEAR_HPP */