LifeV
ZeroDimensionalSolver.cpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4 
5  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
6  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7 
8  This file is part of LifeV.
9 
10  LifeV is free software; you can redistribute it and/or modify
11  it under the terms of the GNU Lesser General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  LifeV is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 
27 /*!
28  * @file
29  * @brief Rythmos solver
30  *
31  * @date 16-11-2011
32  * @author Mahmoud Jafargholi
33  *
34  * @contributors Cristiano Malossi <cristiano.malossi@epfl.ch>
35  * @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 
45 ZeroDimensionalSolver::ZeroDimensionalSolver ( Int numCircuitElements,
46  std::shared_ptr< Epetra_Comm > comm,
47  zeroDimensionalCircuitDataPtr_Type circuitData )
48 {
49  M_comm.swap ( comm );
50  M_commRCP.reset();
51  M_commRCP = Teuchos::rcp ( M_comm );
52  M_modelInterface.reset ( new RythmosModelInterface ( numCircuitElements,
53  M_comm.operator ->(),
54  circuitData ) );
55  M_modelInterfaceRCP = Teuchos::rcp ( M_modelInterface );
56 
57  M_solverInterface.reset ( new RythmosSolverInterface ( numCircuitElements,
58  M_commRCP,
59  M_modelInterfaceRCP ) );
60 
61  M_solverInterfaceRCP.reset();
62  M_solverInterfaceRCP = Teuchos::rcp ( M_solverInterface );
63 
64  M_out = Teuchos::VerboseObjectBase::getDefaultOStream();
65 
66 }
67 void ZeroDimensionalSolver::setup ( const ZeroDimensionalData::solverData_Type& data )
68 {
69  std::string commandLine = "--linear-solver-params-used-file=";
70  commandLine.append ( data.linearSolverParamsFile );
71  char* argv[1];
72  argv[0] = new char[commandLine.size()];
73  for ( Int i = 0; i < ( (Int) commandLine.size() ); ++i )
74  {
75  ( argv[0] ) [i] = commandLine[i];
76  }
77  Int argc = 1;
78  bool success = true; // determine if the run was successfull
79  try
80  {
81  // catch exceptions
82  if (data.fixTimeStep)
83  {
84  M_stepMethod = STEP_METHOD_FIXED;
85  }
86  else
87  {
88  M_stepMethod = STEP_METHOD_VARIABLE;
89  }
90 
91  Int numElements = M_modelInterface->numCircuitElements(); // number of elements in vector
92  EMethod method_val = METHOD_BE;
93  Real maxError;
94  Real reltol;
95  Real abstol;
96  Int maxOrder;
97  bool useNOX;
98 
99  if (!data.method.compare ("BE") )
100  {
101  method_val = METHOD_BE;
102  }
103  if (!data.method.compare ("BDF") )
104  {
105  method_val = METHOD_BDF;
106  }
107  if (!data.method.compare ("IRK") )
108  {
109  method_val = METHOD_IRK;
110  }
111 
112  M_numberTimeStep = data.numberTimeStep;
113  maxError = data.maxError;
114  reltol = data.reltol;
115  abstol = data.abstol;
116  maxOrder = data.maxOrder;
117  M_outputLevel = data.verboseLevel;
118  M_outputLevel = min (max (M_outputLevel, -1), 4);
119  useNOX = data.useNOX;
120 
121  switch (M_outputLevel)
122  {
123  case -1:
124  M_outputLevelTeuchos = Teuchos::VERB_DEFAULT;
125  break;
126  case 0:
127  M_outputLevelTeuchos = Teuchos::VERB_NONE;
128  break;
129  case 1:
130  M_outputLevelTeuchos = Teuchos::VERB_LOW;
131  break;
132  case 2:
133  M_outputLevelTeuchos = Teuchos::VERB_MEDIUM;
134  break;
135  case 3:
136  M_outputLevelTeuchos = Teuchos::VERB_HIGH;
137  break;
138  case 4:
139  M_outputLevelTeuchos = Teuchos::VERB_EXTREME;
140  break;
141  default:
142  break;
143  }
144  std::string extraLSParamsFile = data.extraLSParamsFile;
145 
146  // Parse the command-line options:
147 
148  Teuchos::CommandLineProcessor clp (false); // Don't throw exceptions
149  clp.addOutputSetupOptions (true);
150  Stratimikos::DefaultLinearSolverBuilder lowsfCreator;
151  lowsfCreator.setupCLP (&clp);
152  Teuchos::CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse (argc, argv);
153  delete argv [0];
154  if ( parse_return != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL )
155  {
156  return;
157  }
158  lowsfCreator.readParameters (M_out.get() );
159  if (extraLSParamsFile.length() )
160  {
161  Teuchos::updateParametersFromXmlFile ( "./" + extraLSParamsFile, &*lowsfCreator.getNonconstParameterList() );
162  }
163  *M_out << "\nThe parameter list after being read in:\n";
164  lowsfCreator.getParameterList()->print (*M_out, 2, true, false);
165 
166  // Set up the parameter list for the application:
167  Teuchos::ParameterList params;
168  params.set ( "NumElements", numElements );
169  // Create the factory for the LinearOpWithSolveBase object
170  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<Real> >
171  W_factory;
172  W_factory = lowsfCreator.createLinearSolveStrategy ("");
173  *M_out
174  << "\nCreated a LinearOpWithSolveFactory described as:\n"
175  << Teuchos::describe (*W_factory, Teuchos::VERB_MEDIUM);
176 
177  // create interface to problem
178  Teuchos::RCP<Thyra::ModelEvaluator<Real> >
179  model = Teuchos::rcp (new Thyra::EpetraModelEvaluator (M_solverInterfaceRCP, W_factory) );
180 
181  Thyra::ModelEvaluatorBase::InArgs<Real> model_ic = model->getNominalValues();
182 
183  std::string method;
184  if (method_val == METHOD_BE)
185  {
186  Teuchos::RCP<Thyra::NonlinearSolverBase<Real> > nonlinearSolver;
187  if (useNOX)
188  {
189  Teuchos::RCP<Thyra::NOXNonlinearSolver> _nonlinearSolver = Teuchos::rcp (new Thyra::NOXNonlinearSolver);
190  nonlinearSolver = _nonlinearSolver;
191  }
192  else
193  {
194  Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Real> >
195  _nonlinearSolver = Teuchos::rcp (new Rythmos::TimeStepNonlinearSolver<Real>() );
196  Teuchos::RCP<Teuchos::ParameterList>
197  nonlinearSolverPL = Teuchos::parameterList();
198  nonlinearSolverPL->set ("Default Tol", Real (maxError) );
199  _nonlinearSolver->setParameterList (nonlinearSolverPL);
200  nonlinearSolver = _nonlinearSolver;
201  }
202  model->setDefaultVerbLevel (M_outputLevelTeuchos);
203  nonlinearSolver->setVerbLevel (M_outputLevelTeuchos);
204  M_stepperPtr = Teuchos::rcp (new Rythmos::BackwardEulerStepper<Real> (model, nonlinearSolver) );
205  M_stepperPtr->setVerbLevel (M_outputLevelTeuchos);
206  model->describe (*M_out, M_outputLevelTeuchos);
207  model->description();
208  nonlinearSolver->describe (*M_out, M_outputLevelTeuchos);
209  nonlinearSolver->description();
210  M_stepperPtr->describe (*M_out, M_outputLevelTeuchos);
211  M_stepperPtr->description();
212  method = "Backward Euler";
213  }
214  else if (method_val == METHOD_BDF)
215  {
216  Teuchos::RCP<Thyra::NonlinearSolverBase<Real> > nonlinearSolver;
217  if (useNOX)
218  {
219  Teuchos::RCP<Thyra::NOXNonlinearSolver> _nonlinearSolver = Teuchos::rcp (new Thyra::NOXNonlinearSolver);
220  nonlinearSolver = _nonlinearSolver;
221  }
222  else
223  {
224  Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Real> >
225  _nonlinearSolver = Teuchos::rcp (new Rythmos::TimeStepNonlinearSolver<Real>() );
226  Teuchos::RCP<Teuchos::ParameterList>
227  nonlinearSolverPL = Teuchos::parameterList();
228  nonlinearSolverPL->set ("Default Tol", Real (maxError) );
229  _nonlinearSolver->setParameterList (nonlinearSolverPL);
230  nonlinearSolver = _nonlinearSolver;
231  }
232  nonlinearSolver->setVerbLevel (M_outputLevelTeuchos);
233  Teuchos::RCP<Teuchos::ParameterList> BDFparams = Teuchos::rcp (new Teuchos::ParameterList);
234  Teuchos::RCP<Teuchos::ParameterList> BDFStepControlPL = Teuchos::sublist (BDFparams, "Step Control Settings");
235  BDFStepControlPL->set ( "maxOrder", maxOrder );
236  BDFStepControlPL->set ( "relErrTol", reltol );
237  BDFStepControlPL->set ( "absErrTol", abstol );
238  M_stepperPtr = Teuchos::rcp (new Rythmos::ImplicitBDFStepper<Real> (model, nonlinearSolver, BDFparams) );
239  method = "Implicit BDF";
240  }
241  else if (method_val == METHOD_IRK)
242  {
243  Teuchos::RCP<Thyra::NonlinearSolverBase<Real> > nonlinearSolver;
244  if (useNOX)
245  {
246  Teuchos::RCP<Thyra::NOXNonlinearSolver> _nonlinearSolver = Teuchos::rcp (new Thyra::NOXNonlinearSolver);
247  nonlinearSolver = _nonlinearSolver;
248  }
249  else
250  {
251  Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Real> >
252  _nonlinearSolver = Teuchos::rcp (new Rythmos::TimeStepNonlinearSolver<Real>() );
253  Teuchos::RCP<Teuchos::ParameterList>
254  nonlinearSolverPL = Teuchos::parameterList();
255  nonlinearSolverPL->set ("Default Tol", Real (maxError) );
256  _nonlinearSolver->setParameterList (nonlinearSolverPL);
257  nonlinearSolver = _nonlinearSolver;
258  }
259  nonlinearSolver->setVerbLevel (M_outputLevelTeuchos);
260  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<Real> > irk_W_factory
261  = lowsfCreator.createLinearSolveStrategy ("");
262  Teuchos::RCP<Rythmos::RKButcherTableauBase<Real> > rkbt = Rythmos::createRKBT<Real> ("Implicit 3 Stage 6th order Gauss");
263  M_stepperPtr = Rythmos::implicitRKStepper<Real> (model, nonlinearSolver, irk_W_factory, rkbt);
264  method = "Implicit RK";
265  }
266  else
267  {
268  TEST_FOR_EXCEPT (true);
269  }
270 
271  M_stepperPtr->setInitialCondition (model_ic);
272 
273  M_method = method;
274  }
275  TEUCHOS_STANDARD_CATCH_STATEMENTS (true, *M_out, success)
276  return;
277 
278 }
279 void
280 ZeroDimensionalSolver::takeStep (Real t0, Real t1)
281 {
282 
283  M_finalTime = t1;
284  M_startTime = t0;
285  Real dt = (t1 - t0) / M_numberTimeStep;
286  Real time = t0;
287  Int numSteps = 0;
288  Thyra::ModelEvaluatorBase::InArgs< Real > myIn = M_stepperPtr->getInitialCondition();
289  myIn.describe (*M_out, M_outputLevelTeuchos);
290 
291  if (M_stepMethod == STEP_METHOD_FIXED)
292  {
293  // Integrate forward with fixed step sizes:
294  for (Int i = 1; i <= M_numberTimeStep; ++i)
295  {
296  Real dt_taken = M_stepperPtr->takeStep (dt, Rythmos::STEP_TYPE_FIXED);
297  numSteps++;
298  if (dt_taken != dt)
299  {
300  cerr << "Error, M_stepper took step of dt = " << dt_taken << " when asked to take step of dt = " << dt << std::endl;
301  break;
302  }
303  time += dt_taken;
304  }
305  }
306  else // M_stepMethod == STEP_METHOD_VARIABLE
307 
308  {
309  while (time < M_finalTime)
310  {
311  Real dt_taken = M_stepperPtr->takeStep (0.0, Rythmos::STEP_TYPE_VARIABLE);
312  numSteps++;
313  if (dt_taken < 0)
314  {
315  cerr << "Error, M_stepper failed for some reason with step taken = " << dt_taken << endl;
316  break;
317  }
318  time += dt_taken;
319  *M_out << "Took stepsize of: " << dt_taken << " time = " << time << endl;
320  }
321  }
322  // Get solution M_out of M_stepper:
323  Rythmos::TimeRange< Real > timeRange = M_stepperPtr->getTimeRange();
324  Rythmos::Array< Real > time_vec;
325  Real timeToEvaluate = timeRange.upper();
326  time_vec.push_back (timeToEvaluate);
327  Rythmos::Array< Teuchos::RCP< const Thyra::VectorBase< Real > > > x_vec;
328  Rythmos::Array< Teuchos::RCP< const Thyra::VectorBase< Real > > > xdot_vec;
329  Rythmos::Array< Teuchos::ScalarTraits<Real>::magnitudeType > accuracy_vec;
330 
331  M_stepperPtr->getPoints (time_vec, &x_vec, &xdot_vec, &accuracy_vec
332  );
333  const Rythmos::StepStatus<Real> stepStatus = M_stepperPtr->getStepStatus();
334  Teuchos::RCP<const Thyra::VectorBase<Real> >& x_computed_thyra_ptr = x_vec[0];
335  Teuchos::RCP<const Thyra::VectorBase<Real> >& x_dot_computed_thyra_ptr = xdot_vec[0];
336 
337  // Convert Thyra::VectorBase to Epetra_Vector
338  Teuchos::RCP<const Epetra_Vector> x_computed_ptr = Thyra::get_Epetra_Vector (* (M_solverInterfaceRCP->get_x_map() ), x_computed_thyra_ptr);
339  const Epetra_Vector& x_computed = *x_computed_ptr;
340 
341  Teuchos::RCP<const Epetra_Vector> x_dot_computed_ptr = Thyra::get_Epetra_Vector (* (M_solverInterfaceRCP->get_x_map() ), x_dot_computed_thyra_ptr);
342  const Epetra_Vector& x_dot_computed = *x_dot_computed_ptr;
343 
344  //---------------Replace the solution as initial solution for the next iteration
345 #ifdef HAVE_LIFEV_DEBUG
346  *M_out << "X_computed" << endl;
347  x_computed.Print (*M_out);
348  *M_out << "X_dot_computed" << endl;
349  x_dot_computed.Print (*M_out);
350 #endif
351 
352  M_modelInterface->initializeSolnY (x_computed);
353  M_modelInterface->initializeSolnYp (x_dot_computed);
354  M_modelInterface->extractSolution (time, x_computed , x_dot_computed);
355 }
356 
357 #endif /* HAVE_NOX_THYRA && HAVE_TRILINOS_RYTHMOS */
358 
359 } // LifeV namespace
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405