Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_AuxiliaryIntegralModelEvaluator_impl.hpp
Go to the documentation of this file.
1// @HEADER
2// ****************************************************************************
3// Tempus: Copyright (2017) Sandia Corporation
4//
5// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6// ****************************************************************************
7// @HEADER
8
9#ifndef Tempus_AuxiliaryIntegralModelEvaluator_impl_hpp
10#define Tempus_AuxiliaryIntegralModelEvaluator_impl_hpp
11
14#include "Thyra_VectorStdOps.hpp"
15
16namespace Tempus {
17
18template <typename Scalar>
21 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & model,
22 const int g_index) :
23 model_(model),
24 g_index_(g_index),
25 t_interp_(Teuchos::ScalarTraits<Scalar>::rmax())
26{
27 typedef Thyra::ModelEvaluatorBase MEB;
28
29 space_ = model_->get_g_space(g_index);
30
31 MEB::InArgs<Scalar> me_inArgs = model_->createInArgs();
32 MEB::InArgsSetup<Scalar> inArgs;
33 inArgs.setModelEvalDescription(this->description());
34 inArgs.setSupports(MEB::IN_ARG_x);
35 inArgs.setSupports(MEB::IN_ARG_t);
36 if (me_inArgs.supports(MEB::IN_ARG_x_dot))
37 inArgs.setSupports(MEB::IN_ARG_x_dot);
38 inArgs.setSupports(MEB::IN_ARG_alpha);
39 inArgs.setSupports(MEB::IN_ARG_beta);
40 inArgs.set_Np(me_inArgs.Np());
41 prototypeInArgs_ = inArgs;
42
43 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
44 MEB::OutArgsSetup<Scalar> outArgs;
45 outArgs.setModelEvalDescription(this->description());
46 outArgs.set_Np_Ng(me_inArgs.Np(),0);
47 outArgs.setSupports(MEB::OUT_ARG_f);
48 outArgs.setSupports(MEB::OUT_ARG_W_op);
49 prototypeOutArgs_ = outArgs;
50}
51
52template <typename Scalar>
53void
56 const Teuchos::RCP<const Tempus::SolutionHistory<Scalar> >& sh)
57{
58 sh_ = sh;
59 t_interp_ = Teuchos::ScalarTraits<Scalar>::rmax();
60 forward_state_ = Teuchos::null;
61}
62
63template <typename Scalar>
64Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
66get_p_space(int p) const
67{
68 TEUCHOS_ASSERT(p < model_->Np());
69 return model_->get_p_space(p);
70}
71
72template <typename Scalar>
73Teuchos::RCP<const Teuchos::Array<std::string> >
75get_p_names(int p) const
76{
77 TEUCHOS_ASSERT(p < model_->Np());
78 return model_->get_p_names(p);
79}
80
81template <typename Scalar>
82Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
84get_x_space() const
85{
86 return space_;
87}
88
89template <typename Scalar>
90Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
92get_f_space() const
93{
94 return space_;
95}
96
97template <typename Scalar>
98Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
100create_W_op() const
101{
102 return Thyra::scaledIdentity(space_, Scalar(1.0));
103}
104
105template <typename Scalar>
106Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
108get_W_factory() const
109{
110 return Thyra::scaledIdentitySolveFactory(space_, Scalar(1.0));
111}
112
113template <typename Scalar>
114Thyra::ModelEvaluatorBase::InArgs<Scalar>
116createInArgs() const
117{
118 return prototypeInArgs_;
119}
120
121template <typename Scalar>
122Thyra::ModelEvaluatorBase::InArgs<Scalar>
124getNominalValues() const
125{
126 typedef Thyra::ModelEvaluatorBase MEB;
127 using Teuchos::RCP;
128 using Teuchos::rcp_dynamic_cast;
129
130 MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
131 MEB::InArgs<Scalar> nominal = this->createInArgs();
132
133 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
134
135 // Set initial x, x_dot
136 RCP< Thyra::VectorBase<Scalar> > x = Thyra::createMember(*space_);
137 Thyra::assign(x.ptr(), zero);
138 nominal.set_x(x);
139
140 if (me_nominal.supports(MEB::IN_ARG_x_dot)) {
141 RCP< Thyra::VectorBase<Scalar> > x_dot = Thyra::createMember(*space_);
142 Thyra::assign(x_dot.ptr(), zero);
143 nominal.set_x_dot(x_dot);
144 }
145
146 const int np = model_->Np();
147 for (int i=0; i<np; ++i)
148 nominal.set_p(i, me_nominal.get_p(i));
149
150 return nominal;
151}
152
153template <typename Scalar>
154Thyra::ModelEvaluatorBase::OutArgs<Scalar>
156createOutArgsImpl() const
157{
158 return prototypeOutArgs_;
159}
160
161template <typename Scalar>
162void
164evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
165 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
166{
167 typedef Thyra::ModelEvaluatorBase MEB;
168 using Teuchos::RCP;
169 using Teuchos::rcp_dynamic_cast;
170
171 // Interpolate forward solution at supplied time, reusing previous
172 // interpolation if possible
173 TEUCHOS_ASSERT(sh_ != Teuchos::null);
174 const Scalar t = inArgs.get_t();
175 if (forward_state_ == Teuchos::null || t_interp_ != t) {
176 if (forward_state_ == Teuchos::null)
177 forward_state_ = sh_->interpolateState(t);
178 else
179 sh_->interpolateState(t, forward_state_.get());
180 t_interp_ = t;
181 }
182
183 // Residual
184 RCP<Thyra::VectorBase<Scalar> > f = outArgs.get_f();
185 if (f != Teuchos::null) {
186 MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
187 me_inArgs.set_x(forward_state_->getX());
188 if (me_inArgs.supports(MEB::IN_ARG_x_dot))
189 me_inArgs.set_x_dot(forward_state_->getXDot());
190 if (me_inArgs.supports(MEB::IN_ARG_t))
191 me_inArgs.set_t(t);
192 const int np = me_inArgs.Np();
193 for (int i=0; i<np; ++i)
194 me_inArgs.set_p(i, inArgs.get_p(i));
195
196 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
197 me_outArgs.set_g(g_index_, f);
198
199 model_->evalModel(me_inArgs, me_outArgs);
200
201 // For explicit form, f = g and we are done
202 // For implicit form, f = dz/dt - g
203 if (inArgs.supports(MEB::IN_ARG_x_dot)) {
204 RCP<const Thyra::VectorBase<Scalar> > x_dot = inArgs.get_x_dot();
205 if (x_dot != Teuchos::null)
206 Thyra::V_VmV(f.ptr(), *x_dot, *f);
207 }
208 }
209
210 // W = alpha*df/z_dot + beta*df/dz = alpha * I
211 RCP<Thyra::LinearOpBase<Scalar> > op = outArgs.get_W_op();
212 if (op != Teuchos::null) {
213 const Scalar alpha = inArgs.get_alpha();
214 RCP<Thyra::ScaledIdentityLinearOpWithSolve<Scalar> > si_op =
215 rcp_dynamic_cast<Thyra::ScaledIdentityLinearOpWithSolve<Scalar> >(op);
216 si_op->setScale(alpha);
217 }
218}
219
220} // namespace Tempus
221
222#endif
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
AuxiliaryIntegralModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const int g_index)
Constructor.
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int p) const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int p) const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > space_
Thyra::ModelEvaluatorBase::InArgs< Scalar > prototypeInArgs_
Thyra::ModelEvaluatorBase::OutArgs< Scalar > prototypeOutArgs_
void setForwardSolutionHistory(const Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > &sh)
Set solution history from forward evaluation.
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > model_
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...