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_ =
127 Teuchos::rcp(new StepperOperatorSplitModifierDefault<Scalar>());
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
StepperOperatorSplitAppAction class for StepperOperatorSplit.
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),...
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.