Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_DiagonalImplicitRKModelEvaluator.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
30#ifndef RYTHMOS_DIAGONALIMPLICITRK_MODEL_EVALUATOR_HPP
31#define RYTHMOS_DIAGONALIMPLICITRK_MODEL_EVALUATOR_HPP
32
33
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"
47
48
49namespace Rythmos {
50
51
53template<class Scalar>
54class DiagonalImplicitRKModelEvaluator : virtual public Thyra::StateFuncModelEvaluatorBase<Scalar>
55{
56public:
57
60
63
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
70 );
71
74 const RCP<const Thyra::VectorBase<Scalar> > &x_old,
75 const Scalar &t_old,
76 const Scalar &delta_t
77 );
78
81 int stageNumber,
82 const Thyra::VectorBase<Scalar>& stage_solution
83 );
84
86 void setCurrentStage(
87 int currentStage
88 );
89
91
94
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;
104 Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
106 Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
107
109
110private:
111
114
116 Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
118 void evalModelImpl(
119 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
120 const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
121 ) const;
122
124
125
126private:
127
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_;
132
133 Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs_;
134 Thyra::ModelEvaluatorBase::OutArgs<Scalar> outArgs_;
135 Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
136
137 RCP<const Thyra::ProductVectorSpaceBase<Scalar> > stage_space_;
138 RCP<Thyra::ProductVectorBase<Scalar> > stage_derivatives_;
139
140 bool setTimeStepPointCalled_;
141 RCP<const Thyra::VectorBase<Scalar> > x_old_;
142 Scalar t_old_;
143 Scalar delta_t_;
144 int currentStage_;
145
146 bool isInitialized_;
147
148};
149
150
155template<class Scalar>
156RCP<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
162 )
163{
164 RCP<DiagonalImplicitRKModelEvaluator<Scalar> >
165 dirkModel = rcp(new DiagonalImplicitRKModelEvaluator<Scalar>());
166 dirkModel->initializeDIRKModel(daeModel,basePoint,dirk_W_factory,irkButcherTableau);
167 return dirkModel;
168}
169
170
171// ///////////////////////
172// Definition
173
174
175// Constructors/initializers/accessors
176
177
178template<class Scalar>
180{
181 setTimeStepPointCalled_ = false;
182 isInitialized_ = false;
183 currentStage_ = -1;
184}
185
186
187// Overridden from ImplicitRKModelEvaluatorBase
188
189
190template<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
196 )
197{
198 typedef ScalarTraits<Scalar> ST;
199 // ToDo: Assert input arguments!
200 // How do I verify the basePoint is an authentic InArgs from daeModel?
201 TEUCHOS_TEST_FOR_EXCEPTION(
202 is_null(basePoint.get_x()),
203 std::logic_error,
204 "Error! The basepoint x vector is null!"
205 );
206 TEUCHOS_TEST_FOR_EXCEPTION(
207 is_null(daeModel),
208 std::logic_error,
209 "Error! The model evaluator pointer is null!"
210 );
211 TEUCHOS_TEST_FOR_EXCEPTION(
212 !daeModel->get_x_space()->isCompatible(*(basePoint.get_x()->space())),
213 std::logic_error,
214 "Error! The basepoint input arguments are incompatible with the model evaluator vector space!"
215 );
216 //TEUCHOS_TEST_FOR_EXCEPT(is_null(dirk_W_factory));
217
218 daeModel_ = daeModel;
219 basePoint_ = basePoint;
220 dirk_W_factory_ = dirk_W_factory;
221 dirkButcherTableau_ = irkButcherTableau;
222
223 const int numStages = dirkButcherTableau_->numStages();
224
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());
230
231 // Set up prototypical InArgs
232 {
233 typedef Thyra::ModelEvaluatorBase MEB;
234 MEB::InArgsSetup<Scalar> inArgs;
235 inArgs.setModelEvalDescription(this->description());
236 inArgs.setSupports(MEB::IN_ARG_x);
237 inArgs.setSupports(MEB::IN_ARG_step_size);
238 inArgs.setSupports(MEB::IN_ARG_stage_number);
239 inArgs_ = inArgs;
240 }
241 // Set up prototypical OutArgs
242 {
243 typedef Thyra::ModelEvaluatorBase MEB;
244 MEB::OutArgsSetup<Scalar> outArgs;
245 outArgs.setModelEvalDescription(this->description());
246 outArgs.setSupports(MEB::OUT_ARG_f);
247 outArgs.setSupports(MEB::OUT_ARG_W_op);
248 outArgs_ = outArgs;
249 }
250 // Set up nominal values
251 nominalValues_ = inArgs_;
252
253 isInitialized_ = true;
254}
255
256
257template<class Scalar>
259 const RCP<const Thyra::VectorBase<Scalar> > &x_old,
260 const Scalar &t_old,
261 const Scalar &delta_t
262 )
263{
264 typedef ScalarTraits<Scalar> ST;
265 TEUCHOS_TEST_FOR_EXCEPT( is_null(x_old) );
266 TEUCHOS_TEST_FOR_EXCEPT( delta_t <= ST::zero() );
267 // Verify x_old is compatible with the vector space in the DAE Model.
268 TEUCHOS_TEST_FOR_EXCEPTION(!daeModel_->get_x_space()->isCompatible(*(x_old->space())), std::logic_error,
269 "Error! The incoming VectorBase object is not compatible with the DAE model that was provided at initialization!");
270 x_old_ = x_old;
271 t_old_ = t_old;
272 delta_t_ = delta_t;
273 setTimeStepPointCalled_ = true;
274}
275
276template<class Scalar>
278 int stageNumber,
279 const Thyra::VectorBase<Scalar>& stage_solution
280 )
281{
282 TEUCHOS_TEST_FOR_EXCEPT(stageNumber < 0);
283 TEUCHOS_TEST_FOR_EXCEPT(stageNumber >= dirkButcherTableau_->numStages());
284 Thyra::V_V(stage_derivatives_->getNonconstVectorBlock(stageNumber).ptr(),stage_solution);
285}
286
287template<class Scalar>
289 int currentStage
290 )
291{
292 TEUCHOS_TEST_FOR_EXCEPT(currentStage < 0);
293 TEUCHOS_TEST_FOR_EXCEPT(currentStage >= dirkButcherTableau_->numStages());
294 currentStage_ = currentStage;
295}
296
297
298
299// Overridden from ModelEvaluator
300
301
302template<class Scalar>
303RCP<const Thyra::VectorSpaceBase<Scalar> >
305{
306 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
307 "Error, initializeDIRKModel must be called first!\n"
308 );
309 return daeModel_->get_x_space();
310}
311
312
313template<class Scalar>
314RCP<const Thyra::VectorSpaceBase<Scalar> >
316{
317 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
318 "Error, initializeDIRKModel must be called first!\n"
319 );
320 return daeModel_->get_f_space();
321}
322
323
324template<class Scalar>
325RCP<Thyra::LinearOpBase<Scalar> >
327{
328 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
329 "Error, initializeDIRKModel must be called first!\n"
330 );
331 return daeModel_->create_W_op();
332}
333
334
335template<class Scalar>
336RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
338{
339 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
340 "Error, initializeDIRKModel must be called first!\n"
341 );
342 return daeModel_->get_W_factory();
343}
344
345
346template<class Scalar>
347Thyra::ModelEvaluatorBase::InArgs<Scalar>
349{
350 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
351 "Error, initializeDIRKModel must be called first!\n"
352 );
353 return nominalValues_;
354}
355
356
357template<class Scalar>
358Thyra::ModelEvaluatorBase::InArgs<Scalar>
360{
361 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
362 "Error, initializeDIRKModel must be called first!\n"
363 );
364
365 return inArgs_;
366}
367
368
369// Private functions overridden from ModelEvaluatorDefaultBase
370
371
372template<class Scalar>
373Thyra::ModelEvaluatorBase::OutArgs<Scalar>
375{
376 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
377 "Error, initializeDIRKModel must be called first!\n"
378 );
379 return outArgs_;
380}
381
382
383template<class Scalar>
384void DiagonalImplicitRKModelEvaluator<Scalar>::evalModelImpl(
385 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs_stage,
386 const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs_stage
387 ) const
388{
389
390
391 typedef ScalarTraits<Scalar> ST;
392 typedef Thyra::ModelEvaluatorBase MEB;
393 typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag TScalarMag;
394
395 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
396 "Error! initializeDIRKModel must be called before evalModel\n"
397 );
398
399 TEUCHOS_TEST_FOR_EXCEPTION( !setTimeStepPointCalled_, std::logic_error,
400 "Error! setTimeStepPoint must be called before evalModel"
401 );
402
403 TEUCHOS_TEST_FOR_EXCEPTION( currentStage_ == -1, std::logic_error,
404 "Error! setCurrentStage must be called before evalModel"
405 );
406
407 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
408 "Rythmos::DiagonalImplicitRKModelEvaluator",inArgs_stage,outArgs_stage,daeModel_
409 );
410
411 //
412 // A) Unwrap the inArgs and outArgs
413 //
414
415 const RCP<const Thyra::VectorBase<Scalar> > x_in = inArgs_stage.get_x();
416 const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs_stage.get_f();
417 const RCP<Thyra::LinearOpBase<Scalar> > W_op_out = outArgs_stage.get_W_op();
418
419 //
420 // B) Assemble f_out and W_op_out for given stage
421 //
422
423 MEB::InArgs<Scalar> daeInArgs = daeModel_->createInArgs();
424 MEB::OutArgs<Scalar> daeOutArgs = daeModel_->createOutArgs();
425 const RCP<Thyra::VectorBase<Scalar> > x_i = createMember(daeModel_->get_x_space());
426 daeInArgs.setArgs(basePoint_);
427
428 // B.1) Setup the DAE's inArgs for stage f(currentStage_) ...
429 V_V(stage_derivatives_->getNonconstVectorBlock(currentStage_).ptr(),*x_in);
430 assembleIRKState( currentStage_, dirkButcherTableau_->A(), delta_t_, *x_old_, *stage_derivatives_, outArg(*x_i) );
431 daeInArgs.set_x( x_i );
432 daeInArgs.set_x_dot( x_in );
433 daeInArgs.set_t( t_old_ + dirkButcherTableau_->c()(currentStage_) * delta_t_ );
434 daeInArgs.set_alpha(ST::one());
435 daeInArgs.set_beta( delta_t_ * dirkButcherTableau_->A()(currentStage_,currentStage_) );
436
437 if (daeInArgs.supports(MEB::IN_ARG_step_size)) {
438 TScalarMag scaled_dt;
439
440 if (currentStage_ == 0 && dirkButcherTableau_->A()(currentStage_,currentStage_) == 0) {
441 scaled_dt = Scalar( delta_t_ /dirkButcherTableau_->numStages() );
442 } else {
443 scaled_dt = dirkButcherTableau_->c()(currentStage_) * delta_t_;
444 }
445
446 daeInArgs.set_step_size( scaled_dt );
447 }
448
449 if (daeInArgs.supports(MEB::IN_ARG_stage_number)) {
450 daeInArgs.set_stage_number( dirkButcherTableau_->c()(currentStage_) );
451 }
452
453 // B.2) Setup the DAE's outArgs for stage f(i) ...
454 if (!is_null(f_out))
455 daeOutArgs.set_f( f_out );
456 if (!is_null(W_op_out))
457 daeOutArgs.set_W_op(W_op_out);
458
459 // B.3) Compute f_out(i) and/or W_op_out ...
460 daeModel_->evalModel( daeInArgs, daeOutArgs );
461 daeOutArgs.set_f(Teuchos::null);
462 daeOutArgs.set_W_op(Teuchos::null);
463
464 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
465
466}
467
468
469} // namespace Rythmos
470
471
472#endif // RYTHMOS_DIAGONALIMPLICITRK_MODEL_EVALUATOR_HPP
RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
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< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
void setTimeStepPoint(const RCP< const Thyra::VectorBase< Scalar > > &x_old, const Scalar &t_old, const Scalar &delta_t)
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
void setStageSolution(int stageNumber, const Thyra::VectorBase< Scalar > &stage_solution)
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() 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)