ZeroDimensionalSolver.cpp 13.5 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 solver
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/ZeroDimensionalSolver.hpp>
39

40
namespace LifeV {
41

42
43
#if ( defined(HAVE_NOX_THYRA) && defined(HAVE_TRILINOS_RYTHMOS) )

44
ZeroDimensionalSolver::ZeroDimensionalSolver( Int numCircuitElements,
45
46
                boost::shared_ptr< Epetra_Comm > comm,
                zeroDimensionalCircuitDataPtr_Type circuitData )
47
48
{
    M_comm.swap( comm );
49
    M_commRCP.reset();
50
    M_commRCP = Teuchos::rcp( M_comm );
51
    M_modelInterface.reset( new RythmosModelInterface( numCircuitElements,
52
53
                                    M_comm.operator ->(),
                                    circuitData ) );
54
    M_modelInterfaceRCP = Teuchos::rcp( M_modelInterface );
55

56
    M_solverInterface.reset( new RythmosSolverInterface( numCircuitElements,
57
58
                                    M_commRCP,
                                    M_modelInterfaceRCP ) );
59
60

    M_solverInterfaceRCP.reset();
61
    M_solverInterfaceRCP = Teuchos::rcp( M_solverInterface );
62

63
    M_out = Teuchos::VerboseObjectBase::getDefaultOStream();
64
65

}
66
void ZeroDimensionalSolver::setup( const ZeroDimensionalData::solverData_Type& data )
67
{
68
    std::string commandLine = "--linear-solver-params-used-file=";
69
    commandLine.append( data.linearSolverParamsFile );
70
    char * argv[1];
71
72
73
74
    argv[0] = new char[commandLine.size()];
    for ( Int i = 0; i < ( (Int) commandLine.size() ); ++i )
    {
        ( argv[0] )[i] = commandLine[i];
75
    }
76
    Int argc = 1;
77
    bool success = true; // determine if the run was successfull
78
79
80
81
    try
    { // catch exceptions
        if (data.fixTimeStep)
        {
82
            M_stepMethod = STEP_METHOD_FIXED;
83
84
85
        }
        else
        {
86
            M_stepMethod = STEP_METHOD_VARIABLE;
87
        }
88

89
        Int numElements = M_modelInterface->numCircuitElements(); // number of elements in vector
90
91
92
93
94
95
        EMethod method_val = METHOD_BE;
        Real maxError;
        Real reltol;
        Real abstol;
        Int maxOrder;
        bool useNOX;
96

97
        if (!data.method.compare("BE"))
98
        method_val = METHOD_BE;
99
        if (!data.method.compare("BDF"))
100
        method_val = METHOD_BDF;
101
        if (!data.method.compare("IRK"))
102
103
        method_val = METHOD_IRK;

104
105
106
107
108
        M_numberTimeStep = data.numberTimeStep;
        maxError = data.maxError;
        reltol = data.reltol;
        abstol = data.abstol;
        maxOrder = data.maxOrder;
109
        M_outputLevel = data.verboseLevel;
110
111
        M_outputLevel = min(max(M_outputLevel,-1),4);
        useNOX = data.useNOX;
112

113
114
115
        switch (M_outputLevel)
        {
            case -1:
116
117
            M_outputLevelTeuchos = Teuchos::VERB_DEFAULT;
            break;
118
            case 0:
119
            M_outputLevelTeuchos = Teuchos::VERB_NONE;
120
            break;
121
            case 1:
122
            M_outputLevelTeuchos = Teuchos::VERB_LOW;
123
            break;
124
            case 2:
125
            M_outputLevelTeuchos = Teuchos::VERB_MEDIUM;
126
            break;
127
            case 3:
128
            M_outputLevelTeuchos = Teuchos::VERB_HIGH;
129
            break;
130
            case 4:
131
            M_outputLevelTeuchos = Teuchos::VERB_EXTREME;
132
            break;
133
            default:
134
            break;
135
        }
136
137
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
        std::string extraLSParamsFile = data.extraLSParamsFile;

        // Parse the command-line options:

        Teuchos::CommandLineProcessor clp(false); // Don't throw exceptions
        clp.addOutputSetupOptions(true);
        Stratimikos::DefaultLinearSolverBuilder lowsfCreator;
        lowsfCreator.setupCLP(&clp);
        Teuchos::CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
        delete argv [0];
        if( parse_return != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL ) return;
        lowsfCreator.readParameters(M_out.get());
        if(extraLSParamsFile.length())
        Teuchos::updateParametersFromXmlFile( "./"+extraLSParamsFile, &*lowsfCreator.getNonconstParameterList() );
        *M_out << "\nThe parameter list after being read in:\n";
        lowsfCreator.getParameterList()->print(*M_out,2,true,false);

        // Set up the parameter list for the application:
        Teuchos::ParameterList params;
        params.set( "NumElements", numElements );
        // Create the factory for the LinearOpWithSolveBase object
        Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<Real> >
        W_factory;
        W_factory = lowsfCreator.createLinearSolveStrategy("");
        *M_out
        << "\nCreated a LinearOpWithSolveFactory described as:\n"
        << Teuchos::describe(*W_factory,Teuchos::VERB_MEDIUM);

        // create interface to problem
        Teuchos::RCP<Thyra::ModelEvaluator<Real> >
        model = Teuchos::rcp(new Thyra::EpetraModelEvaluator(M_solverInterfaceRCP,W_factory));

        Thyra::ModelEvaluatorBase::InArgs<Real> model_ic = model->getNominalValues();

        std::string method;
        if (method_val == METHOD_BE)
172
        {
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
            Teuchos::RCP<Thyra::NonlinearSolverBase<Real> > nonlinearSolver;
            if (useNOX)
            {
                Teuchos::RCP<Thyra::NOXNonlinearSolver> _nonlinearSolver = Teuchos::rcp(new Thyra::NOXNonlinearSolver);
                nonlinearSolver = _nonlinearSolver;
            }
            else
            {
                Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Real> >
                _nonlinearSolver = Teuchos::rcp(new Rythmos::TimeStepNonlinearSolver<Real>());
                Teuchos::RCP<Teuchos::ParameterList>
                nonlinearSolverPL = Teuchos::parameterList();
                nonlinearSolverPL->set("Default Tol",Real(maxError));
                _nonlinearSolver->setParameterList(nonlinearSolverPL);
                nonlinearSolver = _nonlinearSolver;
            }
            model->setDefaultVerbLevel(M_outputLevelTeuchos);
            nonlinearSolver->setVerbLevel(M_outputLevelTeuchos);
            M_stepperPtr = Teuchos::rcp(new Rythmos::BackwardEulerStepper<Real>(model,nonlinearSolver));
            M_stepperPtr->setVerbLevel(M_outputLevelTeuchos);
            model->describe(*M_out,M_outputLevelTeuchos);
            model->description();
            nonlinearSolver->describe(*M_out,M_outputLevelTeuchos);
            nonlinearSolver->description();
            M_stepperPtr->describe(*M_out,M_outputLevelTeuchos);
            M_stepperPtr->description();
            method = "Backward Euler";
200
        }
201
        else if (method_val == METHOD_BDF)
202
        {
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
            Teuchos::RCP<Thyra::NonlinearSolverBase<Real> > nonlinearSolver;
            if (useNOX)
            {
                Teuchos::RCP<Thyra::NOXNonlinearSolver> _nonlinearSolver = Teuchos::rcp(new Thyra::NOXNonlinearSolver);
                nonlinearSolver = _nonlinearSolver;
            }
            else
            {
                Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Real> >
                _nonlinearSolver = Teuchos::rcp(new Rythmos::TimeStepNonlinearSolver<Real>());
                Teuchos::RCP<Teuchos::ParameterList>
                nonlinearSolverPL = Teuchos::parameterList();
                nonlinearSolverPL->set("Default Tol",Real(maxError));
                _nonlinearSolver->setParameterList(nonlinearSolverPL);
                nonlinearSolver = _nonlinearSolver;
            }
            nonlinearSolver->setVerbLevel(M_outputLevelTeuchos);
            Teuchos::RCP<Teuchos::ParameterList> BDFparams = Teuchos::rcp(new Teuchos::ParameterList);
            Teuchos::RCP<Teuchos::ParameterList> BDFStepControlPL = Teuchos::sublist(BDFparams, "Step Control Settings");
            BDFStepControlPL->set( "maxOrder", maxOrder );
            BDFStepControlPL->set( "relErrTol", reltol );
            BDFStepControlPL->set( "absErrTol", abstol );
            M_stepperPtr = Teuchos::rcp(new Rythmos::ImplicitBDFStepper<Real>(model,nonlinearSolver,BDFparams));
            method = "Implicit BDF";
227
        }
228
        else if (method_val == METHOD_IRK)
229
        {
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
            Teuchos::RCP<Thyra::NonlinearSolverBase<Real> > nonlinearSolver;
            if (useNOX)
            {
                Teuchos::RCP<Thyra::NOXNonlinearSolver> _nonlinearSolver = Teuchos::rcp(new Thyra::NOXNonlinearSolver);
                nonlinearSolver = _nonlinearSolver;
            }
            else
            {
                Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Real> >
                _nonlinearSolver = Teuchos::rcp(new Rythmos::TimeStepNonlinearSolver<Real>());
                Teuchos::RCP<Teuchos::ParameterList>
                nonlinearSolverPL = Teuchos::parameterList();
                nonlinearSolverPL->set("Default Tol",Real(maxError));
                _nonlinearSolver->setParameterList(nonlinearSolverPL);
                nonlinearSolver = _nonlinearSolver;
            }
            nonlinearSolver->setVerbLevel(M_outputLevelTeuchos);
            Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<Real> > irk_W_factory
            = lowsfCreator.createLinearSolveStrategy("");
            Teuchos::RCP<Rythmos::RKButcherTableauBase<Real> > rkbt = Rythmos::createRKBT<Real>("Implicit 3 Stage 6th order Gauss");
            M_stepperPtr = Rythmos::implicitRKStepper<Real>(model,nonlinearSolver,irk_W_factory,rkbt);
            method = "Implicit RK";
252
253
254
        }
        else
        {
255
            TEST_FOR_EXCEPT(true);
256
        }
257

258
        M_stepperPtr->setInitialCondition(model_ic);
259

260
261
262
263
        M_method = method;
    }
    TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*M_out,success)
    return;
264
265
266

}
void
267
268
ZeroDimensionalSolver::takeStep(Real t0,Real t1)
{
269

270
271
272
273
274
275
276
    M_finalTime = t1;
    M_startTime = t0;
    Real dt = (t1-t0)/M_numberTimeStep;
    Real time = t0;
    Int numSteps = 0;
    Thyra::ModelEvaluatorBase::InArgs< Real > myIn = M_stepperPtr->getInitialCondition();
    myIn.describe(*M_out, M_outputLevelTeuchos);
277

278
    if (M_stepMethod == STEP_METHOD_FIXED)
279
    {
280
281
        // Integrate forward with fixed step sizes:
        for (Int i=1; i<=M_numberTimeStep; ++i)
282
        {
283
284
285
286
287
288
289
290
            Real dt_taken = M_stepperPtr->takeStep(dt,Rythmos::STEP_TYPE_FIXED);
            numSteps++;
            if (dt_taken != dt)
            {
                cerr << "Error, M_stepper took step of dt = " << dt_taken << " when asked to take step of dt = " << dt << std::endl;
                break;
            }
            time += dt_taken;
291
292
        }
    }
293
    else // M_stepMethod == STEP_METHOD_VARIABLE
294

295
    {
296
        while (time < M_finalTime)
297
        {
298
299
300
301
302
303
304
305
306
            Real dt_taken = M_stepperPtr->takeStep(0.0,Rythmos::STEP_TYPE_VARIABLE);
            numSteps++;
            if (dt_taken < 0)
            {
                cerr << "Error, M_stepper failed for some reason with step taken = " << dt_taken << endl;
                break;
            }
            time += dt_taken;
            *M_out << "Took stepsize of: " << dt_taken << " time = " << time << endl;
307
308
        }
    }
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
    // Get solution M_out of M_stepper:
    Rythmos::TimeRange< Real > timeRange =M_stepperPtr->getTimeRange();
    Rythmos::Array< Real > time_vec;
    Real timeToEvaluate = timeRange.upper();
    time_vec.push_back(timeToEvaluate);
    Rythmos::Array< Teuchos::RCP< const Thyra::VectorBase< Real > > > x_vec;
    Rythmos::Array< Teuchos::RCP< const Thyra::VectorBase< Real > > > xdot_vec;
    Rythmos::Array< Teuchos::ScalarTraits<Real>::magnitudeType > accuracy_vec;

    M_stepperPtr->getPoints(time_vec, &x_vec, &xdot_vec, &accuracy_vec
    );
    const Rythmos::StepStatus<Real> stepStatus = M_stepperPtr->getStepStatus();
    Teuchos::RCP<const Thyra::VectorBase<Real> > &x_computed_thyra_ptr = x_vec[0];
    Teuchos::RCP<const Thyra::VectorBase<Real> > &x_dot_computed_thyra_ptr = xdot_vec[0];

    // Convert Thyra::VectorBase to Epetra_Vector
    Teuchos::RCP<const Epetra_Vector> x_computed_ptr = Thyra::get_Epetra_Vector(*(M_solverInterfaceRCP->get_x_map()),x_computed_thyra_ptr);
    const Epetra_Vector &x_computed = *x_computed_ptr;

    Teuchos::RCP<const Epetra_Vector> x_dot_computed_ptr = Thyra::get_Epetra_Vector(*(M_solverInterfaceRCP->get_x_map()),x_dot_computed_thyra_ptr);
    const Epetra_Vector &x_dot_computed = *x_dot_computed_ptr;

    //---------------Replace the solution as initial solution for the next iteration
332
333
334
335
336
337
#ifdef HAVE_LIFEV_DEBUG
    *M_out << "X_computed" << endl;
    x_computed.Print(*M_out);
    *M_out << "X_dot_computed" << endl;
    x_dot_computed.Print(*M_out);
#endif
338

339
340
341
    M_modelInterface->initializeSolnY(x_computed);
    M_modelInterface->initializeSolnYp(x_dot_computed);
    M_modelInterface->extractSolution(time,x_computed , x_dot_computed);
342
}
343

344
345
#endif /* HAVE_NOX_THYRA && HAVE_TRILINOS_RYTHMOS */

346
} // LifeV namespace
347
348