Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_WrapperModelEvaluatorSecondOrder_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_WrapperModelEvaluatorSecondOrder_impl_hpp
10#define Tempus_WrapperModelEvaluatorSecondOrder_impl_hpp
11
12namespace Tempus {
13
14
15template <typename Scalar>
16Thyra::ModelEvaluatorBase::InArgs<Scalar>
18createInArgs() const
19{
20#ifdef VERBOSE_DEBUG_OUTPUT
21 *out_ << "DEBUG: " << __PRETTY_FUNCTION__ << "\n";
22#endif
23 typedef Thyra::ModelEvaluatorBase MEB;
24
25 MEB::InArgsSetup<Scalar> inArgs;
26 inArgs.setModelEvalDescription(this->description());
27 inArgs.set_Np(appModel_->Np());
28 inArgs.setSupports(MEB::IN_ARG_x);
29
30 return std::move(inArgs);
31}
32
33
34template <typename Scalar>
35Thyra::ModelEvaluatorBase::OutArgs<Scalar>
38{
39#ifdef VERBOSE_DEBUG_OUTPUT
40 *out_ << "DEBUG: " << __PRETTY_FUNCTION__ << "\n";
41#endif
42 typedef Thyra::ModelEvaluatorBase MEB;
43
44 MEB::OutArgsSetup<Scalar> outArgs;
45 outArgs.setModelEvalDescription(this->description());
46 outArgs.set_Np_Ng(appModel_->Np(),0);
47 outArgs.setSupports(MEB::OUT_ARG_f);
48 outArgs.setSupports(MEB::OUT_ARG_W_op);
49
50 return std::move(outArgs);
51}
52
53
54template <typename Scalar>
55void
57evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
58 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
59{
60#ifdef VERBOSE_DEBUG_OUTPUT
61 *out_ << "DEBUG: " << __PRETTY_FUNCTION__ << "\n";
62#endif
63 typedef Thyra::ModelEvaluatorBase MEB;
64 using Teuchos::RCP;
65
66 //Setup initial condition
67 //Create and populate inArgs
68 MEB::InArgs<Scalar> appInArgs = appModel_->createInArgs();
69 MEB::OutArgs<Scalar> appOutArgs = appModel_->createOutArgs();
70
71 switch (schemeType_)
72 {
73 case NEWMARK_IMPLICIT_AFORM: {
74 // Specific for the Newmark Implicit a-Form stepper. May want to
75 // redesign this for a generic second-order scheme to not have an
76 // if statement here...
77 // IKT, 3/14/17: this is effectively the same as the
78 // Piro::NewmarkDecorator::evalModel function.
79
80 // The solution variable in NOX is the acceleration, a_{n+1}
81 appInArgs.set_x_dot_dot(inArgs.get_x());
82
83 // Compute the velocity.
84 // v_n = velocity_{pred} + \gamma dt a_n
85 auto velocity = Thyra::createMember(inArgs.get_x()->space());
86 Thyra::V_StVpStV(velocity.ptr(),
87 Scalar(1.0), *v_pred_, delta_t_*gamma_, *inArgs.get_x());
88 appInArgs.set_x_dot(velocity);
89
90 // Compute the displacement.
91 // d_n = displacement_{pred} + \beta dt^2 a_n
92 auto displacement = Thyra::createMember(inArgs.get_x()->space());
93 Thyra::V_StVpStV(displacement.ptr(),
94 Scalar(1.0), *d_pred_, beta_*delta_t_*delta_t_, *inArgs.get_x());
95 appInArgs.set_x(displacement);
96
97 appInArgs.set_W_x_dot_dot_coeff(Scalar(1.0)); // da/da
98 appInArgs.set_alpha(gamma_*delta_t_); // dv/da
99 appInArgs.set_beta(beta_*delta_t_*delta_t_); // dd/da
100
101 appInArgs.set_t(t_);
102 for (int i=0; i<appModel_->Np(); ++i) {
103 if (inArgs.get_p(i) != Teuchos::null)
104 appInArgs.set_p(i, inArgs.get_p(i));
105 }
106
107 //Setup output condition
108 //Create and populate outArgs
109 appOutArgs.set_f(outArgs.get_f());
110 appOutArgs.set_W_op(outArgs.get_W_op());
111
112 // build residual and jacobian
113 appModel_->evalModel(appInArgs,appOutArgs);
114 break;
115 }
116
117 case NEWMARK_IMPLICIT_DFORM: {
118 // Setup initial condition
119 // Populate inArgs
120 RCP<Thyra::VectorBase<Scalar> const>
121 d = inArgs.get_x();
122
123 RCP<Thyra::VectorBase<Scalar>>
124 v = Thyra::createMember(inArgs.get_x()->space());
125
126 RCP<Thyra::VectorBase<Scalar>>
127 a = Thyra::createMember(inArgs.get_x()->space());
128
129#ifdef DEBUG_OUTPUT
130 Teuchos::Range1D range;
131
132 *out_ << "\n*** d_bef ***\n";
133 RTOpPack::ConstSubVectorView<Scalar> dov;
134 d->acquireDetachedView(range, &dov);
135 auto doa = dov.values();
136 for (auto i = 0; i < doa.size(); ++i) *out_ << doa[i] << " ";
137 *out_ << "\n*** d_bef ***\n";
138
139 *out_ << "\n*** v_bef ***\n";
140 RTOpPack::ConstSubVectorView<Scalar> vov;
141 v->acquireDetachedView(range, &vov);
142 auto voa = vov.values();
143 for (auto i = 0; i < voa.size(); ++i) *out_ << voa[i] << " ";
144 *out_ << "\n*** v_bef ***\n";
145
146 *out_ << "\n*** a_bef ***\n";
147 RTOpPack::ConstSubVectorView<Scalar> aov;
148 a->acquireDetachedView(range, &aov);
149 auto aoa = aov.values();
150 for (auto i = 0; i < aoa.size(); ++i) *out_ << aoa[i] << " ";
151 *out_ << "\n*** a_bef ***\n";
152#endif
153
154 Scalar const
155 c = 1.0 / beta_ / delta_t_ / delta_t_;
156
157 // compute acceleration
158 // a_{n+1} = (d_{n+1} - d_pred) / dt / dt / beta
159 Thyra::V_StVpStV(Teuchos::ptrFromRef(*a), c, *d, -c, *d_pred_);
160
161 // compute velocity
162 // v_{n+1} = v_pred + \gamma dt a_{n+1}
163 Thyra::V_StVpStV(
164 Teuchos::ptrFromRef(*v), 1.0, *v_pred_, delta_t_ * gamma_, *a);
165
166 appInArgs.set_x(d);
167 appInArgs.set_x_dot(v);
168 appInArgs.set_x_dot_dot(a);
169
170 appInArgs.set_W_x_dot_dot_coeff(c); // da/dd
171 appInArgs.set_alpha(gamma_ / delta_t_ / beta_); // dv/dd
172 appInArgs.set_beta(1.0); // dd/dd
173
174 appInArgs.set_t(t_);
175 for (int i = 0; i < appModel_->Np(); ++i) {
176 if (inArgs.get_p(i) != Teuchos::null)
177 appInArgs.set_p(i, inArgs.get_p(i));
178 }
179
180 // Setup output condition
181 // Create and populate outArgs
182 appOutArgs.set_f(outArgs.get_f());
183 appOutArgs.set_W_op(outArgs.get_W_op());
184
185 // build residual and jacobian
186 appModel_->evalModel(appInArgs, appOutArgs);
187
188 // compute acceleration
189 // a_{n+1} = (d_{n+1} - d_pred) / dt / dt / beta
190 Thyra::V_StVpStV(Teuchos::ptrFromRef(*a), c, *d, -c, *d_pred_);
191
192 // compute velocity
193 // v_{n+1} = v_pred + \gamma dt a_{n+1}
194 Thyra::V_StVpStV(
195 Teuchos::ptrFromRef(*v), 1.0, *v_pred_, delta_t_ * gamma_, *a);
196
197 appInArgs.set_x(d);
198 appInArgs.set_x_dot(v);
199 appInArgs.set_x_dot_dot(a);
200
201#ifdef DEBUG_OUTPUT
202 *out_ << "\n*** d_aft ***\n";
203 RTOpPack::ConstSubVectorView<Scalar> dnv;
204 d->acquireDetachedView(range, &dnv);
205 auto dna = dnv.values();
206 for (auto i = 0; i < dna.size(); ++i) *out_ << dna[i] << " ";
207 *out_ << "\n*** d_aft ***\n";
208
209 *out_ << "\n*** v_aft ***\n";
210 RTOpPack::ConstSubVectorView<Scalar> vnv;
211 v->acquireDetachedView(range, &vnv);
212 auto vna = vnv.values();
213 for (auto i = 0; i < vna.size(); ++i) *out_ << vna[i] << " ";
214 *out_ << "\n*** v_aft ***\n";
215
216 *out_ << "\n*** a_aft ***\n";
217 RTOpPack::ConstSubVectorView<Scalar> anv;
218 a->acquireDetachedView(range, &anv);
219 auto ana = anv.values();
220 for (auto i = 0; i < ana.size(); ++i) *out_ << ana[i] << " ";
221 *out_ << "\n*** a_aft ***\n";
222#endif
223 break;
224 }
225 }
226}
227
228
229} // namespace Tempus
230
231#endif // Tempus_WrapperModelEvaluatorSecondOrder_impl_hpp
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const