30 #ifndef RYTHMOS_DIAGONALIMPLICITRK_MODEL_EVALUATOR_HPP 31 #define RYTHMOS_DIAGONALIMPLICITRK_MODEL_EVALUATOR_HPP 34 #include "Rythmos_Types.hpp" 35 #include "Rythmos_RKButcherTableauHelpers.hpp" 36 #include "Thyra_StateFuncModelEvaluatorBase.hpp" 37 #include "Thyra_ModelEvaluatorHelpers.hpp" 38 #include "Thyra_ModelEvaluatorDelegatorBase.hpp" 39 #include "Thyra_DefaultProductVectorSpace.hpp" 40 #include "Thyra_DefaultBlockedLinearOp.hpp" 41 #include "Thyra_VectorStdOps.hpp" 42 #include "Thyra_TestingTools.hpp" 43 #include "Teuchos_as.hpp" 44 #include "Teuchos_SerialDenseMatrix.hpp" 45 #include "Teuchos_SerialDenseVector.hpp" 46 #include "Teuchos_Assert.hpp" 53 template<
class Scalar>
66 const RCP<
const Thyra::ModelEvaluator<Scalar> >& daeModel,
67 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& basePoint,
68 const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >& dirk_W_factory,
69 const RCP<
const RKButcherTableauBase<Scalar> >& irkButcherTableau
74 const RCP<
const Thyra::VectorBase<Scalar> > &x_old,
82 const Thyra::VectorBase<Scalar>& stage_solution
96 RCP<const Thyra::VectorSpaceBase<Scalar> >
get_x_space()
const;
98 RCP<const Thyra::VectorSpaceBase<Scalar> >
get_f_space()
const;
100 RCP<Thyra::LinearOpBase<Scalar> >
create_W_op()
const;
102 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
get_W_factory()
const;
106 Thyra::ModelEvaluatorBase::InArgs<Scalar>
createInArgs()
const;
116 Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl()
const;
119 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
120 const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
128 RCP<const Thyra::ModelEvaluator<Scalar> > daeModel_;
129 Thyra::ModelEvaluatorBase::InArgs<Scalar> basePoint_;
130 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > dirk_W_factory_;
131 RCP<const RKButcherTableauBase<Scalar> > dirkButcherTableau_;
133 Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs_;
134 Thyra::ModelEvaluatorBase::OutArgs<Scalar> outArgs_;
135 Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
137 RCP<const Thyra::ProductVectorSpaceBase<Scalar> > stage_space_;
138 RCP<Thyra::ProductVectorBase<Scalar> > stage_derivatives_;
140 bool setTimeStepPointCalled_;
141 RCP<const Thyra::VectorBase<Scalar> > x_old_;
155 template<
class Scalar>
156 RCP<DiagonalImplicitRKModelEvaluator<Scalar> >
158 const RCP<
const Thyra::ModelEvaluator<Scalar> >& daeModel,
159 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
160 const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &dirk_W_factory,
161 const RCP<
const RKButcherTableauBase<Scalar> > &irkButcherTableau
164 RCP<DiagonalImplicitRKModelEvaluator<Scalar> >
166 dirkModel->initializeDIRKModel(daeModel,basePoint,dirk_W_factory,irkButcherTableau);
178 template<
class Scalar>
181 setTimeStepPointCalled_ =
false;
182 isInitialized_ =
false;
190 template<
class Scalar>
192 const RCP<
const Thyra::ModelEvaluator<Scalar> >& daeModel,
193 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& basePoint,
194 const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >& dirk_W_factory,
195 const RCP<
const RKButcherTableauBase<Scalar> >& irkButcherTableau
198 typedef ScalarTraits<Scalar> ST;
201 TEUCHOS_TEST_FOR_EXCEPTION(
202 is_null(basePoint.get_x()),
204 "Error! The basepoint x vector is null!" 206 TEUCHOS_TEST_FOR_EXCEPTION(
209 "Error! The model evaluator pointer is null!" 211 TEUCHOS_TEST_FOR_EXCEPTION(
212 !daeModel->get_x_space()->isCompatible(*(basePoint.get_x()->space())),
214 "Error! The basepoint input arguments are incompatible with the model evaluator vector space!" 218 daeModel_ = daeModel;
219 basePoint_ = basePoint;
220 dirk_W_factory_ = dirk_W_factory;
221 dirkButcherTableau_ = irkButcherTableau;
223 const int numStages = dirkButcherTableau_->numStages();
225 using Teuchos::rcp_dynamic_cast;
226 stage_space_ = productVectorSpace(daeModel_->get_x_space(),numStages);
227 RCP<const Thyra::VectorSpaceBase<Scalar> > vs = rcp_dynamic_cast<
const Thyra::VectorSpaceBase<Scalar> >(stage_space_,
true);
228 stage_derivatives_ = rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(createMember(vs),
true);
229 Thyra::V_S( rcp_dynamic_cast<Thyra::VectorBase<Scalar> >(stage_derivatives_).ptr(),ST::zero());
233 typedef Thyra::ModelEvaluatorBase MEB;
234 MEB::InArgsSetup<Scalar> inArgs;
235 inArgs.setModelEvalDescription(this->description());
236 inArgs.setSupports(MEB::IN_ARG_x);
241 typedef Thyra::ModelEvaluatorBase MEB;
242 MEB::OutArgsSetup<Scalar> outArgs;
243 outArgs.setModelEvalDescription(this->description());
244 outArgs.setSupports(MEB::OUT_ARG_f);
245 outArgs.setSupports(MEB::OUT_ARG_W_op);
249 nominalValues_ = inArgs_;
251 isInitialized_ =
true;
255 template<
class Scalar>
257 const RCP<
const Thyra::VectorBase<Scalar> > &x_old,
259 const Scalar &delta_t
262 typedef ScalarTraits<Scalar> ST;
263 TEUCHOS_TEST_FOR_EXCEPT( is_null(x_old) );
264 TEUCHOS_TEST_FOR_EXCEPT( delta_t <= ST::zero() );
266 TEUCHOS_TEST_FOR_EXCEPTION(!daeModel_->get_x_space()->isCompatible(*(x_old->space())), std::logic_error,
267 "Error! The incoming VectorBase object is not compatible with the DAE model that was provided at initialization!");
271 setTimeStepPointCalled_ =
true;
274 template<
class Scalar>
277 const Thyra::VectorBase<Scalar>& stage_solution
280 TEUCHOS_TEST_FOR_EXCEPT(stageNumber < 0);
281 TEUCHOS_TEST_FOR_EXCEPT(stageNumber >= dirkButcherTableau_->numStages());
282 Thyra::V_V(stage_derivatives_->getNonconstVectorBlock(stageNumber).ptr(),stage_solution);
285 template<
class Scalar>
290 TEUCHOS_TEST_FOR_EXCEPT(currentStage < 0);
291 TEUCHOS_TEST_FOR_EXCEPT(currentStage >= dirkButcherTableau_->numStages());
292 currentStage_ = currentStage;
300 template<
class Scalar>
301 RCP<const Thyra::VectorSpaceBase<Scalar> >
304 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
305 "Error, initializeDIRKModel must be called first!\n" 307 return daeModel_->get_x_space();
311 template<
class Scalar>
312 RCP<const Thyra::VectorSpaceBase<Scalar> >
315 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
316 "Error, initializeDIRKModel must be called first!\n" 318 return daeModel_->get_f_space();
322 template<
class Scalar>
323 RCP<Thyra::LinearOpBase<Scalar> >
326 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
327 "Error, initializeDIRKModel must be called first!\n" 329 return daeModel_->create_W_op();
333 template<
class Scalar>
334 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
337 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
338 "Error, initializeDIRKModel must be called first!\n" 340 return daeModel_->get_W_factory();
344 template<
class Scalar>
345 Thyra::ModelEvaluatorBase::InArgs<Scalar>
348 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
349 "Error, initializeDIRKModel must be called first!\n" 351 return nominalValues_;
355 template<
class Scalar>
356 Thyra::ModelEvaluatorBase::InArgs<Scalar>
359 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
360 "Error, initializeDIRKModel must be called first!\n" 369 template<
class Scalar>
370 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
373 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
374 "Error, initializeDIRKModel must be called first!\n" 380 template<
class Scalar>
382 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs_stage,
383 const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs_stage
387 typedef ScalarTraits<Scalar> ST;
388 typedef Thyra::ModelEvaluatorBase MEB;
390 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
391 "Error! initializeDIRKModel must be called before evalModel\n" 394 TEUCHOS_TEST_FOR_EXCEPTION( !setTimeStepPointCalled_, std::logic_error,
395 "Error! setTimeStepPoint must be called before evalModel" 398 TEUCHOS_TEST_FOR_EXCEPTION( currentStage_ == -1, std::logic_error,
399 "Error! setCurrentStage must be called before evalModel" 402 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
403 "Rythmos::DiagonalImplicitRKModelEvaluator",inArgs_stage,outArgs_stage,daeModel_
410 const RCP<const Thyra::VectorBase<Scalar> > x_in = inArgs_stage.get_x();
411 const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs_stage.get_f();
412 const RCP<Thyra::LinearOpBase<Scalar> > W_op_out = outArgs_stage.get_W_op();
418 MEB::InArgs<Scalar> daeInArgs = daeModel_->createInArgs();
419 MEB::OutArgs<Scalar> daeOutArgs = daeModel_->createOutArgs();
420 const RCP<Thyra::VectorBase<Scalar> > x_i = createMember(daeModel_->get_x_space());
421 daeInArgs.setArgs(basePoint_);
424 V_V(stage_derivatives_->getNonconstVectorBlock(currentStage_).ptr(),*x_in);
425 assembleIRKState( currentStage_, dirkButcherTableau_->A(), delta_t_, *x_old_, *stage_derivatives_, outArg(*x_i) );
426 daeInArgs.set_x( x_i );
427 daeInArgs.set_x_dot( x_in );
428 daeInArgs.set_t( t_old_ + dirkButcherTableau_->c()(currentStage_) * delta_t_ );
429 daeInArgs.set_alpha(ST::one());
430 daeInArgs.set_beta( delta_t_ * dirkButcherTableau_->A()(currentStage_,currentStage_) );
434 daeOutArgs.set_f( f_out );
435 if (!is_null(W_op_out))
436 daeOutArgs.set_W_op(W_op_out);
439 daeModel_->evalModel( daeInArgs, daeOutArgs );
440 daeOutArgs.set_f(Teuchos::null);
441 daeOutArgs.set_W_op(Teuchos::null);
443 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
451 #endif // RYTHMOS_DIAGONALIMPLICITRK_MODEL_EVALUATOR_HPP Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
void initializeDIRKModel(const RCP< const Thyra::ModelEvaluator< Scalar > > &daeModel, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &basePoint, const RCP< Thyra::LinearOpWithSolveFactoryBase< Scalar > > &dirk_W_factory, const RCP< const RKButcherTableauBase< Scalar > > &irkButcherTableau)
void setTimeStepPoint(const RCP< const Thyra::VectorBase< Scalar > > &x_old, const Scalar &t_old, const Scalar &delta_t)
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void setStageSolution(int stageNumber, const Thyra::VectorBase< Scalar > &stage_solution)
void setCurrentStage(int currentStage)
DiagonalImplicitRKModelEvaluator()
RCP< DiagonalImplicitRKModelEvaluator< Scalar > > diagonalImplicitRKModelEvaluator(const RCP< const Thyra::ModelEvaluator< Scalar > > &daeModel, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &basePoint, const RCP< Thyra::LinearOpWithSolveFactoryBase< Scalar > > &dirk_W_factory, const RCP< const RKButcherTableauBase< Scalar > > &irkButcherTableau)
Nonmember constructor function.
RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const