Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_ThetaStepper_def.hpp
1//@HEADER
2// ***********************************************************************
3//
4// Rythmos Package
5// Copyright (2006) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// 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
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25//
26// ***********************************************************************
27//@HEADER
28
29#ifndef Rythmos_THETA_STEPPER_DEF_H
30#define Rythmos_THETA_STEPPER_DEF_H
31
32#include "Rythmos_ConfigDefs.h"
33#ifdef HAVE_RYTHMOS_EXPERIMENTAL
34
35#include "Rythmos_ThetaStepper_decl.hpp"
36
37namespace Rythmos {
38
39
44template<class Scalar>
45RCP<ThetaStepper<Scalar> >
46thetaStepper(
47 const RCP<Thyra::ModelEvaluator<Scalar> >& model,
48 const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
49 RCP<Teuchos::ParameterList>& parameterList
50 )
51{
52 Teuchos::RCP<ThetaStepper<Scalar> > stepper =
53 Teuchos::rcp(new ThetaStepper<Scalar>());
54 stepper->setParameterList(parameterList);
55 stepper->setModel(model);
56 stepper->setSolver(solver);
57
58 return stepper;
59}
60
61// ////////////////////////////
62// Defintions
63
64
65// Constructors, intializers, Misc.
66
67
68template<class Scalar>
70{
71 typedef Teuchos::ScalarTraits<Scalar> ST;
72 this->defaultInitializeAll_();
73 Scalar zero = ST::zero();
74 t_ = -ST::one();
75 t_old_ = zero;
76 dt_ = zero;
77 dt_old_ = zero;
78 numSteps_ = 0;
79}
80
81template<class Scalar>
82void ThetaStepper<Scalar>::defaultInitializeAll_()
83{
84 typedef Teuchos::ScalarTraits<Scalar> ST;
85 isInitialized_ = false;
86 haveInitialCondition_ = false;
87 model_ = Teuchos::null;
88 solver_ = Teuchos::null;
89 //basePoint_;
90 x_ = Teuchos::null;
91 x_old_ = Teuchos::null;
92 x_pre_ = Teuchos::null;
93 x_dot_ = Teuchos::null;
94 x_dot_old_ = Teuchos::null;
95 x_dot_really_old_ = Teuchos::null;
96 x_dot_base_ = Teuchos::null;
97 t_ = ST::nan();
98 t_old_ = ST::nan();
99 dt_ = ST::nan();
100 dt_old_ = ST::nan();
101 numSteps_ = -1;
102 thetaStepperType_ = INVALID_THETA_STEPPER_TYPE;
103 theta_ = ST::nan();
104 neModel_ = Teuchos::null;
105 parameterList_ = Teuchos::null;
106 interpolator_ = Teuchos::null;
107 predictor_corrector_begin_after_step_ = -1;
108 default_predictor_order_ = -1;
109}
110
111template<class Scalar>
112bool ThetaStepper<Scalar>::isImplicit() const
113{
114 return true;
115}
116
117
118template<class Scalar>
119void ThetaStepper<Scalar>::setInterpolator(
120 const RCP<InterpolatorBase<Scalar> >& interpolator
121 )
122{
123#ifdef HAVE_RYTHMOS_DEBUG
124 TEUCHOS_TEST_FOR_EXCEPT(is_null(interpolator));
125#endif
126 interpolator_ = interpolator;
127 isInitialized_ = false;
128}
129
130template<class Scalar>
131RCP<InterpolatorBase<Scalar> >
132 ThetaStepper<Scalar>::getNonconstInterpolator()
133{
134 return interpolator_;
135}
136
137template<class Scalar>
138RCP<const InterpolatorBase<Scalar> >
139 ThetaStepper<Scalar>::getInterpolator() const
140{
141 return interpolator_;
142}
143
144
145template<class Scalar>
146RCP<InterpolatorBase<Scalar> >
147ThetaStepper<Scalar>::unSetInterpolator()
148{
149 RCP<InterpolatorBase<Scalar> > temp_interpolator = interpolator_;
150 interpolator_ = Teuchos::null;
151 //isInitialized_ = false;
152 return temp_interpolator;
153}
154
155
156// Overridden from SolverAcceptingStepperBase
157
158
159template<class Scalar>
160void ThetaStepper<Scalar>::setSolver(
161 const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
162 )
163{
164 using Teuchos::as;
165
166 TEUCHOS_TEST_FOR_EXCEPTION(solver == Teuchos::null, std::logic_error,
167 "Error! Thyra::NonlinearSolverBase RCP passed in through ThetaStepper::setSolver is null!"
168 );
169
170 RCP<Teuchos::FancyOStream> out = this->getOStream();
171 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
172 Teuchos::OSTab ostab(out,1,"TS::setSolver");
173 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
174 *out << "solver = " << solver->description() << std::endl;
175 }
176
177 solver_ = solver;
178
179 isInitialized_ = false;
180
181}
182
183
184template<class Scalar>
185RCP<Thyra::NonlinearSolverBase<Scalar> >
186ThetaStepper<Scalar>::getNonconstSolver()
187{
188 return solver_;
189}
190
191
192template<class Scalar>
193RCP<const Thyra::NonlinearSolverBase<Scalar> >
194ThetaStepper<Scalar>::getSolver() const
195{
196 return solver_;
197}
198
199
200// Overridden from StepperBase
201
202
203template<class Scalar>
204bool ThetaStepper<Scalar>::supportsCloning() const
205{
206 return true;
207}
208
209template<class Scalar>
210RCP<StepperBase<Scalar> >
211ThetaStepper<Scalar>::cloneStepperAlgorithm() const
212{
213 RCP<ThetaStepper<Scalar> >
214 stepper = Teuchos::rcp(new ThetaStepper<Scalar>);
215 stepper->isInitialized_ = isInitialized_;
216 stepper->model_ = model_; // Model is stateless so shallow copy is okay!
217
218 if (!is_null(solver_))
219 stepper->solver_ = solver_->cloneNonlinearSolver().assert_not_null();
220
221 stepper->basePoint_ = basePoint_;
222
223 if (!is_null(x_))
224 stepper->x_ = x_->clone_v().assert_not_null();
225 if (!is_null(x_old_))
226 stepper->x_old_ = x_old_->clone_v().assert_not_null();
227 if (!is_null(x_pre_))
228 stepper->x_pre_ = x_pre_->clone_v().assert_not_null();
229
230 if (!is_null(x_dot_))
231 stepper->x_dot_ = x_dot_->clone_v().assert_not_null();
232 if (!is_null(x_dot_old_))
233 stepper->x_dot_old_ = x_dot_old_->clone_v().assert_not_null();
234 if (!is_null(x_dot_really_old_))
235 stepper->x_dot_really_old_ = x_dot_really_old_->clone_v().assert_not_null();
236 if (!is_null(x_dot_base_))
237 stepper->x_dot_base_ = x_dot_base_->clone_v().assert_not_null();
238
239 stepper->t_ = t_;
240 stepper->t_old_ = t_old_;
241
242 stepper->dt_ = dt_;
243 stepper->dt_old_ = dt_old_;
244
245 stepper->numSteps_ = numSteps_;
246
247 stepper->thetaStepperType_ = thetaStepperType_;
248 stepper->theta_ = theta_;
249 stepper->predictor_corrector_begin_after_step_ = predictor_corrector_begin_after_step_;
250 stepper->default_predictor_order_ = default_predictor_order_;
251
252 if (!is_null(neModel_))
253 stepper->neModel_
254 = Teuchos::rcp(new Rythmos::SingleResidualModelEvaluator<Scalar>);
255
256 if (!is_null(parameterList_))
257 stepper->parameterList_ = parameterList(*parameterList_);
258
259 if (!is_null(interpolator_))
260 stepper->interpolator_
261 = interpolator_->cloneInterpolator().assert_not_null(); // ToDo: Implement cloneInterpolator()
262
263 return stepper;
264}
265
266template<class Scalar>
267void ThetaStepper<Scalar>::setModel(
268 const RCP<const Thyra::ModelEvaluator<Scalar> >& model
269 )
270{
271
272 using Teuchos::as;
273
274 TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
275 assertValidModel( *this, *model );
276
277 RCP<Teuchos::FancyOStream> out = this->getOStream();
278 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
279 Teuchos::OSTab ostab(out,1,"TS::setModel");
280 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
281 *out << "model = " << model->description() << std::endl;
282 }
283 model_ = model;
284
285 // Wipe out x. This will either be set thorugh setInitialCondition(...) or
286 // it will be taken from the model's nominal vlaues!
287 x_ = Teuchos::null;
288 x_old_ = Teuchos::null;
289 x_pre_ = Teuchos::null;
290
291 x_dot_ = Teuchos::null;
292 x_dot_old_ = Teuchos::null;
293 x_dot_really_old_ = Teuchos::null;
294 x_dot_base_ = Teuchos::null;
295
296 isInitialized_ = false;
297 haveInitialCondition_ = setDefaultInitialConditionFromNominalValues<Scalar>(
298 *model_, Teuchos::ptr(this) );
299
300}
301
302template<class Scalar>
303void ThetaStepper<Scalar>::setNonconstModel(
304 const RCP<Thyra::ModelEvaluator<Scalar> >& model
305 )
306{
307 this->setModel(model); // TODO 09/09/09 tscoffe: use ConstNonconstObjectContainer!
308}
309
310
311template<class Scalar>
312RCP<const Thyra::ModelEvaluator<Scalar> >
313ThetaStepper<Scalar>::getModel() const
314{
315 return model_;
316}
317
318
319template<class Scalar>
320RCP<Thyra::ModelEvaluator<Scalar> >
321ThetaStepper<Scalar>::getNonconstModel()
322{
323 TEUCHOS_TEST_FOR_EXCEPT(true);
324 return Teuchos::null;
325}
326
327
328template<class Scalar>
329void ThetaStepper<Scalar>::setInitialCondition(
330 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
331 )
332{
333
334 typedef Teuchos::ScalarTraits<Scalar> ST;
335
336 basePoint_ = initialCondition;
337
338 // x
339
340 RCP<const Thyra::VectorBase<Scalar> >
341 x_init = initialCondition.get_x();
342
343#ifdef HAVE_RYTHMOS_DEBUG
344 TEUCHOS_TEST_FOR_EXCEPTION(
345 is_null(x_init), std::logic_error,
346 "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
347 "then x can not be null!" );
348#endif
349
350 x_ = x_init->clone_v();
351
352 // x_dot
353
354 RCP<const Thyra::VectorBase<Scalar> >
355 x_dot_init = initialCondition.get_x_dot();
356
357 if (!is_null(x_dot_init)) {
358 x_dot_ = x_dot_init->clone_v();
359 }
360 else {
361 x_dot_ = createMember(x_->space());
362 assign(x_dot_.ptr(),ST::zero());
363 }
364
365 // t
366
367 t_ = initialCondition.get_t();
368
369 t_old_ = t_;
370
371 dt_old_ = 0.0;
372
373 // x pre
374
375 x_pre_ = x_->clone_v();
376
377 // x old
378
379 x_old_ = x_->clone_v();
380
381 // x dot base
382
383 x_dot_base_ = x_->clone_v();
384
385 // x dot old
386
387 x_dot_old_ = x_dot_->clone_v();
388
389 // x dot really old
390
391 x_dot_really_old_ = x_dot_->clone_v();
392
393 haveInitialCondition_ = true;
394
395}
396
397
398template<class Scalar>
399Thyra::ModelEvaluatorBase::InArgs<Scalar>
400ThetaStepper<Scalar>::getInitialCondition() const
401{
402 return basePoint_;
403}
404
405
406template<class Scalar>
407Scalar ThetaStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepSizeType)
408{
409
410 using Teuchos::as;
411 using Teuchos::incrVerbLevel;
412 typedef Teuchos::ScalarTraits<Scalar> ST;
413 typedef Thyra::NonlinearSolverBase<Scalar> NSB;
414 typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
415
416 initialize_();
417
418 // DEBUG
419 //this->setOverridingVerbLevel(Teuchos::VERB_EXTREME);
420 //this->setVerbLevel(Teuchos::VERB_EXTREME);
421
422 RCP<Teuchos::FancyOStream> out = this->getOStream();
423 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
424 Teuchos::OSTab ostab(out,1,"TS::takeStep");
425 VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
426
427 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
428 *out
429 << "\nEntering " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
430 << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n";
431 }
432
433 dt_ = dt;
434 V_StV( x_old_.ptr(), Scalar(ST::one()), *x_ );
435 V_StV( x_dot_really_old_.ptr(), Scalar(ST::one()), *x_dot_old_ );
436 V_StV( x_dot_old_.ptr(), Scalar(ST::one()), *x_dot_ );
437
438 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
439 *out << "\nSetting dt_ and old data ..." << std::endl;
440 *out << "\ndt_ = " << dt_;
441 *out << "\nx_old_ = " << *x_old_;
442 *out << "\nx_dot_old_ = " << *x_dot_old_;
443 *out << "\nx_dot_really_old_ = " << *x_dot_really_old_;
444 }
445
446 if ((stepSizeType == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
447 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
448 *out << "\nThe arguments to takeStep are not valid for ThetaStepper at this time." << std::endl;
449 return(Scalar(-ST::one()));
450 }
451 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
452 *out << "\ndt = " << dt << std::endl;
453 *out << "\nstepSizeType = " << stepSizeType << std::endl;
454 }
455
456 // compute predictor
457 obtainPredictor_();
458
459 //
460 // Setup the nonlinear equations:
461 //
462 // substitute:
463 //
464 // x_dot = ( 1/(theta*dt) )*x + ( -1/(theta*dt) )*x_old + ( -(1-theta)/theta )*x_dot_old
465 //
466 // f( x_dot, x, t ) = 0
467 //
468
469 const double theta = (numSteps_+1>=predictor_corrector_begin_after_step_) ? theta_ : 1.0;
470
471 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
472 *out << "\nSetting x_dot_base_ ..." << std::endl;
473 *out << "\ntheta = " << theta;
474 *out << "\nx_ = " << *x_;
475 *out << "\nx_dot_old_ = " << *x_dot_old_;
476 }
477
478 const Scalar coeff_x_dot = Scalar(ST::one()/(theta*dt));
479 const Scalar coeff_x = ST::one();
480
481 const Scalar x_coeff = Scalar(-coeff_x_dot);
482 const Scalar x_dot_old_coeff = Scalar( -(ST::one()-theta)/theta);
483
484 V_StV( x_dot_base_.ptr(), x_coeff, *x_old_ );
485 Vp_StV( x_dot_base_.ptr(), x_dot_old_coeff, *x_dot_old_);
486
487 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
488 *out << "\nx_dot_base_ = " << *x_dot_base_;
489 }
490
491 if(!neModel_.get()) {
492 neModel_ = Teuchos::rcp(new Rythmos::SingleResidualModelEvaluator<Scalar>());
493 }
494
495 neModel_->initializeSingleResidualModel(
496 model_,
497 basePoint_,
498 coeff_x_dot,
499 x_dot_base_,
500 coeff_x,
501 Teuchos::null, // x_base
502 t_+dt, // t_base
503 Teuchos::null // x_bar_init
504 );
505 if( solver_->getModel().get() != neModel_.get() ) {
506 solver_->setModel(neModel_);
507 }
508
509 solver_->setVerbLevel(this->getVerbLevel());
510
511 //
512 // Solve the implicit nonlinear system to a tolerance of ???
513 //
514
515 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
516 *out << "\nSolving the implicit theta-stepper timestep equation"
517 << " with theta = " << theta << "\n";
518 }
519
520 Thyra::SolveStatus<Scalar>
521 neSolveStatus = solver_->solve(&*x_);
522
523 // In the above solve, on input *x_ is the old value of x for the previous
524 // time step which is used as the initial guess for the solver. On output,
525 // *x_ is the converged timestep solution.
526
527 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
528 *out << "\nOutput status of nonlinear solve:\n" << neSolveStatus;
529 }
530
531 // 2007/05/18: rabartl: ToDo: Above, get the solve status from the above
532 // solve and at least print warning message if the solve fails! Actually,
533 // you should most likely thrown an exception if the solve fails or return
534 // false if appropriate
535
536 //
537 // Update the step data
538 //
539
540 // x_dot = ( 1/(theta*dt) )*x + ( -1/(theta*dt) )*x_old + ( -(1-theta)/theta )*x_dot_old
541
542 V_StV( x_dot_.ptr(), Scalar( ST::one()/(theta*dt)), *x_ );
543 Vp_StV( x_dot_.ptr(), Scalar(-ST::one()/(theta*dt)), *x_old_ );
544 Vp_StV( x_dot_.ptr(), Scalar( -(ST::one()-theta)/theta), *x_dot_old_ );
545
546 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
547 *out << "\nUpdating x_dot_ ...\n";
548 *out << "\nx_dot_ = " << *x_dot_ << std::endl;
549 }
550
551 t_old_ = t_;
552 dt_old_ = dt;
553 t_ += dt;
554
555 numSteps_++;
556
557 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
558 *out << "\nt_old_ = " << t_old_ << std::endl;
559 *out << "\nt_ = " << t_ << std::endl;
560 }
561
562#ifdef HAVE_RYTHMOS_DEBUG
563
564 if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
565 *out << "\nChecking to make sure that solution and the interpolated solution are the same! ...\n";
566
567 {
568
569 using SMT = ScalarTraits<ScalarMag>;
570
571 Teuchos::OSTab tab(out);
572
573 const StepStatus<Scalar> stepStatus = this->getStepStatus();
574
575 RCP<const Thyra::VectorBase<Scalar> >
576 x = stepStatus.solution,
577 xdot = stepStatus.solutionDot;
578
579 Array<Scalar> time_vec = Teuchos::tuple(stepStatus.time);
580 Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
581 this->getPoints(time_vec,&x_vec,&xdot_vec,0);
582
583 RCP<const Thyra::VectorBase<Scalar> >
584 x_interp = x_vec[0],
585 xdot_interp = xdot_vec[0];
586
587 TEUCHOS_TEST_FOR_EXCEPT(
588 !Thyra::testRelNormDiffErr(
589 "x", *x, "x_interp", *x_interp,
590 "2*epsilon", ScalarMag(100.0*SMT::eps()),
591 "2*epsilon", ScalarMag(100.0*SMT::eps()),
592 includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
593 )
594 );
595
596 TEUCHOS_TEST_FOR_EXCEPT(
597 !Thyra::testRelNormDiffErr(
598 "xdot", *xdot, "xdot_interp", *xdot_interp,
599 "2*epsilon", ScalarMag(100.0*SMT::eps()),
600 "2*epsilon", ScalarMag(100.0*SMT::eps()),
601 includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
602 )
603 );
604
605 }
606
607 // 2007/07/25: rabartl: ToDo: Move the above test into a helper function so
608 // that it can be used from lots of different places!
609
610#endif // HAVE_RYTHMOS_DEBUG
611
612 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
613 *out
614 << "\nLeaving " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
615 << "::takeStep(...) ...\n";
616 }
617
618 return(dt);
619
620}
621
622
623template<class Scalar>
624const StepStatus<Scalar> ThetaStepper<Scalar>::getStepStatus() const
625{
626
627 StepStatus<Scalar> stepStatus; // Defaults to unknown status
628
629 if (!isInitialized_) {
630 stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
631 }
632 else if (numSteps_ > 0) {
633 stepStatus.stepStatus = STEP_STATUS_CONVERGED;
634 }
635 // else unknown
636
637 stepStatus.stepSize = dt_;
638 stepStatus.order = 1;
639 stepStatus.time = t_;
640 stepStatus.solution = x_;
641 stepStatus.solutionDot = x_dot_;
642
643 return(stepStatus);
644
645}
646
647
648// Overridden from InterpolationBufferBase
649
650
651template<class Scalar>
652RCP<const Thyra::VectorSpaceBase<Scalar> >
653ThetaStepper<Scalar>::get_x_space() const
654{
655 return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
656}
657
658
659template<class Scalar>
660void ThetaStepper<Scalar>::addPoints(
661 const Array<Scalar>& time_vec,
662 const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
663 const Array<RCP<const Thyra::VectorBase<Scalar> > >& /* xdot_vec */
664 )
665{
666
667 typedef Teuchos::ScalarTraits<Scalar> ST;
668 using Teuchos::as;
669
670 initialize_();
671
672#ifdef HAVE_RYTHMOS_DEBUG
673 TEUCHOS_TEST_FOR_EXCEPTION(
674 time_vec.size() == 0, std::logic_error,
675 "Error, addPoints called with an empty time_vec array!\n");
676#endif // HAVE_RYTHMOS_DEBUG
677
678 RCP<Teuchos::FancyOStream> out = this->getOStream();
679 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
680 Teuchos::OSTab ostab(out,1,"TS::setPoints");
681
682 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
683 *out << "time_vec = " << std::endl;
684 for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
685 *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
686 }
687 }
688 else if (time_vec.size() == 1) {
689 int n = 0;
690 t_ = time_vec[n];
691 t_old_ = t_;
692 Thyra::V_V(x_.ptr(),*x_vec[n]);
693 Thyra::V_V(x_dot_base_.ptr(),*x_);
694 }
695 else {
696 int n = time_vec.size()-1;
697 int nm1 = time_vec.size()-2;
698 t_ = time_vec[n];
699 t_old_ = time_vec[nm1];
700 Thyra::V_V(x_.ptr(),*x_vec[n]);
701 Scalar dt = t_ - t_old_;
702 Thyra::V_StV(x_dot_base_.ptr(),Scalar(-ST::one()/dt),*x_vec[nm1]);
703 }
704}
705
706
707template<class Scalar>
708TimeRange<Scalar> ThetaStepper<Scalar>::getTimeRange() const
709{
710 if ( !isInitialized_ && haveInitialCondition_ )
711 return timeRange<Scalar>(t_,t_);
712 if ( !isInitialized_ && !haveInitialCondition_ )
713 return invalidTimeRange<Scalar>();
714 return timeRange<Scalar>(t_old_,t_);
715}
716
717
718template<class Scalar>
719void ThetaStepper<Scalar>::getPoints(
720 const Array<Scalar>& time_vec,
721 Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
722 Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
723 Array<ScalarMag>* accuracy_vec
724 ) const
725{
726 using Teuchos::constOptInArg;
727 using Teuchos::ptr;
728 typedef Teuchos::ScalarTraits<Scalar> ST;
729
730 TEUCHOS_ASSERT(haveInitialCondition_);
731
732 RCP<Thyra::VectorBase<Scalar> > x_temp = x_;
733 if (compareTimeValues(t_old_,t_)!=0) {
734 Scalar dt = t_ - t_old_;
735 x_temp = x_dot_base_->clone_v();
736 Thyra::Vt_S(x_temp.ptr(),Scalar(-ST::one()*dt)); // undo the scaling
737 }
738 defaultGetPoints<Scalar>(
739 t_old_, constOptInArg(*x_temp), constOptInArg(*x_dot_old_),
740 t_, constOptInArg(*x_), constOptInArg(*x_dot_),
741 time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
742 ptr(interpolator_.get())
743 );
744
745 /*
746 using Teuchos::as;
747 typedef Teuchos::ScalarTraits<Scalar> ST;
748 typedef typename ST::magnitudeType ScalarMag;
749 typedef Teuchos::ScalarTraits<ScalarMag> SMT;
750 typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
751 typename DataStore<Scalar>::DataStoreVector_t ds_out;
752
753#ifdef HAVE_RYTHMOS_DEBUG
754 TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
755 TEUCHOS_TEST_FOR_EXCEPT( 0 == x_vec );
756#endif
757
758 RCP<Teuchos::FancyOStream> out = this->getOStream();
759 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
760 Teuchos::OSTab ostab(out,1,"TS::getPoints");
761 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
762 *out
763 << "\nEntering " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
764 << "::getPoints(...) ...\n";
765 }
766 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
767 for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
768 *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
769 }
770 *out << "I can interpolate in the interval [" << t_old_ << "," << t_ << "]." << std::endl;
771 }
772
773 if (t_old_ != t_) {
774 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
775 *out << "Passing two points to interpolator: " << t_old_ << " and " << t_ << std::endl;
776 }
777 DataStore<Scalar> ds_temp;
778 Scalar dt = t_ - t_old_;
779#ifdef HAVE_RYTHMOS_DEBUG
780 TEUCHOS_TEST_FOR_EXCEPT(
781 !Thyra::testRelErr(
782 "dt", dt, "dt_", dt_,
783 "1e+4*epsilon", ScalarMag(1e+4*SMT::eps()),
784 "1e+2*epsilon", ScalarMag(1e+2*SMT::eps()),
785 as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) ? out.get() : 0
786 )
787 );
788#endif
789 ds_temp.time = t_old_;
790 ds_temp.x = x_old_;
791 ds_temp.xdot = x_dot_old_;
792 ds_temp.accuracy = ScalarMag(dt);
793 ds_nodes.push_back(ds_temp);
794 ds_temp.time = t_;
795 ds_temp.x = x_;
796 ds_temp.xdot = x_dot_;
797 ds_temp.accuracy = ScalarMag(dt);
798 ds_nodes.push_back(ds_temp);
799 }
800 else {
801 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
802 *out << "Passing one point to interpolator: " << t_ << std::endl;
803 }
804 DataStore<Scalar> ds_temp;
805 ds_temp.time = t_;
806 ds_temp.x = x_;
807 ds_temp.xdot = x_dot_;
808 ds_temp.accuracy = ScalarMag(ST::zero());
809 ds_nodes.push_back(ds_temp);
810 }
811 interpolate<Scalar>(*interpolator_,rcp(&ds_nodes,false),time_vec,&ds_out);
812 Array<Scalar> time_out;
813 dataStoreVectorToVector(ds_out,&time_out,x_vec,xdot_vec,accuracy_vec);
814 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
815 *out << "Passing out the interpolated values:" << std::endl;
816 for (int i=0; i<Teuchos::as<int>(time_out.size()) ; ++i) {
817 if (x_vec) {
818 if ( (*x_vec)[i] == Teuchos::null) {
819 *out << "x_vec[" << i << "] = Teuchos::null" << std::endl;
820 }
821 else {
822 *out << "time[" << i << "] = " << time_out[i] << std::endl;
823 *out << "x_vec[" << i << "] = " << std::endl;
824 (*x_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
825 }
826 }
827 if (xdot_vec) {
828 if ( (*xdot_vec)[i] == Teuchos::null) {
829 *out << "xdot_vec[" << i << "] = Teuchos::null" << std::endl;
830 }
831 else {
832 *out << "xdot_vec[" << i << "] = " << std::endl;
833 (*xdot_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
834 }
835 }
836 if(accuracy_vec) {
837 *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
838 }
839 }
840 }
841 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
842 *out
843 << "Leaving " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
844 << "::getPoints(...) ...\n";
845 }
846 */
847
848}
849
850
851template<class Scalar>
852void ThetaStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
853{
854 using Teuchos::as;
855
856 TEUCHOS_ASSERT( time_vec != NULL );
857
858 time_vec->clear();
859 if (!haveInitialCondition_) {
860 return;
861 }
862
863 time_vec->push_back(t_old_);
864 if (numSteps_ > 0) {
865 time_vec->push_back(t_);
866 }
867
868 RCP<Teuchos::FancyOStream> out = this->getOStream();
869 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
870 Teuchos::OSTab ostab(out,1,"TS::getNodes");
871 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
872 *out << this->description() << std::endl;
873 for (int i=0 ; i<Teuchos::as<int>(time_vec->size()) ; ++i) {
874 *out << "time_vec[" << i << "] = " << (*time_vec)[i] << std::endl;
875 }
876 }
877}
878
879
880template<class Scalar>
881void ThetaStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
882{
883 initialize_();
884 using Teuchos::as;
885 RCP<Teuchos::FancyOStream> out = this->getOStream();
886 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
887 Teuchos::OSTab ostab(out,1,"TS::removeNodes");
888 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
889 *out << "time_vec = " << std::endl;
890 for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
891 *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
892 }
893 }
894 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ThetaStepper at this time.\n");
895 // TODO:
896 // if any time in time_vec matches t_ or t_old_, then do the following:
897 // remove t_old_: set t_old_ = t_ and set x_dot_base_ = x_
898 // remove t_: set t_ = t_old_ and set x_ = -dt*x_dot_base_
899}
900
901
902template<class Scalar>
903int ThetaStepper<Scalar>::getOrder() const
904{
905 return (thetaStepperType_==ImplicitEuler) ? 1 : 2;
906}
907
908
909// Overridden from Teuchos::ParameterListAcceptor
910
911
912template <class Scalar>
913void ThetaStepper<Scalar>::setParameterList(
914 RCP<Teuchos::ParameterList> const& paramList
915 )
916{
917 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
918 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
919 parameterList_ = paramList;
920 Teuchos::readVerboseObjectSublist(&*parameterList_,this);
921
922 RCP<ParameterList> pl_theta = Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
923
924 std::string thetaStepperTypeString =
925 Teuchos::getParameter<std::string>(*pl_theta, ThetaStepperType_name);
926
927 if (thetaStepperTypeString == "Implicit Euler")
928 thetaStepperType_ = ImplicitEuler;
929 else if (thetaStepperTypeString == "Trapezoid")
930 thetaStepperType_ = Trapezoid;
931 else
932 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
933 "Value of " << ThetaStepperType_name << " = " << thetaStepperTypeString
934 << " is invalid for Rythmos::ThetaStepper");
935
936 default_predictor_order_ =
937 Teuchos::getParameter<int>(*pl_theta, PredictorOrder_name);
938
939 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
940 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
941 if ( Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
942 *out << ThetaStepperType_name << " = " << thetaStepperTypeString << std::endl;
943 }
944}
945
946
947template <class Scalar>
948RCP<Teuchos::ParameterList>
949ThetaStepper<Scalar>::getNonconstParameterList()
950{
951 return(parameterList_);
952}
953
954
955template <class Scalar>
956RCP<Teuchos::ParameterList>
957ThetaStepper<Scalar>::unsetParameterList()
958{
959 RCP<Teuchos::ParameterList>
960 temp_param_list = parameterList_;
961 parameterList_ = Teuchos::null;
962 return(temp_param_list);
963}
964
965
966template<class Scalar>
967RCP<const Teuchos::ParameterList>
968ThetaStepper<Scalar>::getValidParameters() const
969{
970 using Teuchos::ParameterList;
971
972 static RCP<const ParameterList> validPL;
973
974 if (is_null(validPL)) {
975
976 RCP<ParameterList> pl_top_level = Teuchos::parameterList();
977
978 RCP<ParameterList> pl = Teuchos::sublist(pl_top_level, RythmosStepControlSettings_name);
979
980 pl->set<std::string> ( ThetaStepperType_name, ThetaStepperType_default,
981 "Name of Stepper Type in Theta Stepper"
982 );
983
984 pl->set<int> ( PredictorOrder_name, PredictorOrder_default,
985 "Order of Predictor in Theta Stepper, can be 0, 1, 2"
986 );
987
988 Teuchos::setupVerboseObjectSublist(&*pl_top_level);
989 validPL = pl_top_level;
990 }
991 return validPL;
992}
993
994
995// Overridden from Teuchos::Describable
996
997
998template<class Scalar>
999void ThetaStepper<Scalar>::describe(
1000 Teuchos::FancyOStream &out,
1001 const Teuchos::EVerbosityLevel verbLevel
1002 ) const
1003{
1004 using Teuchos::as;
1005 Teuchos::OSTab tab(out);
1006 if (!isInitialized_) {
1007 out << this->description() << " : This stepper is not initialized yet" << std::endl;
1008 return;
1009 }
1010 if (
1011 as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
1012 ||
1013 as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
1014 )
1015 {
1016 out << this->description() << "::describe:" << std::endl;
1017 out << "model = " << model_->description() << std::endl;
1018 out << "solver = " << solver_->description() << std::endl;
1019 if (neModel_ == Teuchos::null) {
1020 out << "neModel = Teuchos::null" << std::endl;
1021 } else {
1022 out << "neModel = " << neModel_->description() << std::endl;
1023 }
1024 }
1025 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
1026 out << "t_ = " << t_ << std::endl;
1027 out << "t_old_ = " << t_old_ << std::endl;
1028 }
1029 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
1030 }
1031 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
1032 out << "model_ = " << std::endl;
1033 model_->describe(out,verbLevel);
1034 out << "solver_ = " << std::endl;
1035 solver_->describe(out,verbLevel);
1036 if (neModel_ == Teuchos::null) {
1037 out << "neModel = Teuchos::null" << std::endl;
1038 } else {
1039 out << "neModel = " << std::endl;
1040 neModel_->describe(out,verbLevel);
1041 }
1042 out << "x_ = " << std::endl;
1043 x_->describe(out,verbLevel);
1044 out << "x_dot_base_ = " << std::endl;
1045 x_dot_base_->describe(out,verbLevel);
1046 }
1047}
1048
1049
1050// private
1051
1052
1053template <class Scalar>
1054void ThetaStepper<Scalar>::initialize_()
1055{
1056
1057 if (isInitialized_)
1058 return;
1059
1060 TEUCHOS_TEST_FOR_EXCEPT(is_null(model_));
1061 TEUCHOS_TEST_FOR_EXCEPT(is_null(solver_));
1062 TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
1063
1064#ifdef HAVE_RYTHMOS_DEBUG
1065 THYRA_ASSERT_VEC_SPACES(
1066 "Rythmos::ThetaStepper::setInitialCondition(...)",
1067 *x_->space(), *model_->get_x_space() );
1068#endif // HAVE_RYTHMOS_DEBUG
1069
1070 if ( is_null(interpolator_) ) {
1071 // If an interpolator has not been explicitly set, then just create
1072 // a default linear interpolator.
1073 interpolator_ = Teuchos::rcp(new LinearInterpolator<Scalar> );
1074 // 2007/05/18: rabartl: ToDo: Replace this with a Hermete interplator
1075 // when it is implementated!
1076 }
1077
1078 if (thetaStepperType_ == ImplicitEuler)
1079 {
1080 theta_ = 1.0;
1081 predictor_corrector_begin_after_step_ = 2;
1082 }
1083 else
1084 {
1085 theta_ = 0.5;
1086 predictor_corrector_begin_after_step_ = 3;
1087 }
1088
1089 isInitialized_ = true;
1090
1091}
1092
1093template<class Scalar>
1094void ThetaStepper<Scalar>::obtainPredictor_()
1095{
1096
1097 using Teuchos::as;
1098 typedef Teuchos::ScalarTraits<Scalar> ST;
1099
1100 if (!isInitialized_) {
1101 return;
1102 }
1103
1104 RCP<Teuchos::FancyOStream> out = this->getOStream();
1105 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1106 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1107 *out << "Obtaining predictor..." << std::endl;
1108 }
1109
1110 const int preferred_predictor_order = std::min(default_predictor_order_, thetaStepperType_ + 1);
1111 const int max_predictor_order_at_this_timestep = std::max(0, numSteps_);
1112
1113 const int predictor_order = std::min(preferred_predictor_order, max_predictor_order_at_this_timestep);
1114
1115 switch (predictor_order)
1116 {
1117 case 0:
1118 V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1119 break;
1120 case 1:
1121 {
1122 V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1123
1124 TEUCHOS_TEST_FOR_EXCEPT (dt_ <= 0.0);
1125
1126 Vp_StV(x_pre_.ptr(), dt_, *x_dot_old_);
1127 }
1128 break;
1129 case 2:
1130 {
1131 V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1132
1133 TEUCHOS_TEST_FOR_EXCEPT (dt_ <= 0.0);
1134 TEUCHOS_TEST_FOR_EXCEPT (dt_old_ <= 0.0);
1135
1136 const Scalar coeff_x_dot_old = (0.5 * dt_) * (2.0 + dt_/dt_old_);
1137 const Scalar coeff_x_dot_really_old = - (0.5 * dt_) * (dt_/dt_old_);
1138
1139 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1140 *out << "x_dot_old_ = " << *x_dot_old_ << std::endl;
1141 *out << "x_dot_really_old_ = " << *x_dot_really_old_ << std::endl;
1142 }
1143
1144 Vp_StV( x_pre_.ptr(), coeff_x_dot_old, *x_dot_old_);
1145 Vp_StV( x_pre_.ptr(), coeff_x_dot_really_old, *x_dot_really_old_);
1146 }
1147 break;
1148 default:
1149 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
1150 "Invalid predictor order " << predictor_order << ". Valid values are 0, 1, and 2.");
1151 }
1152
1153 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1154 *out << "x_pre_ = " << *x_pre_ << std::endl;
1155 }
1156
1157 // copy to current solution
1158 V_StV(x_.ptr(), Scalar(ST::one()), *x_pre_);
1159}
1160
1161//
1162// Explicit Instantiation macro
1163//
1164// Must be expanded from within the Rythmos namespace!
1165//
1166
1167#define RYTHMOS_THETA_STEPPER_INSTANT(SCALAR) \
1168 \
1169 template class ThetaStepper< SCALAR >; \
1170 \
1171 template RCP< ThetaStepper< SCALAR > > \
1172 thetaStepper( \
1173 const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
1174 const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \
1175 RCP<Teuchos::ParameterList>& parameterList \
1176 );
1177
1178
1179} // namespace Rythmos
1180
1181#endif // HAVE_RYTHMOS_EXPERIMENTAL
1182
1183#endif //Rythmos_THETA_STEPPER_DEF_H