ZeroDimensionalData.cpp 11.3 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 File containing a class for 0D model data handling.
30
31
32
33
 *
 *  @date 16-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/ZeroDimensionalData.hpp>
39

40
namespace LifeV {
41
42
43
44

// ===================================================
// Constructors
// ===================================================
45
46
47
48
49
ZeroDimensionalData::ZeroDimensionalData() :
    M_time(), M_outPutFormat( OutPutFormat( "20",
                                            "12",
                                            "   ",
                                            5000 ) ), M_circuitData( new ZeroDimensionalCircuitData )
50
51
52
53
54
55
56
57
58
59
60
61
{
}
ZeroDimensionalData::~ZeroDimensionalData()
{
    M_voltageFileStream.close();
    M_currentFileStream.close();
    M_balanceFileStream.close();
}

// ===================================================
// Methods
// ===================================================
62
63
64
void ZeroDimensionalData::setup( const GetPot& dataFile,
                                 bcPtr_Type bc,
                                 const std::string & section )
65
66
{
    if ( !M_time.get() )
67
68
        M_time.reset( new time_Type( dataFile,
                                     section + "/time_discretization" ) );
69

70
71
    std::string circuirDataFile = dataFile( ( section + "/CircuitDataFile" ).data(),
                                            "./inputFilexx.dat" );
72
    GetPot dataFileCircuit( circuirDataFile );
73
74
75
76
77
78
79
80
81
82
83
84
85
86
    std::string circuitFile = dataFileCircuit( "Files/InputFiles/CircuitFile",
                                               "./circuitFilexx.dat" );
    std::string folder = dataFileCircuit( "Files/OutputFiles/Folder",
                                          "outputxx" );//TODO create folder output

    std::string voltageFile = folder + "/" + dataFileCircuit( "Files/OutputFiles/VoltageFile",
                                                              "voltagexx.txt" );
    std::string currentFile = folder + "/" + dataFileCircuit( "Files/OutputFiles/CurrentFile",
                                                              "currentxx.txt" );
    std::string balanceFile = folder + "/" + dataFileCircuit( "Files/OutputFiles/BalanceCurrentinNode",
                                                              "balancexx.txt" );
    M_voltageFileStream.open( voltageFile.c_str() );
    M_currentFileStream.open( currentFile.c_str() );
    M_balanceFileStream.open( balanceFile.c_str() );
87
88

    //Build the Circuit
89
90
    M_circuitData->buildCircuit( circuitFile.c_str(),
                                 bc );
91
92
93
94

    //assign varible index
    assignVaribleIndex();

95
96
97
98
99
100
101
    M_solverData.method = dataFile( ( section + "/Solver/method" ).data(), "IRK" );
    M_solverData.numberTimeStep = dataFile( ( section + "/Solver/numberTimeStep" ).data(), 1 );
    M_solverData.maxError = dataFile( ( section + "/Solver/maxError" ).data(), 1.0 );
    M_solverData.reltol = dataFile( ( section + "/Solver/reltol" ).data(), 0.1 );
    M_solverData.abstol = dataFile( ( section + "/Solver/abstol" ).data(), 1.0 );
    M_solverData.maxOrder = dataFile( ( section + "/Solver/maxOrder" ).data(), 1 );
    M_solverData.verbose = dataFile( ( section + "/Solver/verbose" ).data(), false );
102
    M_solverData.verboseLevel = dataFile( ( section + "/Solver/verboseLevel" ).data(), 0 );
103
104
105
106
107
108
109
    M_solverData.useNOX = dataFile( ( section + "/Solver/useNOX" ).data(), false );
    M_solverData.fixTimeStep = dataFile( ( section + "/Solver/fixTimeStep" ).data(), true );
    M_solverData.extraLSParamsFile = dataFile( ( section + "/Solver/extraLinearSolverParamsFile" ).data(), "./Extra_AztecOO_Params.xml" );
    M_solverData.linearSolverParamsFile = dataFile( ( section + "/Solver/linearSolverParamsUsedFile" ).data(), "./lowsf.aztecoo.used.xml" );

    //Set zero initial condition
    // TODO: change to general initial condition
110
111
112
113
114
115
    initializeSolution();

    //Wrire Header OutputFiles
    writeHeaders();
}

116
void ZeroDimensionalData::initializeSolution()
117
{
118
119
    ptrVecZeroDimensionalElementPtr_Type elementList = M_circuitData->Elements() ->elementList();
    for ( iterZeroDimensionalElement_Type theElement = elementList->begin(); theElement != elementList->end(); theElement++ )
120
    {
121
122
        ( *theElement )->setCurrent( 0.0 );
        ( *theElement )->setDeltaCurrent( 0.0 );
123
124
    }

125
126
    ptrVecZeroDimensionalNodePtr_Type nodeList = M_circuitData->Nodes() ->nodeList();
    for ( iterZeroDimensionalNode_Type theNode = nodeList->begin(); theNode != nodeList->end(); theNode++ )
127
    {
128
129
        ( *theNode )->setVoltage( 0.0 );
        ( *theNode )->setDeltaVoltage( 0.0 );
130
131
132
    }
}

133
134
//! update source elements
void ZeroDimensionalData::updateBC()
135
{
136
137
    Real time = M_time->time();
    ptrVecZeroDimensionalElementCurrentSourcePtr_Type currentElementList = M_circuitData->Elements()-> currentSourceList();
138

139
    for ( iterZeroDimensionalElementCurrentSource_Type theElement = currentElementList->begin(); theElement != currentElementList->end(); theElement++ )
140
    {
141
142
        ( *theElement )->setCurrentByTime( time );
        std::cout << ( *theElement )->current() << std::endl;
143
144
    }

145
146
    ptrVecZeroDimensionalElementVoltageSourcePtr_Type voltageElementList = M_circuitData->Elements()-> voltageSourceList();
    for ( iterZeroDimensionalElementVoltageSourcePtr_Type theElement = voltageElementList->begin(); theElement != voltageElementList->end(); theElement++ )
147
    {
148
149
        ( *theElement )->setVoltageByTime( time );
        std::cout << ( *theElement )->voltage() << std::endl;
150
151
152
    }
}

153
154
155
156
157
158
void ZeroDimensionalData::saveSolution()
{
    ptrVecZeroDimensionalElementPtr_Type elementList = M_circuitData->Elements() ->elementList();
    M_outPutFormat.writeDataFormat( M_time->time(),
                                    M_currentFileStream,
                                    M_outPutFormat.space );
159
160

    //write current
161
162
163
164
165
    for ( iterZeroDimensionalElement_Type theElement = elementList->begin(); theElement != elementList->end(); theElement++ )
    {
        M_outPutFormat.writeDataFormat( ( *theElement )->current(),
                                        M_currentFileStream,
                                        M_outPutFormat.space );
166
167
168
    }

    //write voltage and current balance at each node
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
    M_outPutFormat.writeNewLine( M_currentFileStream );
    ptrVecZeroDimensionalNodePtr_Type nodeList = M_circuitData->Nodes() ->nodeList();
    M_outPutFormat.writeDataFormat( M_time->time(),
                                    M_voltageFileStream,
                                    M_outPutFormat.space );
    M_outPutFormat.writeDataFormat( M_time->time(),
                                    M_balanceFileStream,
                                    M_outPutFormat.space );

    for ( iterZeroDimensionalNode_Type theNode = nodeList->begin(); theNode != nodeList->end(); theNode++ )
    {
        M_outPutFormat.writeDataFormat( ( *theNode )->voltage(),
                                        M_voltageFileStream,
                                        M_outPutFormat.space );
        M_outPutFormat.writeDataFormat( ( *theNode )->currentBalance(),
                                        M_balanceFileStream,
                                        M_outPutFormat.space );
186
187
    }

188
189
    M_outPutFormat.writeNewLine( M_voltageFileStream );
    M_outPutFormat.writeNewLine( M_balanceFileStream );
190
191
192
}


193
void ZeroDimensionalData::showMeVariables()
194
{
195
196
    std::cout << "This MultiscaleModel0D::ShowVariables----------------------" << std::endl;
    zeroDimensionalNodeSPtr_Type Nodes = M_circuitData->Nodes();
197
198
    zeroDimensionalElementSPtr_Type Elements = M_circuitData->Elements();

199
200
201
    for ( iterZeroDimensionalNodeUnknown_Type theUnknownNode = Nodes->unknownNodeList()->begin(); theUnknownNode != Nodes->unknownNodeList()->end(); theUnknownNode++ )
    {
        ( *theUnknownNode )->showMe( 1 );
202
    }
203

204
205
206
207
    for ( iterZeroDimensionalElementPassiveInductor_Type theInductor = Elements->inductorList() ->begin(); theInductor
                    != Elements->inductorList() ->end(); theInductor++ )
    {
        ( *theInductor )->showMe( 1 );
208
209
210
    }
}

211
void ZeroDimensionalData::assignVaribleIndex()
212
{
213
214
215
216
217
    zeroDimensionalNodeSPtr_Type Nodes = M_circuitData->Nodes();
    ptrVecZeroDimensionalNodeUnknownPtr_Type unKnownList = Nodes->unknownNodeList();
    zeroDimensionalElementSPtr_Type Elements = M_circuitData->Elements();
    ptrVecZeroDimensionalElementPassiveInductorPtr_Type inductorList = Elements->inductorList();
    M_unknownCounter = 0;
218

219
220
    //set variable index for unknown nodes
    for ( iterZeroDimensionalNodeUnknown_Type theUnknownNode = unKnownList->begin(); theUnknownNode != unKnownList->end(); theUnknownNode++ )
221
    {
222
223
        ( *theUnknownNode )->assignVariableIndex( M_unknownCounter );
        M_unknownCounter++;
224
    }
225

226
227
    //set variable index for inductor current
    for ( iterZeroDimensionalElementPassiveInductor_Type theInductor = inductorList->begin(); theInductor != inductorList->end(); theInductor++ )
228
    {
229
230
        ( *theInductor )->assignVariableIndex( M_unknownCounter );
        M_unknownCounter++;
231
232
    }
}
233

234
235
236
237
void ZeroDimensionalData::writeHeaders()
{
    //write header for current file
    ptrVecZeroDimensionalElementPtr_Type elementList = M_circuitData->Elements() ->elementList();
238
    M_outPutFormat.writeDataFormat( "% time", M_currentFileStream, M_outPutFormat.space );
239
240
241
242
243
244
245
246
    for ( iterZeroDimensionalElement_Type theElement = elementList->begin(); theElement != elementList->end(); theElement++ )
    {
        M_outPutFormat.writeDataFormat( ( *theElement )->id(), M_currentFileStream, M_outPutFormat.space );
    }
    M_outPutFormat.writeNewLine( M_currentFileStream );

    //write header for voltage and current balance at each node (voltage file and balance file)
    ptrVecZeroDimensionalNodePtr_Type nodeList = M_circuitData->Nodes() ->nodeList();
247
248
    M_outPutFormat.writeDataFormat( "% time", M_balanceFileStream, M_outPutFormat.space );
    M_outPutFormat.writeDataFormat( "% time", M_voltageFileStream, M_outPutFormat.space );
249
250
251
252
253
254
255
256
257
    for ( iterZeroDimensionalNode_Type theNode = nodeList->begin(); theNode != nodeList->end(); theNode++ )
    {
        M_outPutFormat.writeDataFormat( ( *theNode )->id(), M_voltageFileStream, M_outPutFormat.space );
        M_outPutFormat.writeDataFormat( ( *theNode )->id(), M_balanceFileStream, M_outPutFormat.space );
    }
    M_outPutFormat.writeNewLine( M_voltageFileStream );
    M_outPutFormat.writeNewLine( M_balanceFileStream );
}

258
} // LifeV namespace