ZeroDimensionalRythmosModelInterface.cpp 8.9 KB
Newer Older
1
2
//@HEADER
/*
Cristiano Malossi's avatar
Cristiano Malossi committed
3
*******************************************************************************
4

Cristiano Malossi's avatar
Cristiano Malossi committed
5
6
    Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
    Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7

Cristiano Malossi's avatar
Cristiano Malossi committed
8
    This file is part of LifeV.
9

Cristiano Malossi's avatar
Cristiano Malossi committed
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

Cristiano Malossi's avatar
Cristiano Malossi committed
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

Cristiano Malossi's avatar
Cristiano Malossi committed
20
21
    You should have received a copy of the GNU Lesser General Public License
    along with LifeV.  If not, see <http://www.gnu.org/licenses/>.
22

Cristiano Malossi's avatar
Cristiano Malossi committed
23
24
*******************************************************************************
*/
25
26
27
28
//@HEADER

/*!
 *  @file
Cristiano Malossi's avatar
Cristiano Malossi committed
29
 *  @brief Rythmos Model Interface
30
31
32
33
 *
 *  @date 21-11-2011
 *  @author Mahmoud Jafargholi
 *
34
35
 *  @contributors Cristiano Malossi <cristiano.malossi@epfl.ch>
 *  @mantainer    Cristiano Malossi <cristiano.malossi@epfl.ch>
36
37
 */

38
#include <lifev/zero_dimensional/solver/ZeroDimensionalRythmosModelInterface.hpp>
39

40
namespace LifeV {
41

42
#if ( defined(HAVE_NOX_THYRA) && defined(HAVE_TRILINOS_RYTHMOS) )
43
44
45
// ===================================================
// Constructors
// ===================================================
46
RythmosModelInterface::RythmosModelInterface( Int numCircuitElements,
47
48
                                              Epetra_Comm* comm,
                                              zeroDimensionalCircuitDataPtr_Type circuitData ) :
49
    M_numCircuitElements( numCircuitElements ), M_numMyElements( 0 ), M_myPID( comm->MyPID() ), M_numProc( comm->NumProc() ),
50
51
52
53
54
55

