Rythmos - Transient Integration for Differential Equations  Version of the Day
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 
49 namespace Rythmos {
50 
51 
53 template<class Scalar>
54 class DiagonalImplicitRKModelEvaluator : virtual public Thyra::StateFuncModelEvaluatorBase<Scalar>
55 {
56 public:
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 
73  void setTimeStepPoint(
74  const RCP<const Thyra::VectorBase<Scalar> > &x_old,
75  const Scalar &t_old,
76  const Scalar &delta_t
77  );
78 
80  void setStageSolution(
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 
110 private:
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 
126 private:
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 
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
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 
178 template<class Scalar>
180 {
181  setTimeStepPointCalled_ = false;
182  isInitialized_ = false;
183  currentStage_ = -1;
184 }
185 
186 
187 // Overridden from ImplicitRKModelEvaluatorBase
188 
189 
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
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_ = inArgs;
238  }
239  // Set up prototypical OutArgs
240  {
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);
246  outArgs_ = outArgs;
247  }
248  // Set up nominal values
249  nominalValues_ = inArgs_;
250 
251  isInitialized_ = true;
252 }
253 
254 
255 template<class Scalar>
257  const RCP<const Thyra::VectorBase<Scalar> > &x_old,
258  const Scalar &t_old,
259  const Scalar &delta_t
260  )
261 {
262  typedef ScalarTraits<Scalar> ST;
263  TEUCHOS_TEST_FOR_EXCEPT( is_null(x_old) );
264  TEUCHOS_TEST_FOR_EXCEPT( delta_t <= ST::zero() );
265  // Verify x_old is compatible with the vector space in the DAE Model.
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!");
268  x_old_ = x_old;
269  t_old_ = t_old;
270  delta_t_ = delta_t;
271  setTimeStepPointCalled_ = true;
272 }
273 
274 template<class Scalar>
276  int stageNumber,
277  const Thyra::VectorBase<Scalar>& stage_solution
278  )
279 {
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);
283 }
284 
285 template<class Scalar>
287  int currentStage
288  )
289 {
290  TEUCHOS_TEST_FOR_EXCEPT(currentStage < 0);
291  TEUCHOS_TEST_FOR_EXCEPT(currentStage >= dirkButcherTableau_->numStages());
292  currentStage_ = currentStage;
293 }
294 
295 
296 
297 // Overridden from ModelEvaluator
298 
299 
300 template<class Scalar>
301 RCP<const Thyra::VectorSpaceBase<Scalar> >
303 {
304  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
305  "Error, initializeDIRKModel must be called first!\n"
306  );
307  return daeModel_->get_x_space();
308 }
309 
310 
311 template<class Scalar>
312 RCP<const Thyra::VectorSpaceBase<Scalar> >
314 {
315  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
316  "Error, initializeDIRKModel must be called first!\n"
317  );
318  return daeModel_->get_f_space();
319 }
320 
321 
322 template<class Scalar>
323 RCP<Thyra::LinearOpBase<Scalar> >
325 {
326  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
327  "Error, initializeDIRKModel must be called first!\n"
328  );
329  return daeModel_->create_W_op();
330 }
331 
332 
333 template<class Scalar>
334 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
336 {
337  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
338  "Error, initializeDIRKModel must be called first!\n"
339  );
340  return daeModel_->get_W_factory();
341 }
342 
343 
344 template<class Scalar>
345 Thyra::ModelEvaluatorBase::InArgs<Scalar>
347 {
348  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
349  "Error, initializeDIRKModel must be called first!\n"
350  );
351  return nominalValues_;
352 }
353 
354 
355 template<class Scalar>
356 Thyra::ModelEvaluatorBase::InArgs<Scalar>
358 {
359  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
360  "Error, initializeDIRKModel must be called first!\n"
361  );
362  return inArgs_;
363 }
364 
365 
366 // Private functions overridden from ModelEvaluatorDefaultBase
367 
368 
369 template<class Scalar>
370 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
372 {
373  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
374  "Error, initializeDIRKModel must be called first!\n"
375  );
376  return outArgs_;
377 }
378 
379 
380 template<class Scalar>
382  const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs_stage,
383  const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs_stage
384  ) const
385 {
386 
387  typedef ScalarTraits<Scalar> ST;
388  typedef Thyra::ModelEvaluatorBase MEB;
389 
390  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
391  "Error! initializeDIRKModel must be called before evalModel\n"
392  );
393 
394  TEUCHOS_TEST_FOR_EXCEPTION( !setTimeStepPointCalled_, std::logic_error,
395  "Error! setTimeStepPoint must be called before evalModel"
396  );
397 
398  TEUCHOS_TEST_FOR_EXCEPTION( currentStage_ == -1, std::logic_error,
399  "Error! setCurrentStage must be called before evalModel"
400  );
401 
402  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
403  "Rythmos::DiagonalImplicitRKModelEvaluator",inArgs_stage,outArgs_stage,daeModel_
404  );
405 
406  //
407  // A) Unwrap the inArgs and outArgs
408  //
409 
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();
413 
414  //
415  // B) Assemble f_out and W_op_out for given stage
416  //
417 
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_);
422 
423  // B.1) Setup the DAE's inArgs for stage f(currentStage_) ...
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_) );
431 
432  // B.2) Setup the DAE's outArgs for stage f(i) ...
433  if (!is_null(f_out))
434  daeOutArgs.set_f( f_out );
435  if (!is_null(W_op_out))
436  daeOutArgs.set_W_op(W_op_out);
437 
438  // B.3) Compute f_out(i) and/or W_op_out ...
439  daeModel_->evalModel( daeInArgs, daeOutArgs );
440  daeOutArgs.set_f(Teuchos::null);
441  daeOutArgs.set_W_op(Teuchos::null);
442 
443  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
444 
445 }
446 
447 
448 } // namespace Rythmos
449 
450 
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)
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