Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_StepperOperatorSplit_impl.hpp
Go to the documentation of this file.
1// @HEADER
2// ****************************************************************************
3// Tempus: Copyright (2017) Sandia Corporation
4//
5// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6// ****************************************************************************
7// @HEADER
8
9#ifndef Tempus_StepperOperatorSplit_impl_hpp
10#define Tempus_StepperOperatorSplit_impl_hpp
11
12#include "Tempus_StepperFactory.hpp"
14
15
16namespace Tempus {
17
18
19template<class Scalar>
21{
22 this->setStepperName( "Operator Split");
23 this->setStepperType( "Operator Split");
24 this->setUseFSAL( false);
25 this->setICConsistency( "None");
26 this->setICConsistencyCheck( false);
27
28 this->setOrder (1);
29 this->setOrderMin(1);
30 this->setOrderMax(1);
31 this->setAppAction(Teuchos::null);
32
33 OpSpSolnHistory_ = rcp(new SolutionHistory<Scalar>());
34 OpSpSolnHistory_->setStorageLimit(2);
35 OpSpSolnHistory_->setStorageType(Tempus::STORAGE_TYPE_STATIC);
36}
37
38
39template<class Scalar>
41 std::vector<Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > > appModels,
42 std::vector<Teuchos::RCP<Stepper<Scalar> > > subStepperList,
43 bool useFSAL,
44 std::string ICConsistency,
45 bool ICConsistencyCheck,
46 int order,
47 int orderMin,
48 int orderMax,
49 const Teuchos::RCP<StepperOperatorSplitAppAction<Scalar> >& stepperOSAppAction)
50{
51 this->setStepperName( "Operator Split");
52 this->setStepperType( "Operator Split");
53 this->setUseFSAL( useFSAL);
54 this->setICConsistency( ICConsistency);
55 this->setICConsistencyCheck( ICConsistencyCheck);
56
57 this->setSubStepperList(subStepperList);
58 this->setOrder (order);
59 this->setOrderMin(orderMin);
60 this->setOrderMax(orderMax);
61
62 this->setAppAction(stepperOSAppAction);
63 OpSpSolnHistory_ = rcp(new SolutionHistory<Scalar>());
64 OpSpSolnHistory_->setStorageLimit(2);
65 OpSpSolnHistory_->setStorageType(Tempus::STORAGE_TYPE_STATIC);
66
67 if ( !(appModels.empty()) ) {
68 this->setModels(appModels);
69 this->initialize();
70 }
71}
72
73
74template<class Scalar>
76 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
77{
78 if (appModel != Teuchos::null) {
79 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
80 Teuchos::OSTab ostab(out,1,"StepperOperatorSplit::setModel()");
81 *out << "Warning -- No ModelEvaluator to set for StepperOperatorSplit, "
82 << "because it is a Stepper of Steppers.\n" << std::endl;
83 }
84}
85
86template<class Scalar>
87Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
89{
90 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > model;
91 typename std::vector<Teuchos::RCP<Stepper<Scalar> > >::const_iterator
92 subStepperIter = subStepperList_.begin();
93 for (; subStepperIter < subStepperList_.end(); subStepperIter++) {
94 model = (*subStepperIter)->getModel();
95 if (model != Teuchos::null) break;
96 }
97 if ( model == Teuchos::null ) {
98 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
99 Teuchos::OSTab ostab(out,1,"StepperOperatorSplit::getModel()");
100 *out << "Warning -- StepperOperatorSplit::getModel() "
101 << "Could not find a valid model! Returning null!" << std::endl;
102 }
103
104 return model;
105}
106
107template<class Scalar>
109 Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> > /* solver */)
110{
111 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
112 Teuchos::OSTab ostab(out,1,"StepperOperatorSplit::setSolver()");
113 *out << "Warning -- No solver to set for StepperOperatorSplit "
114 << "because it is a Stepper of Steppers.\n" << std::endl;
115
116 this->isInitialized_ = false;
117}
118
119template<class Scalar>
121 Teuchos::RCP<StepperOperatorSplitAppAction<Scalar> > appAction)
122{
123 if (appAction == Teuchos::null) {
124 // Create default appAction, otherwise keep current.
125 if (stepperOSAppAction_ == Teuchos::null) {
126 stepperOSAppAction_ =
128 }
129 } else {
130 stepperOSAppAction_ =
131 Teuchos::rcp_dynamic_cast<StepperOperatorSplitAppAction<Scalar> > (appAction, true);
132 }
133 this->isInitialized_ = false;
134}
135
136
137template<class Scalar>
139 Teuchos::RCP<Stepper<Scalar> > stepper, bool useFSAL)
140{
141 stepper->setUseFSAL(useFSAL);
142 subStepperList_.push_back(stepper);
143}
144
145
146template<class Scalar>
148 std::vector<Teuchos::RCP<Stepper<Scalar> > > subStepperList)
149{
150 using Teuchos::RCP;
151 using Teuchos::ParameterList;
152
153 subStepperList_ = subStepperList;
154
155 typename std::vector<Teuchos::RCP<Stepper<Scalar> > >::iterator
156 subStepperIter = subStepperList_.begin();
157
158 for (; subStepperIter<subStepperList_.end(); subStepperIter++) {
159 auto subStepper = *(subStepperIter);
160 bool useFSAL = subStepper->getUseFSAL();
161 if (useFSAL) {
162 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
163 Teuchos::OSTab ostab(out,1,"StepperOperatorSplit::setSubStepperList()");
164 *out << "Warning -- subStepper = '"
165 << subStepper->getStepperType() << "' has \n"
166 << " subStepper->getUseFSAL() = " << useFSAL << ".\n"
167 << " subSteppers usually can not use the FSAL priniciple with\n"
168 << " operator splitting. Proceeding with it set to true.\n"
169 << std::endl;
170 }
171 }
172
173 this->isInitialized_ = false;
174}
175
176template<class Scalar>
178 std::vector<Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > > appModels)
179{
180 using Teuchos::RCP;
181 using Teuchos::ParameterList;
182
183 TEUCHOS_TEST_FOR_EXCEPTION(subStepperList_.size() != appModels.size(),
184 std::logic_error, "Error - Number of models and Steppers do not match!\n"
185 << " There are " << appModels.size() << " models.\n"
186 << " There are " << subStepperList_.size() << " steppers.\n");
187
188 typename std::vector<RCP<const Thyra::ModelEvaluator<Scalar> > >::iterator
189 appModelIter = appModels.begin();
190
191 typename std::vector<Teuchos::RCP<Stepper<Scalar> > >::iterator
192 subStepperIter = subStepperList_.begin();
193
194 for (; appModelIter<appModels.end() || subStepperIter<subStepperList_.end();
195 appModelIter++, subStepperIter++)
196 {
197 auto appModel = *(appModelIter);
198 auto subStepper = *(subStepperIter);
199 subStepper->setModel(appModel);
200 }
201
202 this->isInitialized_ = false;
203}
204
205template<class Scalar>
207{
208 if (tempState_ == Teuchos::null) {
209 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >model = this->getModel();
210 TEUCHOS_TEST_FOR_EXCEPTION( model == Teuchos::null, std::logic_error,
211 "Error - StepperOperatorSplit::initialize() Could not find "
212 "a valid model!\n");
213 tempState_ = createSolutionStateME(model, this->getDefaultStepperState());
214 }
215
216 if (!isOneStepMethod() ) {
217 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
218 Teuchos::OSTab ostab(out,1,"StepperOperatorSplit::initialize()");
219 typename std::vector<Teuchos::RCP<Stepper<Scalar> > >::const_iterator
220 subStepperIter = subStepperList_.begin();
221 for (; subStepperIter < subStepperList_.end(); subStepperIter++) {
222 *out << "SubStepper, " << (*subStepperIter)->getStepperType()
223 << ", isOneStepMethod = " << (*subStepperIter)->isOneStepMethod()
224 << std::endl;
225 }
226 TEUCHOS_TEST_FOR_EXCEPTION(!isOneStepMethod(), std::logic_error,
227 "Error - OperatorSplit only works for one-step methods!\n");
228 }
229
230 // Ensure that subSteppers are initialized.
231 typename std::vector<Teuchos::RCP<Stepper<Scalar> > >::const_iterator
232 subStepperIter = subStepperList_.begin();
233 for (; subStepperIter < subStepperList_.end(); subStepperIter++)
234 (*subStepperIter)->initialize();
235
237}
238
239template<class Scalar>
241 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
242{
243 typename std::vector<Teuchos::RCP<Stepper<Scalar> > >::iterator
244 subStepperIter = subStepperList_.begin();
245 for (; subStepperIter < subStepperList_.end(); subStepperIter++)
246 (*subStepperIter)->setInitialConditions(solutionHistory);
247
248 Teuchos::RCP<SolutionState<Scalar> > initialState =
249 solutionHistory->getCurrentState();
250
251 // Check if we need Stepper storage for xDot
252 this->setStepperXDot(initialState->getXDot());
253 if (initialState->getXDot() == Teuchos::null)
254 this->setStepperXDot(initialState->getX()->clone_v());
255
256}
257
258template<class Scalar>
260 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
261{
262 this->checkInitialized();
263
264 using Teuchos::RCP;
265
266 TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperOperatorSplit::takeStep()");
267 {
268 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
269 std::logic_error,
270 "Error - StepperOperatorSplit<Scalar>::takeStep(...)\n"
271 "Need at least two SolutionStates for OperatorSplit.\n"
272 " Number of States = " << solutionHistory->getNumStates() << "\n"
273 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
274 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
275 RCP<StepperOperatorSplit<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
276 stepperOSAppAction_->execute(solutionHistory, thisStepper,
278
279 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
280
281 // Create OperatorSplit SolutionHistory to pass to subSteppers.
282 tempState_->copy(solutionHistory->getCurrentState());
283 OpSpSolnHistory_->clear();
284 OpSpSolnHistory_->addState(tempState_);
285 OpSpSolnHistory_->addWorkingState(workingState, false);
286
287 RCP<SolutionState<Scalar> > currentSubState =
288 OpSpSolnHistory_->getCurrentState();
289 RCP<SolutionState<Scalar> > workingSubState =
290 OpSpSolnHistory_->getWorkingState();
291
292 bool pass = true;
293 typename std::vector<Teuchos::RCP<Stepper<Scalar> > >::iterator
294 subStepperIter = subStepperList_.begin();
295 for (; subStepperIter < subStepperList_.end() && pass; subStepperIter++) {
296
297 stepperOSAppAction_->execute(solutionHistory, thisStepper,
299
300 (*subStepperIter)->takeStep(OpSpSolnHistory_);
301
302 stepperOSAppAction_->execute(solutionHistory, thisStepper,
304
305 if (workingSubState->getSolutionStatus() == Status::FAILED) {
306 pass = false;
307 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
308 Teuchos::OSTab ostab(out,1,"StepperOperatorSplit::takeStep()");
309 *out << "SubStepper, " << (*subStepperIter)->getStepperType()
310 << ", failed!" << std::endl;
311 break;
312 }
313
314 // "promote" workingSubState
315 currentSubState = OpSpSolnHistory_->getCurrentState();
316 currentSubState->copySolutionData(workingSubState);
317 }
318
319 if (pass == true) workingState->setSolutionStatus(Status::PASSED);
320 else workingState->setSolutionStatus(Status::FAILED);
321 workingState->setOrder(this->getOrder());
322 workingState->computeNorms(solutionHistory->getCurrentState());
323 OpSpSolnHistory_->clear();
324
325 stepperOSAppAction_->execute(solutionHistory, thisStepper,
327 }
328 return;
329}
330
331
338template<class Scalar>
339Teuchos::RCP<Tempus::StepperState<Scalar> > StepperOperatorSplit<Scalar>::
341{
342 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
343 rcp(new StepperState<Scalar>(this->getStepperType()));
344 return stepperState;
345}
346
347
348template<class Scalar>
350 Teuchos::FancyOStream &out,
351 const Teuchos::EVerbosityLevel verbLevel) const
352{
353 out.setOutputToRootOnly(0);
354 out << std::endl;
355 Stepper<Scalar>::describe(out, verbLevel);
356
357 out << "--- StepperOperatorSplit ---\n";
358 out << " subStepperList_.size() = " << subStepperList_.size() << std::endl;
359 typename std::vector<Teuchos::RCP<Stepper<Scalar> > >::const_iterator
360 subStepperIter = subStepperList_.begin();
361 for (; subStepperIter < subStepperList_.end(); subStepperIter++) {
362 out << " SubStepper = " << (*subStepperIter)->getStepperType()<<std::endl;
363 out << " = " << (*subStepperIter)->isInitialized() << std::endl;
364 out << " = " << (*subStepperIter) << std::endl;
365 }
366 out << " OpSpSolnHistory_ = " << OpSpSolnHistory_ << std::endl;
367 out << " tempState_ = " << tempState_ << std::endl;
368 out << " stepperOSAppAction_ = " << stepperOSAppAction_ << std::endl;
369 out << " order_ = " << order_ << std::endl;
370 out << " orderMin_ = " << orderMin_ << std::endl;
371 out << " orderMax_ = " << orderMax_ << std::endl;
372 out << "----------------------------" << std::endl;
373}
374
375
376template<class Scalar>
377bool StepperOperatorSplit<Scalar>::isValidSetup(Teuchos::FancyOStream & out) const
378{
379 out.setOutputToRootOnly(0);
380 bool isValidSetup = true;
381
382 if ( !Stepper<Scalar>::isValidSetup(out) ) isValidSetup = false;
383
384 if (subStepperList_.size() == 0) {
385 isValidSetup = false;
386 out << "The substepper list is empty!\n";
387 }
388
389 typename std::vector<Teuchos::RCP<Stepper<Scalar> > >::const_iterator
390 subStepperIter = subStepperList_.begin();
391
392 for (; subStepperIter < subStepperList_.end(); subStepperIter++) {
393 auto subStepper = *(subStepperIter);
394 if ( !subStepper->isInitialized() ) {
395 isValidSetup = false;
396 out << "The subStepper, " << subStepper->description()
397 << ", is not initialized!\n";
398 }
399 }
400 if (stepperOSAppAction_ == Teuchos::null) {
401 isValidSetup = false;
402 out << "The Operator-Split AppAction is not set!\n";
403 }
404
405 return isValidSetup;
406}
407
408
409template<class Scalar>
410Teuchos::RCP<const Teuchos::ParameterList>
412{
413 auto pl = this->getValidParametersBasic();
414 pl->template set<int>("Minimum Order", orderMin_,
415 "Minimum Operator-split order. (default = 1)\n");
416 pl->template set<int>("Order", order_,
417 "Operator-split order. (default = 1)\n");
418 pl->template set<int>("Maximum Order", orderMax_,
419 "Maximum Operator-split order. (default = 1)\n");
420
421 std::ostringstream list;
422 size_t size = subStepperList_.size();
423 for(std::size_t i = 0; i < size-1; ++i) {
424 list << "'" << subStepperList_[i]->getStepperName() << "', ";
425 }
426 list << "'" << subStepperList_[size-1]->getStepperName() << "'";
427 pl->template set<std::string>("Stepper List", list.str(),
428 "Comma deliminated list of single quoted Steppers, e.g., \"'Operator 1', 'Operator 2'\".");
429
430 for(std::size_t i = 0; i < size; ++i) {
431 pl->set(subStepperList_[i]->getStepperName(), *(subStepperList_[i]->getValidParameters()));
432 }
433
434 return pl;
435}
436
437
438template<class Scalar>
440 std::vector<Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > > appModels,
441 Teuchos::RCP<Teuchos::ParameterList> stepperPL)
442{
443 if (stepperPL != Teuchos::null) {
444 using Teuchos::RCP;
445 using Teuchos::ParameterList;
446
447 // Parse Stepper List String
448 std::vector<std::string> stepperListStr;
449 stepperListStr.clear();
450 std::string str = stepperPL->get<std::string>("Stepper List");
451 std::string delimiters(",");
452 // Skip delimiters at the beginning
453 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
454 // Find the first delimiter
455 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
456 while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
457 std::string token = str.substr(lastPos,pos-lastPos);
458 // Strip single quotes
459 std::string::size_type beg = token.find_first_of("'") + 1;
460 std::string::size_type end = token.find_last_of ("'");
461 stepperListStr.push_back(token.substr(beg,end-beg));
462
463 lastPos = str.find_first_not_of(delimiters, pos); // Skip delimiters
464 pos = str.find_first_of(delimiters, lastPos); // Find next delimiter
465 }
466
467 TEUCHOS_TEST_FOR_EXCEPTION(stepperListStr.size() != appModels.size(),
468 std::logic_error, "Error - Number of models and Steppers do not match!\n"
469 << " There are " << appModels.size() << " models.\n"
470 << " There are " << stepperListStr.size() << " steppers.\n"
471 << " " << str << "\n");
472
473 typename
474 std::vector<RCP<const Thyra::ModelEvaluator<Scalar> > >::iterator
475 aMI = appModels.begin();
476 typename std::vector<std::string>::iterator sLSI = stepperListStr.begin();
477
478 for (; aMI<appModels.end() || sLSI<stepperListStr.end(); aMI++, sLSI++) {
479 RCP<ParameterList> subStepperPL = Teuchos::sublist(stepperPL,*sLSI,true);
480 auto name = subStepperPL->name();
481 lastPos = name.rfind("->");
482 std::string newName = name.substr(lastPos+2,name.length());
483 subStepperPL->setName(newName);
484 bool useFSAL = subStepperPL->template get<bool>("Use FSAL",false);
485 auto sf = Teuchos::rcp(new StepperFactory<Scalar>());
486 auto subStepper = sf->createStepper(subStepperPL, *aMI);
487 if (useFSAL) {
488 Teuchos::RCP<Teuchos::FancyOStream> out =
489 Teuchos::VerboseObjectBase::getDefaultOStream();
490 Teuchos::OSTab ostab(out,1,"StepperFactory::createSubSteppers()");
491 *out << "Warning -- subStepper = '"
492 << subStepper->getStepperType() << "' has \n"
493 << " subStepper->getUseFSAL() = " << useFSAL << ".\n"
494 << " subSteppers usually can not use the FSAL priniciple with\n"
495 << " operator splitting. Proceeding with it set to true.\n"
496 << std::endl;
497 }
498 this->addStepper(subStepper, useFSAL);
499 }
500 }
501}
502
503
504// Nonmember constructor - ModelEvaluator and ParameterList
505// ------------------------------------------------------------------------
506template<class Scalar>
507Teuchos::RCP<StepperOperatorSplit<Scalar> >
509 std::vector<Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > > appModels,
510 Teuchos::RCP<Teuchos::ParameterList> pl)
511{
512 auto stepper = Teuchos::rcp(new StepperOperatorSplit<Scalar>());
513
514 if (pl != Teuchos::null) {
515 stepper->setStepperValues(pl);
516 stepper->setOrderMin(pl->get<int>("Minimum Order", 1));
517 stepper->setOrder (pl->get<int>("Order", 1));
518 stepper->setOrderMax(pl->get<int>("Maximum Order", 1));
519 }
520
521 if ( !(appModels.empty()) ) {
522 stepper->createSubSteppers(appModels, pl);
523 stepper->initialize();
524 }
525
526 return stepper;
527}
528
529
530} // namespace Tempus
531#endif // Tempus_StepperOperatorSplit_impl_hpp
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
StepperOperatorSplitAppAction class for StepperOperatorSplit.
OperatorSplit stepper loops through the Stepper list.
virtual void setModels(std::vector< Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > > appModels)
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
virtual void setAppAction(Teuchos::RCP< StepperOperatorSplitAppAction< Scalar > > appAction)
virtual void initialize()
Initialize during construction and after changing input parameters.
virtual void setSolver(Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > solver)
Set solver.
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
virtual void setSubStepperList(std::vector< Teuchos::RCP< Stepper< Scalar > > > subStepperList)
void createSubSteppers(std::vector< Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > > appModels, Teuchos::RCP< Teuchos::ParameterList > pl)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions and make them consistent.
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Get a default (initial) StepperState.
virtual void addStepper(Teuchos::RCP< Stepper< Scalar > > stepper, bool useFSAL=false)
Add Stepper to subStepper list. In most cases, subSteppers cannot use xDotOld (thus the default),...
StepperState is a simple class to hold state information about the stepper.
Thyra Base interface for time steppers.
virtual void initialize()
Initialize after construction and changing input parameters.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
@ STORAGE_TYPE_STATIC
Keep a fix number of states.
Teuchos::RCP< SolutionState< Scalar > > createSolutionStateME(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< StepperState< Scalar > > &stepperState=Teuchos::null, const Teuchos::RCP< PhysicsState< Scalar > > &physicsState=Teuchos::null)
Nonmember constructor from Thyra ModelEvaluator.
Teuchos::RCP< StepperOperatorSplit< Scalar > > createStepperOperatorSplit(std::vector< Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > > appModels, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.