    M_standardMap( 0 ), M_initialSolutionY( 0 ), M_initialSolutionYp( 0 ), M_circuitData( circuitData )
{
    M_comm = comm->Clone();
    ( *M_comm ).PrintInfo( std::cout );
    M_commSharedPtr.reset( M_comm );
56
    M_mapEpetraPtr.reset( new MapEpetra( M_numCircuitElements,
57
58
59
60
                                         M_commSharedPtr ) );
    M_standardMap = ( M_mapEpetraPtr->map( Unique ) ).get();
    M_numMyElements = M_standardMap->NumMyElements();
    M_A.reset( new matrix_Type( *M_mapEpetraPtr,
61
                                M_numCircuitElements ) );
62
    M_B.reset( new matrix_Type( *M_mapEpetraPtr,
63
                                M_numCircuitElements ) );
64

65
    for ( Int i = 0; i < M_numCircuitElements; i++ )
66
    {
67
        for ( Int j = 0; j < M_numCircuitElements; j++ )
68
69
70
71
        {
            M_A->addToCoefficient( i, j, 1.0 );
            M_B->addToCoefficient( i, j, 1.0 );
        }
72
    }
73
74
    M_A->globalAssemble();
    M_B->globalAssemble();
75

76
77
78
    M_graph = new Epetra_CrsGraph( M_A->matrixPtr()->Graph() );
    M_graphSharedPtr.reset( M_graph );
    M_C.reset( new vector_Type( *M_mapEpetraPtr ) );
79

80
81
    M_fA.reset( new vectorEpetra_Type( *M_standardMap ) );
    M_fB.reset( new vectorEpetra_Type( *M_standardMap ) );
82

83
84
    M_Y0.reset( new vectorEpetra_Type( *M_standardMap ) );
    M_Yp0.reset( new vectorEpetra_Type( *M_standardMap ) );
85

86
87
88
89
90
    M_initialSolutionY = M_Y0.get();
    M_initialSolutionYp = M_Yp0.get();
    ;
    initializeSolnY();
    initializeSolnYp();
91
92
93
94
95
96
}

// Destructor
RythmosModelInterface::~RythmosModelInterface()
{

97
98
99
100
101
102
103
104
    M_circuitData.reset();
    M_mapEpetraPtr.reset();
    M_graphSharedPtr.reset();
    M_A.reset();
    M_B.reset();
    M_C.reset();
    M_Y0.reset();
    M_Yp0.reset();
105
106
}

107
108
109
bool RythmosModelInterface::computeF( __attribute__((unused))  const Epetra_Vector& x,
                                      __attribute__((unused))  Epetra_Vector& FVec,
                                      __attribute__((unused))  NOX::Epetra::Interface::Required::FillType fillType )
110
{
111
112
    std::cout << "Warning: I should not be here, 0D, Rythmos Model Interface";
    return false;
113
114
}

115
116
bool RythmosModelInterface::computeJacobian( __attribute__((unused))  const Epetra_Vector& x,
                                             __attribute__((unused))  Epetra_Operator& my_Jac )
117
{
118
119
    std::cout << "Warning: I should not be here, 0D, Rythmos Model Interface";
    return false;
120
121
}

122
bool RythmosModelInterface::computePrecMatrix( __attribute__((unused))  const Epetra_Vector& x )
123
{
124
125
    std::cout << "Warning: I should not be here, 0D, Rythmos Model Interface";
    return false;
126
}
127
128
129
bool RythmosModelInterface::computePreconditioner( __attribute__((unused))   const Epetra_Vector& x,
                                                   __attribute__((unused))   Epetra_Operator& my_Prec,
                                                   __attribute__((unused))   Teuchos::ParameterList* precParams )
130
{
131
132
133
    cout << "ERROR: Interface::preconditionVector()  " << endl;
    std::cout << "Error: I should not be here, 0D, Rythmos Model Interface";
    throw "Interface Error";
134
135
}

136
137
138
bool RythmosModelInterface::evaluate( __attribute__((unused))   Real t,
                                      __attribute__((unused))  const Epetra_Vector* x,
                                      __attribute__((unused))   Epetra_Vector* f )
139
{
140
141
    std::cout << "Warning: I should not be here. 0D, Rythmos Inreface, ::evaluate" << std::endl;
    return true;
142
143
144
145
}

Epetra_Vector& RythmosModelInterface::getSolutionY()
{
146
    return *M_initialSolutionY;
147
148
149
150
}

Epetra_Vector& RythmosModelInterface::getSolutionYp()
{
151
    return *M_initialSolutionYp;
152
153
154
155
}

Epetra_Map& RythmosModelInterface::getMap()
{
156
    return *M_standardMap;
157
158
159
160
}

bool RythmosModelInterface::initializeSolnY()
{
161
162
    M_initialSolutionY->PutScalar( 0 );
    return true;
163
164
}

165
bool RythmosModelInterface::initializeSolnY( const vectorEpetra_Type& y )
166
{
167
168
    ( *M_initialSolutionY ) = y;
    return true;
169
170
}

171
bool RythmosModelInterface::initializeSolnYp( const vectorEpetra_Type& yp )
172
{
173
174
    ( *M_initialSolutionYp ) = yp;
    return true;
175
176
177
178
}

bool RythmosModelInterface::initializeSolnYp()
{
179
180
    M_initialSolutionYp->PutScalar( 0 );
    return true;
181
182
183
}
Epetra_CrsGraph& RythmosModelInterface::getGraph()
{
184
    return *M_graph;
185
186
}

187
188
189
190
191
bool RythmosModelInterface::evaluateFImplicit( const Real& t,
                                               const Epetra_Vector* x,
                                               const Epetra_Vector* x_dot,
                                               Epetra_Vector* f )
{
192

193
    //updateCircuitData from Y and Yp
194
#ifdef HAVE_LIFEV_DEBUG
195
196
    x->Print(std::cout);
    x_dot->Print(std::cout);
197
#endif
198
199
200
201
202
203
204
    M_circuitData->updateCircuitDataFromY( t,
                                           x,
                                           x_dot );
    M_circuitData->updateABC( *M_A,
                              *M_B,
                              *M_C );
    //calculate the sum
205
#ifdef HAVE_LIFEV_DEBUG
206
207
208
    M_A->matrixPtr()->Print(std::cout);
    M_B->matrixPtr()->Print(std::cout);
    M_C->epetraVector().Print(std::cout);
209
210
#endif

211
212
213
    M_A->matrixPtr()->Multiply1( false,
                                 *x_dot,
                                 *M_fA );
214
#ifdef HAVE_LIFEV_DEBUG
215
    M_A->matrixPtr()->Print(std::cout);
216
217
#endif

218
219
220
    M_B->matrixPtr()->Multiply1( false,
                                 *x,
                                 *M_fB );
221
#ifdef HAVE_LIFEV_DEBUG
222
    M_B->matrixPtr()->Print(std::cout);
223
#endif
224
    for ( Int i = 0; i < M_numCircuitElements; i++ )
225
226
227
    {
        ( *f )[i] = ( *M_fA )[i] + ( *M_fB )[i] + ( *M_C )[i];
    }
228
#ifdef HAVE_LIFEV_DEBUG
229
    f->Print(std::cout);
230
#endif
231
    return true;
232
233
}

234
235
236
237
238
239
240
241
bool RythmosModelInterface::evaluateWImplicit( const Real& t,
                                               const Real& alpha,
                                               const Real& beta,
                                               const Epetra_Vector* x,
                                               const Epetra_Vector* x_dot,
                                               Epetra_CrsMatrix* W )
{
    //updateCircuitData from Y and Yp
242
#ifdef HAVE_LIFEV_DEBUG
243
244
    x->Print(std::cout);
    x_dot->Print(std::cout);
245
#endif
246
247
248
249
250
251
    M_circuitData->updateCircuitDataFromY( t,
                                           x,
                                           x_dot );
    M_circuitData->updateABC( *M_A,
                              *M_B,
                              *M_C );
252
#ifdef HAVE_LIFEV_DEBUG
253
254
255
    M_A->matrixPtr()->Print(std::cout);
    M_B->matrixPtr()->Print(std::cout);
    M_C->epetraVector().Print(std::cout);
256
#endif
257
258
259
260
261
    M_A->operator*=( alpha );
    M_B->operator*=( beta );
    M_A->operator+=( *M_B );
    Epetra_CrsMatrix* tmp = M_A->matrixPtr().get();
    ( *W ) = ( *tmp );
262
#ifdef HAVE_LIFEV_DEBUG
263
    W->Print(std::cout);
264
#endif
265
    return true;
266
267
}

268
void RythmosModelInterface::extractSolution( const Real &t,
269
270
271
                                        const vectorEpetra_Type& y,
                                        const vectorEpetra_Type& yp )
{
272
    M_circuitData->extractSolutionFromY( t, y, yp );
273
}
274
275
276

#endif /* HAVE_NOX_THYRA && HAVE_TRILINOS_RYTHMOS */

277
} // LifeV namespace