29 #ifndef RYTHMOS_STEPPER_HELPERS_DEF_HPP 30 #define RYTHMOS_STEPPER_HELPERS_DEF_HPP 32 #include "Rythmos_StepperHelpers_decl.hpp" 33 #include "Rythmos_InterpolationBufferHelpers.hpp" 34 #include "Rythmos_InterpolatorBaseHelpers.hpp" 35 #include "Teuchos_Assert.hpp" 36 #include "Thyra_AssertOp.hpp" 37 #include "Thyra_VectorStdOps.hpp" 42 using Teuchos::ConstNonconstObjectContainer;
44 template<
class Scalar>
45 void assertValidModel(
46 const StepperBase<Scalar>& stepper,
47 const Thyra::ModelEvaluator<Scalar>& model
51 typedef Thyra::ModelEvaluatorBase MEB;
53 TEUCHOS_ASSERT(stepper.acceptsModel());
55 const MEB::InArgs<Scalar> inArgs = model.createInArgs();
56 const MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
59 TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_x));
60 TEUCHOS_ASSERT(outArgs.supports(MEB::OUT_ARG_f));
62 if (stepper.isImplicit()) {
63 TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_x_dot) );
64 TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_alpha) );
65 TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_beta) );
66 TEUCHOS_ASSERT( outArgs.supports(MEB::OUT_ARG_W) );
78 template<
class Scalar>
79 bool setDefaultInitialConditionFromNominalValues(
80 const Thyra::ModelEvaluator<Scalar>& model,
81 const Ptr<StepperBase<Scalar> >& stepper
85 typedef ScalarTraits<Scalar> ST;
86 typedef Thyra::ModelEvaluatorBase MEB;
88 if (isInitialized(*stepper))
91 MEB::InArgs<Scalar> initCond = model.getNominalValues();
93 if (!is_null(initCond.get_x())) {
97 #ifdef HAVE_RYTHMOS_DEBUG 98 THYRA_ASSERT_VEC_SPACES(
"setInitialConditionIfExists(...)",
99 *model.get_x_space(), *initCond.get_x()->space() );
101 if (initCond.supports(MEB::IN_ARG_x_dot)) {
102 if (is_null(initCond.get_x_dot())) {
103 const RCP<Thyra::VectorBase<Scalar> > x_dot =
104 createMember(model.get_x_space());
105 assign(x_dot.ptr(), ST::zero());
108 #ifdef HAVE_RYTHMOS_DEBUG 109 THYRA_ASSERT_VEC_SPACES(
"setInitialConditionIfExists(...)",
110 *model.get_x_space(), *initCond.get_x_dot()->space() );
114 stepper->setInitialCondition(initCond);
125 template<
class Scalar>
126 void restart( StepperBase<Scalar> *stepper )
128 #ifdef HAVE_RYTHMOS_DEBUG 129 TEUCHOS_TEST_FOR_EXCEPT(0==stepper);
130 #endif // HAVE_RYTHMOS_DEBUG 131 typedef Thyra::ModelEvaluatorBase MEB;
133 stepStatus = stepper->getStepStatus();
134 const RCP<const Thyra::ModelEvaluator<Scalar> >
135 model = stepper->getModel();
137 MEB::InArgs<double> initialCondition = model->createInArgs();
138 initialCondition.setArgs(model->getNominalValues());
140 RCP<const Thyra::VectorBase<double> > x, x_dot;
141 Rythmos::get_x_and_x_dot(*stepper,stepStatus.
time,&x,&x_dot);
142 initialCondition.set_x(x);
143 initialCondition.set_x_dot(x_dot);
144 initialCondition.set_t(stepStatus.
time);
147 stepper->setInitialCondition(initialCondition);
150 template<
class Scalar>
151 void eval_model_explicit(
152 const Thyra::ModelEvaluator<Scalar> &model,
153 Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
154 const VectorBase<Scalar>& x_in,
155 const typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag &t_in,
156 const Ptr<VectorBase<Scalar> >& f_out
159 typedef Thyra::ModelEvaluatorBase MEB;
160 MEB::InArgs<Scalar> inArgs = model.createInArgs();
161 MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
162 inArgs.setArgs(basePoint);
163 inArgs.set_x(Teuchos::rcp(&x_in,
false));
164 if (inArgs.supports(MEB::IN_ARG_t)) {
172 if (inArgs.supports(MEB::IN_ARG_x_dot)) {
173 inArgs.set_x_dot(Teuchos::null);
175 outArgs.set_f(Teuchos::rcp(&*f_out,
false));
176 model.evalModel(inArgs,outArgs);
180 #ifdef HAVE_THYRA_ME_POLYNOMIAL 183 template<
class Scalar>
184 void eval_model_explicit_poly(
185 const Thyra::ModelEvaluator<Scalar> &model,
186 Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
187 const Teuchos::Polynomial< VectorBase<Scalar> > &x_poly,
188 const typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag &t,
189 const Ptr<Teuchos::Polynomial<VectorBase<Scalar> > >& f_poly
192 typedef Thyra::ModelEvaluatorBase MEB;
193 MEB::InArgs<Scalar> inArgs = model.createInArgs();
194 MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
195 inArgs.setArgs(basePoint);
196 inArgs.set_x_poly(Teuchos::rcp(&x_poly,
false));
197 if (inArgs.supports(MEB::IN_ARG_t)) {
200 outArgs.set_f_poly(Teuchos::rcp(&*f_poly,
false));
202 model.evalModel(inArgs,outArgs);
206 #endif // HAVE_THYRA_ME_POLYNOMIAL 209 template<
class Scalar>
210 void defaultGetPoints(
212 const Ptr<
const VectorBase<Scalar> >& x_old,
213 const Ptr<
const VectorBase<Scalar> >& xdot_old,
215 const Ptr<
const VectorBase<Scalar> >& x,
216 const Ptr<
const VectorBase<Scalar> >& xdot,
217 const Array<Scalar>& time_vec,
218 const Ptr<Array<Teuchos::RCP<
const Thyra::VectorBase<Scalar> > > >& x_vec,
219 const Ptr<Array<Teuchos::RCP<
const Thyra::VectorBase<Scalar> > > >& xdot_vec,
220 const Ptr<Array<
typename Teuchos::ScalarTraits<Scalar>::magnitudeType> >& accuracy_vec,
221 const Ptr<InterpolatorBase<Scalar> > interpolator
224 typedef Teuchos::ScalarTraits<Scalar> ST;
225 assertTimePointsAreSorted(time_vec);
226 TimeRange<Scalar> tr(t_old, t);
227 TEUCHOS_ASSERT( tr.isValid() );
228 if (!is_null(x_vec)) {
231 if (!is_null(xdot_vec)) {
234 if (!is_null(accuracy_vec)) {
235 accuracy_vec->clear();
237 typename Array<Scalar>::const_iterator time_it = time_vec.begin();
238 RCP<const VectorBase<Scalar> > tmpVec;
239 RCP<const VectorBase<Scalar> > tmpVecDot;
240 for (; time_it != time_vec.end() ; time_it++) {
241 Scalar time = *time_it;
242 asssertInTimeRange(tr, time);
243 Scalar accuracy = ST::zero();
244 if (compareTimeValues(time,t_old)==0) {
245 if (!is_null(x_old)) {
246 tmpVec = x_old->clone_v();
248 if (!is_null(xdot_old)) {
249 tmpVecDot = xdot_old->clone_v();
251 }
else if (compareTimeValues(time,t)==0) {
253 tmpVec = x->clone_v();
255 if (!is_null(xdot)) {
256 tmpVecDot = xdot->clone_v();
259 TEUCHOS_TEST_FOR_EXCEPTION(
260 is_null(interpolator), std::logic_error,
261 "Error, getPoints: This stepper only supports time values on the boundaries!\n" 267 typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
268 typename DataStore<Scalar>::DataStoreVector_t ds_out;
271 DataStore<Scalar> ds;
273 ds.x = rcp(x_old.get(),
false);
274 ds.xdot = rcp(xdot_old.get(),
false);
275 ds_nodes.push_back(ds);
279 DataStore<Scalar> ds;
281 ds.x = rcp(x.get(),
false);
282 ds.xdot = rcp(xdot.get(),
false);
283 ds_nodes.push_back(ds);
285 Array<Scalar> time_vec_in;
286 time_vec_in.push_back(time);
287 interpolate<Scalar>(*interpolator,rcp(&ds_nodes,
false),time_vec_in,&ds_out);
288 Array<Scalar> time_vec_out;
289 Array<RCP<const VectorBase<Scalar> > > x_vec_out;
290 Array<RCP<const VectorBase<Scalar> > > xdot_vec_out;
291 Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> accuracy_vec_out;
292 dataStoreVectorToVector(ds_out,&time_vec_out,&x_vec_out,&xdot_vec_out,&accuracy_vec_out);
293 TEUCHOS_ASSERT( time_vec_out.length()==1 );
294 tmpVec = x_vec_out[0];
295 tmpVecDot = xdot_vec_out[0];
296 accuracy = accuracy_vec_out[0];
298 if (!is_null(x_vec)) {
299 x_vec->push_back(tmpVec);
301 if (!is_null(xdot_vec)) {
302 xdot_vec->push_back(tmpVecDot);
304 if (!is_null(accuracy_vec)) {
305 accuracy_vec->push_back(accuracy);
307 tmpVec = Teuchos::null;
308 tmpVecDot = Teuchos::null;
313 template<
class Scalar>
314 void setStepperModel(
315 const Ptr<StepperBase<Scalar> >& stepper,
316 const RCP<
const Thyra::ModelEvaluator<Scalar> >& model
319 stepper->setModel(model);
322 template<
class Scalar>
323 void setStepperModel(
324 const Ptr<StepperBase<Scalar> >& stepper,
325 const RCP<Thyra::ModelEvaluator<Scalar> >& model
328 stepper->setNonconstModel(model);
331 template<
class Scalar>
332 void setStepperModel(
333 const Ptr<StepperBase<Scalar> >& stepper,
334 ConstNonconstObjectContainer<Thyra::ModelEvaluator<Scalar> >& model
337 if (model.isConst()) {
338 stepper->setModel(model.getConstObj());
341 stepper->setNonconstModel(model.getNonconstObj());
353 #ifdef HAVE_THYRA_ME_POLYNOMIAL 355 #define RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR) \ 356 template void eval_model_explicit_poly( \ 357 const Thyra::ModelEvaluator< SCALAR > &model, \ 358 Thyra::ModelEvaluatorBase::InArgs< SCALAR > &basePoint, \ 359 const Teuchos::Polynomial< VectorBase< SCALAR > > &x_poly, \ 360 const Thyra::ModelEvaluatorBase::InArgs< SCALAR >::ScalarMag &t, \ 361 const Ptr<Teuchos::Polynomial<VectorBase< SCALAR > > >& f_poly \ 364 #else // HAVE_THYRA_ME_POLYNOMIAL 366 #define RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR) 368 #endif // HAVE_THYRA_ME_POLYNOMIAL 371 #define RYTHMOS_STEPPER_HELPERS_INSTANT(SCALAR) \ 373 template void assertValidModel( \ 374 const StepperBase< SCALAR >& stepper, \ 375 const Thyra::ModelEvaluator< SCALAR >& model \ 377 template bool setDefaultInitialConditionFromNominalValues( \ 378 const Thyra::ModelEvaluator< SCALAR >& model, \ 379 const Ptr<StepperBase< SCALAR > >& stepper \ 381 template void restart( StepperBase< SCALAR > *stepper ); \ 383 template void eval_model_explicit( \ 384 const Thyra::ModelEvaluator< SCALAR > &model, \ 385 Thyra::ModelEvaluatorBase::InArgs< SCALAR > &basePoint, \ 386 const VectorBase< SCALAR >& x_in, \ 387 const Thyra::ModelEvaluatorBase::InArgs< SCALAR >::ScalarMag &t_in, \ 388 const Ptr<VectorBase< SCALAR > >& f_out \ 391 RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR) \ 393 template void defaultGetPoints( \ 394 const SCALAR & t_old, \ 395 const Ptr<const VectorBase< SCALAR > >& x_old, \ 396 const Ptr<const VectorBase< SCALAR > >& xdot_old, \ 398 const Ptr<const VectorBase< SCALAR > >& x, \ 399 const Ptr<const VectorBase< SCALAR > >& xdot, \ 400 const Array< SCALAR >& time_vec, \ 401 const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase< SCALAR > > > >& x_vec, \ 402 const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase< SCALAR > > > >& xdot_vec, \ 403 const Ptr<Array<Teuchos::ScalarTraits< SCALAR >::magnitudeType> >& accuracy_vec, \ 404 const Ptr<InterpolatorBase< SCALAR > > interpolator \ 407 template void setStepperModel( \ 408 const Ptr<StepperBase< SCALAR > >& stepper, \ 409 const RCP<const Thyra::ModelEvaluator< SCALAR > >& model \ 412 template void setStepperModel( \ 413 const Ptr<StepperBase< SCALAR > >& stepper, \ 414 const RCP<Thyra::ModelEvaluator< SCALAR > >& model \ 417 template void setStepperModel( \ 418 const Ptr<StepperBase< SCALAR > >& stepper, \ 419 Teuchos::ConstNonconstObjectContainer<Thyra::ModelEvaluator< SCALAR > >& model \ 425 #endif // RYTHMOS_STEPPER_HELPERS_DEF_HPP