Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultStateEliminationModelEvaluator.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) 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// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef THYRA_DEFAULT_STATE_ELIMINATION_MODEL_EVALUATOR_HPP
43#define THYRA_DEFAULT_STATE_ELIMINATION_MODEL_EVALUATOR_HPP
44
45#include "Thyra_ModelEvaluatorDelegatorBase.hpp"
46#include "Thyra_DefaultNominalBoundsOverrideModelEvaluator.hpp"
47#include "Thyra_NonlinearSolverBase.hpp"
48#include "Teuchos_Time.hpp"
49
50
51namespace Thyra {
52
53
61template<class Scalar>
63 : virtual public ModelEvaluatorDelegatorBase<Scalar>
64{
65public:
66
69
72
75 const Teuchos::RCP<ModelEvaluator<Scalar> > &thyraModel,
76 const Teuchos::RCP<NonlinearSolverBase<Scalar> > &stateSolver
77 );
78
80 void initialize(
81 const Teuchos::RCP<ModelEvaluator<Scalar> > &thyraModel,
82 const Teuchos::RCP<NonlinearSolverBase<Scalar> > &stateSolver
83 );
84
86 void uninitialize(
87 Teuchos::RCP<ModelEvaluator<Scalar> > *thyraModel = NULL,
88 Teuchos::RCP<NonlinearSolverBase<Scalar> > *stateSolver = NULL
89 );
90
92
95
97 std::string description() const;
98
100
103
120
122
123private:
124
127
129 ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
131 void evalModelImpl(
134 ) const;
135
137
138
139private:
140
143
145
146 mutable Teuchos::RCP<VectorBase<Scalar> > x_guess_solu_;
147
148};
149
150// /////////////////////////////////
151// Implementations
152
153// Constructors/initializers/accessors/utilities
154
155template<class Scalar>
157{}
158
159template<class Scalar>
161 const Teuchos::RCP<ModelEvaluator<Scalar> > &thyraModel
162 ,const Teuchos::RCP<NonlinearSolverBase<Scalar> > &stateSolver
163 )
164{
165 initialize(thyraModel,stateSolver);
166}
167
168template<class Scalar>
170 const Teuchos::RCP<ModelEvaluator<Scalar> > &thyraModel
171 ,const Teuchos::RCP<NonlinearSolverBase<Scalar> > &stateSolver
172 )
173{
175 TEUCHOS_TEST_FOR_EXCEPT(!stateSolver.get());
176 stateSolver_ = stateSolver;
177 x_guess_solu_ = Teuchos::null; // We will get the guess at the last possible moment!
178 wrappedThyraModel_ = Teuchos::rcp(
180 Teuchos::rcp_const_cast<ModelEvaluator<Scalar> >(thyraModel)
181 ,Teuchos::null
182 )
183 );
184 stateSolver_->setModel(wrappedThyraModel_);
185}
186
187template<class Scalar>
191 )
192{
193 if(thyraModel) *thyraModel = this->getUnderlyingModel();
194 if(stateSolver) *stateSolver = stateSolver_;
196 stateSolver_ = Teuchos::null;
197 wrappedThyraModel_ = Teuchos::null;
198 x_guess_solu_ = Teuchos::null;
199}
200
201// Public functions overridden from Teuchos::Describable
202
203
204template<class Scalar>
206{
208 thyraModel = this->getUnderlyingModel();
209 std::ostringstream oss;
210 oss << "Thyra::DefaultStateEliminationModelEvaluator{";
211 oss << "thyraModel=";
212 if(thyraModel.get())
213 oss << "\'"<<thyraModel->description()<<"\'";
214 else
215 oss << "NULL";
216 oss << ",stateSolver=";
217 if(stateSolver_.get())
218 oss << "\'"<<stateSolver_->description()<<"\'";
219 else
220 oss << "NULL";
221 oss << "}";
222 return oss.str();
223}
224
225// Public functions overridden from ModelEvaulator
226
227template<class Scalar>
230{
231 return Teuchos::null;
232}
233
234template<class Scalar>
237{
238 return Teuchos::null;
239}
240
241template<class Scalar>
244{
245 typedef ModelEvaluatorBase MEB;
247 thyraModel = this->getUnderlyingModel();
248 MEB::InArgsSetup<Scalar> nominalValues(thyraModel->getNominalValues());
249 nominalValues.setModelEvalDescription(this->description());
250 nominalValues.setUnsupportsAndRelated(MEB::IN_ARG_x); // Wipe out x, x_dot ...
251 return nominalValues;
252}
253
254template<class Scalar>
257{
258 typedef ModelEvaluatorBase MEB;
260 thyraModel = this->getUnderlyingModel();
261 MEB::InArgsSetup<Scalar> lowerBounds(thyraModel->getLowerBounds());
262 lowerBounds.setModelEvalDescription(this->description());
263 lowerBounds.setUnsupportsAndRelated(MEB::IN_ARG_x); // Wipe out x, x_dot ...
264 return lowerBounds;
265}
266
267template<class Scalar>
270{
271 typedef ModelEvaluatorBase MEB;
273 thyraModel = this->getUnderlyingModel();
274 MEB::InArgsSetup<Scalar> upperBounds(thyraModel->getUpperBounds());
275 upperBounds.setModelEvalDescription(this->description());
276 upperBounds.setUnsupportsAndRelated(MEB::IN_ARG_x); // Wipe out x, x_dot ...
277 return upperBounds;
278}
279
280template<class Scalar>
283{
284 return Teuchos::null;
285}
286
287template<class Scalar>
290{
291 return Teuchos::null;
292}
293
294
295template<class Scalar>
298{
299 typedef ModelEvaluatorBase MEB;
301 thyraModel = this->getUnderlyingModel();
302 const MEB::InArgs<Scalar> wrappedInArgs = thyraModel->createInArgs();
303 MEB::InArgsSetup<Scalar> inArgs;
304 inArgs.setModelEvalDescription(this->description());
305 inArgs.set_Np(wrappedInArgs.Np());
306 inArgs.setSupports(wrappedInArgs);
307 inArgs.setUnsupportsAndRelated(MEB::IN_ARG_x); // Wipe out x, x_dot ...
308 return inArgs;
309}
310
311
312// Private functions overridden from ModelEvaulatorDefaultBase
313
314
315template<class Scalar>
318{
319 typedef ModelEvaluatorBase MEB;
321 thyraModel = this->getUnderlyingModel();
322 const MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
323 const int Np = wrappedOutArgs.Np(), Ng = wrappedOutArgs.Ng();
324 MEB::OutArgsSetup<Scalar> outArgs;
325 outArgs.setModelEvalDescription(this->description());
326 outArgs.set_Np_Ng(Np,Ng);
327 outArgs.setSupports(wrappedOutArgs);
328 outArgs.setUnsupportsAndRelated(MEB::IN_ARG_x); // wipe out DgDx ...
329 outArgs.setUnsupportsAndRelated(MEB::OUT_ARG_f); // wipe out f, W, DfDp ...
330 return outArgs;
331}
332
333template<class Scalar>
334void DefaultStateEliminationModelEvaluator<Scalar>::evalModelImpl(
335 const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
336 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
337 ) const
338{
339 typedef ModelEvaluatorBase MEB;
340 using Teuchos::rcp;
341 using Teuchos::rcp_const_cast;
342 using Teuchos::rcp_dynamic_cast;
343 using Teuchos::OSTab;
344 using Teuchos::as;
345
346 Teuchos::Time totalTimer(""), timer("");
347 totalTimer.start(true);
348
349 const Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
350 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
351 Teuchos::OSTab tab(out);
352 if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
353 *out << "\nEntering Thyra::DefaultStateEliminationModelEvaluator<Scalar>::evalModel(...) ...\n";
354
356 thyraModel = this->getUnderlyingModel();
357
358 const int Np = outArgs.Np(), Ng = outArgs.Ng();
359
360 // Get the intial state guess if not already gotten
361 if (is_null(x_guess_solu_)) {
362 const ModelEvaluatorBase::InArgs<Scalar>
363 nominalValues = thyraModel->getNominalValues();
364 if(nominalValues.get_x().get()) {
365 x_guess_solu_ = nominalValues.get_x()->clone_v();
366 }
367 else {
368 x_guess_solu_ = createMember(thyraModel->get_x_space());
369 assign(x_guess_solu_.ptr(), as<Scalar>(0.0));
370 }
371 }
372
373 // Reset the nominal values
374 MEB::InArgs<Scalar> wrappedNominalValues = thyraModel->getNominalValues();
375 wrappedNominalValues.setArgs(inArgs,true);
376 wrappedNominalValues.set_x(x_guess_solu_);
377
379 //VOTSME thyraModel_outputTempState(rcp(&wrappedThyraModel,false),out,verbLevel);
380
382 VOTSNSB statSolver_outputTempState(
383 stateSolver_,out
384 ,static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW) ? Teuchos::VERB_LOW : Teuchos::VERB_NONE
385 );
386
387 if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_EXTREME))
388 *out
389 << "\ninArgs =\n" << Teuchos::describe(inArgs,verbLevel)
390 << "\noutArgs on input =\n" << Teuchos::describe(outArgs,Teuchos::VERB_LOW);
391
392 if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
393 *out << "\nSolving f(x,...) for x ...\n";
394
395 wrappedThyraModel_->setNominalValues(
396 rcp(new MEB::InArgs<Scalar>(wrappedNominalValues))
397 );
398
399 SolveStatus<Scalar> solveStatus = stateSolver_->solve(&*x_guess_solu_,NULL);
400
401 if( solveStatus.solveStatus == SOLVE_STATUS_CONVERGED ) {
402
403 if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
404 *out << "\nComputing the output functions at the solved state solution ...\n";
405
406 MEB::InArgs<Scalar> wrappedInArgs = thyraModel->createInArgs();
407 MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
408 wrappedInArgs.setArgs(inArgs,true);
409 wrappedInArgs.set_x(x_guess_solu_);
410 wrappedOutArgs.setArgs(outArgs,true);
411
412 for( int l = 0; l < Np; ++l ) {
413 for( int j = 0; j < Ng; ++j ) {
414 if(
415 outArgs.supports(MEB::OUT_ARG_DgDp,j,l).none()==false
416 && outArgs.get_DgDp(j,l).isEmpty()==false
417 )
418 {
419 // Set DfDp(l) and DgDx(j) to be computed!
420 //wrappedOutArgs.set_DfDp(l,...);
421 //wrappedOutArgs.set_DgDx(j,...);
423 }
424 }
425 }
426
427 thyraModel->evalModel(wrappedInArgs,wrappedOutArgs);
428
429 //
430 // Compute DgDp(j,l) using direct sensitivties
431 //
432 for( int l = 0; l < Np; ++l ) {
433 if(
434 wrappedOutArgs.supports(MEB::OUT_ARG_DfDp,l).none()==false
435 && wrappedOutArgs.get_DfDp(l).isEmpty()==false
436 )
437 {
438 //
439 // Compute: D(l) = -inv(DfDx)*DfDp(l)
440 //
442 for( int j = 0; j < Ng; ++j ) {
443 if(
444 outArgs.supports(MEB::OUT_ARG_DgDp,j,l).none()==false
445 && outArgs.get_DgDp(j,l).isEmpty()==false
446 )
447 {
448 //
449 // Compute: DgDp(j,l) = DgDp(j,l) + DgDx(j)*D
450 //
452 }
453 }
454 }
455 }
456 // ToDo: Add a mode to compute DgDp(l) using adjoint sensitivities?
457
458 }
459 else {
460
461 if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
462 *out << "\nFailed to converge, returning NaNs ...\n";
463 outArgs.setFailed();
464
465 }
466
467 if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_EXTREME))
468 *out
469 << "\noutArgs on output =\n" << Teuchos::describe(outArgs,verbLevel);
470
471 totalTimer.stop();
472 if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
473 *out
474 << "\nTotal evaluation time = "<<totalTimer.totalElapsedTime()<<" sec\n"
475 << "\nLeaving Thyra::DefaultStateEliminationModelEvaluator<Scalar>::evalModel(...) ...\n";
476
477}
478
479
480} // namespace Thyra
481
482
483#endif // THYRA_DEFAULT_STATE_ELIMINATION_MODEL_EVALUATOR_HPP
T * get() const
This class wraps any ModelEvaluator object and allows the client to overide the state contained in th...
This class wraps any ModelEvaluator object along with a NonlinearSolverBase object and eliminates the...
Teuchos::RCP< LinearOpWithSolveBase< Scalar > > create_W() const
Teuchos::RCP< const VectorSpaceBase< Scalar > > get_x_space() const
void initialize(const Teuchos::RCP< ModelEvaluator< Scalar > > &thyraModel, const Teuchos::RCP< NonlinearSolverBase< Scalar > > &stateSolver)
Teuchos::RCP< const VectorSpaceBase< Scalar > > get_f_space() const
Concrete aggregate class for all input arguments computable by a ModelEvaluator subclass object.
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object.
Base subclass for ModelEvaluator that defines some basic types.
This is a base class that delegetes almost all function to a wrapped model evaluator object.
void initialize(const RCP< ModelEvaluator< Scalar > > &model)
Initialize given a non-const model evaluator.
Pure abstract base interface for evaluating a stateless "model" that can be mapped into a number of d...
Base class for all nonlinear equation solvers.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
@ SOLVE_STATUS_CONVERGED
The requested solution criteria has likely been achieved.
TypeTo as(const TypeFrom &t)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)