Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
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
42namespace Rythmos {
43
44
49template<class Scalar>
51 : virtual public Thyra::StateFuncModelEvaluatorBase<Scalar>
52{
53public:
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
95private:
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
110private:
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
131template<class Scalar>
132RCP<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> >
142 model(new 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
155template<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
167template<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
211template<class Scalar>
212RCP<const Thyra::VectorSpaceBase<Scalar> >
217
218
219template<class Scalar>
220RCP<const Thyra::VectorSpaceBase<Scalar> >
225
226
227template<class Scalar>
228RCP<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
245template<class Scalar>
246RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
248{
249 return W_bar_factory_;
250}
251
252
253template<class Scalar>
254Thyra::ModelEvaluatorBase::InArgs<Scalar>
256{
257 typedef Thyra::ModelEvaluatorBase MEB;
258 TEUCHOS_TEST_FOR_EXCEPT(true);
259 return MEB::InArgs<Scalar>();
260}
261
262
263template<class Scalar>
264Thyra::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
278template<class Scalar>
279Thyra::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
293template<class Scalar>
294void 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.
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)
RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const