Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
Rythmos_TimeDiscretizedBackwardEulerModelEvaluator.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_TIME_DISCRETIZED_BACKWARD_EULER_MODEL_EVALUATOR_HPP
31 #define RYTHMOS_TIME_DISCRETIZED_BACKWARD_EULER_MODEL_EVALUATOR_HPP
32 
33 
34 #include "Rythmos_Types.hpp"
35 #include "Thyra_StateFuncModelEvaluatorBase.hpp"
36 #include "Thyra_ProductVectorBase.hpp"
37 #include "Thyra_DefaultProductVectorSpace.hpp"
38 #include "Thyra_DefaultBlockedLinearOp.hpp"
39 #include "Thyra_DefaultBlockedTriangularLinearOpWithSolveFactory.hpp" // Default implementation
40 
41 
42 namespace Rythmos {
43 
44 
49 template<class Scalar>
51  : virtual public Thyra::StateFuncModelEvaluatorBase<Scalar>
52 {
53 public:
54 
57 
60 
62 
65 
67  void initialize(
68  const RCP<const Thyra::ModelEvaluator<Scalar> > &daeModel,
69  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initCond,
70  const Scalar finalTime,
71  const int numTimeSteps,
72  const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory = Teuchos::null
73  );
74 
76 
79 
81  RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
83  RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
85  RCP<Thyra::LinearOpBase<Scalar> > create_W_op() const;
87  RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > get_W_factory() const;
89  Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
91  Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
92 
94 
95 private:
96 
99 
101  Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
103  void evalModelImpl(
104  const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
105  const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
106  ) const;
107 
109 
110 private:
111 
112  RCP<const Thyra::ModelEvaluator<Scalar> > daeModel_;
113  Thyra::ModelEvaluatorBase::InArgs<Scalar> initCond_;
114  Scalar finalTime_;
115  int numTimeSteps_;
116 
117  Scalar initTime_;
118  Scalar delta_t_;
119 
120  RCP<const Thyra::ProductVectorSpaceBase<Scalar> > x_bar_space_;
121  RCP<const Thyra::ProductVectorSpaceBase<Scalar> > f_bar_space_;
122  RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_bar_factory_;
123 
124 };
125 
126 
131 template<class Scalar>
132 RCP<TimeDiscretizedBackwardEulerModelEvaluator<Scalar> >
134  const RCP<const Thyra::ModelEvaluator<Scalar> > &daeModel,
135  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initCond,
136  const Scalar finalTime,
137  const int numTimeSteps,
138  const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory
139  )
140 {
141  RCP<TimeDiscretizedBackwardEulerModelEvaluator<Scalar> >
143  model->initialize(daeModel,initCond,finalTime,numTimeSteps,W_bar_factory);
144  return model;
145 }
146 
147 
148 // ///////////////////////
149 // Definition
150 
151 
152 // Constructors/initializers/accessors
153 
154 
155 template<class Scalar>
157  :finalTime_(-1.0),
158  numTimeSteps_(-1),
159  initTime_(0.0),
160  delta_t_(-1.0) // Flag for uninitialized!
161 {}
162 
163 
164 // Overridden from TimeDiscretizedBackwardEulerModelEvaluatorBase
165 
166 
167 template<class Scalar>
169  const RCP<const Thyra::ModelEvaluator<Scalar> > &daeModel,
170  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initCond,
171  const Scalar finalTime,
172  const int numTimeSteps,
173  const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory
174  )
175 {
176 
177  TEUCHOS_TEST_FOR_EXCEPT(is_null(daeModel));
178  TEUCHOS_TEST_FOR_EXCEPT(is_null(initCond.get_x()));
179  TEUCHOS_TEST_FOR_EXCEPT(is_null(initCond.get_x_dot()));
180  TEUCHOS_TEST_FOR_EXCEPT(finalTime <= initCond.get_t());
181  TEUCHOS_TEST_FOR_EXCEPT(numTimeSteps <= 0);
182  // ToDo: Validate that daeModel is of the right form!
183 
184  daeModel_ = daeModel;
185  initCond_ = initCond;
186  finalTime_ = finalTime;
187  numTimeSteps_ = numTimeSteps;
188 
189  initTime_ = initCond.get_t();
190  delta_t_ = (finalTime_ - initTime_) / numTimeSteps_;
191 
192  x_bar_space_ = productVectorSpace(daeModel_->get_x_space(),numTimeSteps_);
193  f_bar_space_ = productVectorSpace(daeModel_->get_f_space(),numTimeSteps_);
194 
195  if (!is_null(W_bar_factory)) {
196  W_bar_factory_ = W_bar_factory;
197  }
198  else {
199  W_bar_factory_ =
200  Thyra::defaultBlockedTriangularLinearOpWithSolveFactory<Scalar>(
201  daeModel_->get_W_factory()
202  );
203  }
204 
205 }
206 
207 
208 // Public functions overridden from ModelEvaluator
209 
210 
211 template<class Scalar>
212 RCP<const Thyra::VectorSpaceBase<Scalar> >
214 {
215  return x_bar_space_;
216 }
217 
218 
219 template<class Scalar>
220 RCP<const Thyra::VectorSpaceBase<Scalar> >
222 {
223  return f_bar_space_;
224 }
225 
226 
227 template<class Scalar>
228 RCP<Thyra::LinearOpBase<Scalar> >
230 {
231  // Create the block structure for W_op_bar right away!
232  RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar> >
233  W_op_bar = Thyra::defaultBlockedLinearOp<Scalar>();
234  W_op_bar->beginBlockFill( f_bar_space_, x_bar_space_ );
235  for ( int k = 0; k < numTimeSteps_; ++k ) {
236  W_op_bar->setNonconstBlock( k, k, daeModel_->create_W_op() );
237  if (k > 0)
238  W_op_bar->setNonconstBlock( k, k-1, daeModel_->create_W_op() );
239  }
240  W_op_bar->endBlockFill();
241  return W_op_bar;
242 }
243 
244 
245 template<class Scalar>
246 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
248 {
249  return W_bar_factory_;
250 }
251 
252 
253 template<class Scalar>
254 Thyra::ModelEvaluatorBase::InArgs<Scalar>
256 {
257  typedef Thyra::ModelEvaluatorBase MEB;
258  TEUCHOS_TEST_FOR_EXCEPT(true);
259  return MEB::InArgs<Scalar>();
260 }
261 
262 
263 template<class Scalar>
264 Thyra::ModelEvaluatorBase::InArgs<Scalar>
266 {
267  typedef Thyra::ModelEvaluatorBase MEB;
268  MEB::InArgsSetup<Scalar> inArgs;
269  inArgs.setModelEvalDescription(this->description());
270  inArgs.setSupports(MEB::IN_ARG_x);
271  return inArgs;
272 }
273 
274 
275 // Private functions overridden from ModelEvaluatorDefaultBase
276 
277 
278 template<class Scalar>
279 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
281 {
282  typedef Thyra::ModelEvaluatorBase MEB;
283  MEB::OutArgs<Scalar> daeOutArgs = daeModel_->createOutArgs();
284  MEB::OutArgsSetup<Scalar> outArgs;
285  outArgs.setModelEvalDescription(this->description());
286  outArgs.setSupports(MEB::OUT_ARG_f);
287  outArgs.setSupports(MEB::OUT_ARG_W_op);
288  outArgs.set_W_properties(daeOutArgs.get_W_properties());
289  return outArgs;
290 }
291 
292 
293 template<class Scalar>
294 void TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::evalModelImpl(
295  const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs_bar,
296  const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs_bar
297  ) const
298 {
299 
300 
301  using Teuchos::rcp_dynamic_cast;
302  // typedef ScalarTraits<Scalar> ST; // unused
303  typedef Thyra::ModelEvaluatorBase MEB;
304  typedef Thyra::VectorBase<Scalar> VB;
305  typedef Thyra::ProductVectorBase<Scalar> PVB;
306  typedef Thyra::BlockedLinearOpBase<Scalar> BLWB;
307 
308 /*
309  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
310  "Rythmos::ImplicitRKModelEvaluator",inArgs_bar,outArgs_bar,daeModel_
311  );
312 */
313 
314  TEUCHOS_TEST_FOR_EXCEPTION( delta_t_ <= 0.0, std::logic_error,
315  "Error, you have not initialized this object correctly!" );
316 
317  //
318  // A) Unwrap the inArgs and outArgs to get at product vectors and block op
319  //
320 
321  const RCP<const PVB> x_bar = rcp_dynamic_cast<const PVB>(inArgs_bar.get_x(), true);
322  const RCP<PVB> f_bar = rcp_dynamic_cast<PVB>(outArgs_bar.get_f(), true);
323  RCP<BLWB> W_op_bar = rcp_dynamic_cast<BLWB>(outArgs_bar.get_W_op(), true);
324 
325  //
326  // B) Assemble f_bar and W_op_bar by looping over stages
327  //
328 
329  MEB::InArgs<Scalar> daeInArgs = daeModel_->createInArgs();
330  MEB::OutArgs<Scalar> daeOutArgs = daeModel_->createOutArgs();
331  const RCP<VB> x_dot_i = createMember(daeModel_->get_x_space());
332  daeInArgs.setArgs(initCond_);
333 
334  Scalar t_i = initTime_; // ToDo: Define t_init!
335 
336  const Scalar oneOverDeltaT = 1.0/delta_t_;
337 
338  for ( int i = 0; i < numTimeSteps_; ++i ) {
339 
340  // B.1) Setup the DAE's inArgs for time step eqn f(i) ...
341  const RCP<const Thyra::VectorBase<Scalar> >
342  x_i = x_bar->getVectorBlock(i),
343  x_im1 = ( i==0 ? initCond_.get_x() : x_bar->getVectorBlock(i-1) );
344  V_VmV( x_dot_i.ptr(), *x_i, *x_im1 ); // x_dot_i = 1/dt * ( x[i] - x[i-1] )
345  Vt_S( x_dot_i.ptr(), oneOverDeltaT ); // ...
346  daeInArgs.set_x_dot( x_dot_i );
347  daeInArgs.set_x( x_i );
348  daeInArgs.set_t( t_i );
349  daeInArgs.set_alpha( oneOverDeltaT );
350  daeInArgs.set_beta( 1.0 );
351 
352  // B.2) Setup the DAE's outArgs for f(i) and/or W(i,i) ...
353  if (!is_null(f_bar))
354  daeOutArgs.set_f( f_bar->getNonconstVectorBlock(i) );
355  if (!is_null(W_op_bar))
356  daeOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,i).assert_not_null());
357 
358  // B.3) Compute f_bar(i) and/or W_op_bar(i,i) ...
359  daeModel_->evalModel( daeInArgs, daeOutArgs );
360  daeOutArgs.set_f(Teuchos::null);
361  daeOutArgs.set_W_op(Teuchos::null);
362 
363  // B.4) Evaluate W_op_bar(i,i-1)
364  if ( !is_null(W_op_bar) && i > 0 ) {
365  daeInArgs.set_alpha( -oneOverDeltaT );
366  daeInArgs.set_beta( 0.0 );
367  daeOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,i-1).assert_not_null());
368  daeModel_->evalModel( daeInArgs, daeOutArgs );
369  daeOutArgs.set_W_op(Teuchos::null);
370  }
371 
372  //
373  t_i += delta_t_;
374 
375  }
376 
377 /*
378  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
379 */
380 
381 }
382 
383 
384 } // namespace Rythmos
385 
386 
387 #endif // RYTHMOS_TIME_DISCRETIZED_BACKWARD_EULER_MODEL_EVALUATOR_HPP
RCP< TimeDiscretizedBackwardEulerModelEvaluator< Scalar > > timeDiscretizedBackwardEulerModelEvaluator(const RCP< const Thyra::ModelEvaluator< Scalar > > &daeModel, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &initCond, const Scalar finalTime, const int numTimeSteps, const RCP< Thyra::LinearOpWithSolveFactoryBase< Scalar > > &W_bar_factory)
Non-member constructor.
RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
void initialize(const RCP< const Thyra::ModelEvaluator< Scalar > > &daeModel, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &initCond, const Scalar finalTime, const int numTimeSteps, const RCP< Thyra::LinearOpWithSolveFactoryBase< Scalar > > &W_bar_factory=Teuchos::null)