Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_ImplicitRKModelEvaluator.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_IMPLICITRK_MODEL_EVALUATOR_HPP
31#define RYTHMOS_IMPLICITRK_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 ImplicitRKModelEvaluator : 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> >& irk_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
80
83
85 RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
87 RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
89 RCP<Thyra::LinearOpBase<Scalar> > create_W_op() const;
91 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > get_W_factory() const;
93 Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
95 Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
96
98
99private:
100
103
105 Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
107 void evalModelImpl(
108 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
109 const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
110 ) const;
111
113
114
115private:
116
117 RCP<const Thyra::ModelEvaluator<Scalar> > daeModel_;
118 Thyra::ModelEvaluatorBase::InArgs<Scalar> basePoint_;
119 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > irk_W_factory_;
120 RCP<const RKButcherTableauBase<Scalar> > irkButcherTableau_;
121
122 Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs_;
123 Thyra::ModelEvaluatorBase::OutArgs<Scalar> outArgs_;
124 Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
125
126 RCP<const Thyra::ProductVectorSpaceBase<Scalar> > x_bar_space_;
127 RCP<const Thyra::ProductVectorSpaceBase<Scalar> > f_bar_space_;
128 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_bar_factory_;
129
130 bool setTimeStepPointCalled_;
131 RCP<const Thyra::VectorBase<Scalar> > x_old_;
132 Scalar t_old_;
133 Scalar delta_t_;
134
135 bool isInitialized_;
136
137};
138
139
144template<class Scalar>
145RCP<ImplicitRKModelEvaluator<Scalar> >
147 const RCP<const Thyra::ModelEvaluator<Scalar> >& daeModel,
148 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& basePoint,
149 const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >& irk_W_factory,
150 const RCP<const RKButcherTableauBase<Scalar> >& irkButcherTableau
151 )
152{
153 RCP<ImplicitRKModelEvaluator<Scalar> >
154 irkModel = rcp(new ImplicitRKModelEvaluator<Scalar>());
155 irkModel->initializeIRKModel(daeModel,basePoint,irk_W_factory,irkButcherTableau);
156 return irkModel;
157}
158
159
160// ///////////////////////
161// Definition
162
163
164// Constructors/initializers/accessors
165
166
167template<class Scalar>
169{
170 setTimeStepPointCalled_ = false;
171 isInitialized_ = false;
172}
173
174
175// Overridden from ImplicitRKModelEvaluatorBase
176
177
178template<class Scalar>
180 const RCP<const Thyra::ModelEvaluator<Scalar> >& daeModel,
181 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& basePoint,
182 const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >& irk_W_factory,
183 const RCP<const RKButcherTableauBase<Scalar> >& irkButcherTableau
184 )
185{
186 // ToDo: Assert input arguments!
187 // How do I verify the basePoint is an authentic InArgs from daeModel?
188 TEUCHOS_TEST_FOR_EXCEPTION(
189 is_null(basePoint.get_x()),
190 std::logic_error,
191 "Error! The basepoint x vector is null!"
192 );
193 TEUCHOS_TEST_FOR_EXCEPTION(
194 is_null(daeModel),
195 std::logic_error,
196 "Error! The model evaluator pointer is null!"
197 );
198 TEUCHOS_TEST_FOR_EXCEPTION(
199 !daeModel->get_x_space()->isCompatible(*(basePoint.get_x()->space())),
200 std::logic_error,
201 "Error! The basepoint input arguments are incompatible with the model evaluator vector space!"
202 );
203 TEUCHOS_TEST_FOR_EXCEPT(is_null(irk_W_factory));
204
205 daeModel_ = daeModel;
206 basePoint_ = basePoint;
207 irk_W_factory_ = irk_W_factory;
208 irkButcherTableau_ = irkButcherTableau;
209
210 const int numStages = irkButcherTableau_->numStages();
211
212 x_bar_space_ = productVectorSpace(daeModel_->get_x_space(),numStages);
213 f_bar_space_ = productVectorSpace(daeModel_->get_f_space(),numStages);
214
215 // HACK! Remove the preconditioner factory for now!
216 if (irk_W_factory_->acceptsPreconditionerFactory())
217 irk_W_factory_->unsetPreconditionerFactory();
218
219 // ToDo: create the block diagonal preconditioner factory and set this on
220 // irk_W_factory_!
221
222 // Set up prototypical InArgs
223 {
224 typedef Thyra::ModelEvaluatorBase MEB;
225 MEB::InArgsSetup<Scalar> inArgs;
226 inArgs.setModelEvalDescription(this->description());
227 inArgs.setSupports(MEB::IN_ARG_x);
228 inArgs_ = inArgs;
229 }
230 // Set up prototypical OutArgs
231 {
232 typedef Thyra::ModelEvaluatorBase MEB;
233 MEB::OutArgsSetup<Scalar> outArgs;
234 outArgs.setModelEvalDescription(this->description());
235 outArgs.setSupports(MEB::OUT_ARG_f);
236 outArgs.setSupports(MEB::OUT_ARG_W_op);
237 outArgs_ = outArgs;
238 }
239 // Set up nominal values
240 nominalValues_ = inArgs_;
241
242 isInitialized_ = true;
243}
244
245
246template<class Scalar>
248 const RCP<const Thyra::VectorBase<Scalar> > &x_old,
249 const Scalar &t_old,
250 const Scalar &delta_t
251 )
252{
253 typedef ScalarTraits<Scalar> ST;
254 TEUCHOS_TEST_FOR_EXCEPT( is_null(x_old) );
255 TEUCHOS_TEST_FOR_EXCEPT( delta_t <= ST::zero() );
256 // Verify x_old is compatible with the vector space in the DAE Model.
257 TEUCHOS_TEST_FOR_EXCEPTION(!daeModel_->get_x_space()->isCompatible(*(x_old->space())), std::logic_error,
258 "Error! The incoming VectorBase object is not compatible with the DAE model that was provided at initialization!");
259 x_old_ = x_old;
260 t_old_ = t_old;
261 delta_t_ = delta_t;
262 setTimeStepPointCalled_ = true;
263}
264
265
266// Overridden from ModelEvaluator
267
268
269template<class Scalar>
270RCP<const Thyra::VectorSpaceBase<Scalar> >
272{
273 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
274 "Error, initializeIRKModel must be called first!\n"
275 );
276 return x_bar_space_;
277}
278
279
280template<class Scalar>
281RCP<const Thyra::VectorSpaceBase<Scalar> >
283{
284 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
285 "Error, initializeIRKModel must be called first!\n"
286 );
287 return f_bar_space_;
288}
289
290
291template<class Scalar>
292RCP<Thyra::LinearOpBase<Scalar> >
294{
295 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
296 "Error, initializeIRKModel must be called first!\n"
297 );
298 // Create the block structure for W_op_bar right away!
299 const int numStages = irkButcherTableau_->numStages();
300 RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar> >
301 W_op_bar = Thyra::defaultBlockedLinearOp<Scalar>();
302 W_op_bar->beginBlockFill( f_bar_space_, x_bar_space_ );
303 for ( int i = 0; i < numStages; ++i )
304 for ( int j = 0; j < numStages; ++j )
305 W_op_bar->setNonconstBlock( i, j, daeModel_->create_W_op() );
306 W_op_bar->endBlockFill();
307 return W_op_bar;
308}
309
310
311template<class Scalar>
312RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
314{
315 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
316 "Error, initializeIRKModel must be called first!\n"
317 );
318 return irk_W_factory_;
319}
320
321
322template<class Scalar>
323Thyra::ModelEvaluatorBase::InArgs<Scalar>
325{
326 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
327 "Error, initializeIRKModel must be called first!\n"
328 );
329 return nominalValues_;
330}
331
332
333template<class Scalar>
334Thyra::ModelEvaluatorBase::InArgs<Scalar>
336{
337 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
338 "Error, initializeIRKModel must be called first!\n"
339 );
340 return inArgs_;
341}
342
343
344// Private functions overridden from ModelEvaluatorDefaultBase
345
346
347template<class Scalar>
348Thyra::ModelEvaluatorBase::OutArgs<Scalar>
350{
351 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
352 "Error, initializeIRKModel must be called first!\n"
353 );
354 return outArgs_;
355}
356
357
358template<class Scalar>
359void ImplicitRKModelEvaluator<Scalar>::evalModelImpl(
360 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs_bar,
361 const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs_bar
362 ) const
363{
364
365 using Teuchos::rcp_dynamic_cast;
366 typedef ScalarTraits<Scalar> ST;
367 typedef Thyra::ModelEvaluatorBase MEB;
368 typedef Thyra::VectorBase<Scalar> VB;
369 typedef Thyra::ProductVectorBase<Scalar> PVB;
370 typedef Thyra::BlockedLinearOpBase<Scalar> BLWB;
371
372 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
373 "Error! initializeIRKModel must be called before evalModel\n"
374 );
375
376 TEUCHOS_TEST_FOR_EXCEPTION( !setTimeStepPointCalled_, std::logic_error,
377 "Error! setTimeStepPoint must be called before evalModel"
378 );
379
380 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
381 "Rythmos::ImplicitRKModelEvaluator",inArgs_bar,outArgs_bar,daeModel_
382 );
383
384 //
385 // A) Unwrap the inArgs and outArgs to get at product vectors and block op
386 //
387
388 const RCP<const PVB> x_bar = rcp_dynamic_cast<const PVB>(inArgs_bar.get_x(), true);
389 const RCP<PVB> f_bar = rcp_dynamic_cast<PVB>(outArgs_bar.get_f(), true);
390 const RCP<BLWB> W_op_bar = rcp_dynamic_cast<BLWB>(outArgs_bar.get_W_op(), true);
391
392 //
393 // B) Assemble f_bar and W_op_bar by looping over stages
394 //
395
396 MEB::InArgs<Scalar> daeInArgs = daeModel_->createInArgs();
397 MEB::OutArgs<Scalar> daeOutArgs = daeModel_->createOutArgs();
398 const RCP<VB> x_i = createMember(daeModel_->get_x_space());
399 daeInArgs.setArgs(basePoint_);
400
401 const int numStages = irkButcherTableau_->numStages();
402
403 for ( int i = 0; i < numStages; ++i ) {
404
405 // B.1) Setup the DAE's inArgs for stage f(i) ...
406 assembleIRKState( i, irkButcherTableau_->A(), delta_t_, *x_old_, *x_bar, outArg(*x_i) );
407 daeInArgs.set_x( x_i );
408 daeInArgs.set_x_dot( x_bar->getVectorBlock(i) );
409 daeInArgs.set_t( t_old_ + irkButcherTableau_->c()(i) * delta_t_ );
410 Scalar alpha = ST::zero();
411 if (i == 0) {
412 alpha = ST::one();
413 } else {
414 alpha = ST::zero();
415 }
416 Scalar beta = delta_t_ * irkButcherTableau_->A()(i,0);
417 daeInArgs.set_alpha( alpha );
418 daeInArgs.set_beta( beta );
419
420 // B.2) Setup the DAE's outArgs for stage f(i) ...
421 if (!is_null(f_bar))
422 daeOutArgs.set_f( f_bar->getNonconstVectorBlock(i) );
423 if (!is_null(W_op_bar)) {
424 daeOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,0));
425 }
426
427 // B.3) Compute f_bar(i) and/or W_op_bar(i,0) ...
428 daeModel_->evalModel( daeInArgs, daeOutArgs );
429 daeOutArgs.set_f(Teuchos::null);
430 daeOutArgs.set_W_op(Teuchos::null);
431
432 // B.4) Evaluate the rest of the W_op_bar(i,j=1...numStages-1) ...
433 if (!is_null(W_op_bar)) {
434 for ( int j = 1; j < numStages; ++j ) {
435 alpha = ST::zero();
436 if (i == j) {
437 alpha = ST::one();
438 } else {
439 alpha = ST::zero();
440 }
441 beta = delta_t_ * irkButcherTableau_->A()(i,j);
442 daeInArgs.set_alpha( alpha );
443 daeInArgs.set_beta( beta );
444 daeOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,j));
445 daeModel_->evalModel( daeInArgs, daeOutArgs );
446 daeOutArgs.set_W_op(Teuchos::null);
447 }
448 }
449
450 }
451
452 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
453
454}
455
456
457} // namespace Rythmos
458
459
460#endif // RYTHMOS_IMPLICITRK_MODEL_EVALUATOR_HPP
RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
void initializeIRKModel(const RCP< const Thyra::ModelEvaluator< Scalar > > &daeModel, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &basePoint, const RCP< Thyra::LinearOpWithSolveFactoryBase< Scalar > > &irk_W_factory, const RCP< const RKButcherTableauBase< Scalar > > &irkButcherTableau)
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
RCP< ImplicitRKModelEvaluator< Scalar > > implicitRKModelEvaluator(const RCP< const Thyra::ModelEvaluator< Scalar > > &daeModel, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &basePoint, const RCP< Thyra::LinearOpWithSolveFactoryBase< Scalar > > &irk_W_factory, const RCP< const RKButcherTableauBase< Scalar > > &irkButcherTableau)
Nonmember constructor function.
void setTimeStepPoint(const RCP< const Thyra::VectorBase< Scalar > > &x_old, const Scalar &t_old, const Scalar &delta_t)