Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
HarmonicOscillatorModel_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_TEST_HARMONIC_OSCILLATOR_MODEL_IMPL_HPP
10#define TEMPUS_TEST_HARMONIC_OSCILLATOR_MODEL_IMPL_HPP
11
12#include "Teuchos_StandardParameterEntryValidators.hpp"
13
14#include "Thyra_DefaultSpmdVectorSpace.hpp"
15#include "Thyra_DetachedVectorView.hpp"
16#include "Thyra_DetachedMultiVectorView.hpp"
17#include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
18#include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
19#include "Thyra_DefaultLinearOpSource.hpp"
20#include "Thyra_VectorStdOps.hpp"
21#include <iostream>
22
23
24namespace Tempus_Test {
25
26template<class Scalar>
28HarmonicOscillatorModel(Teuchos::RCP<Teuchos::ParameterList> pList_, const bool use_accel_IC):
29 out_(Teuchos::VerboseObjectBase::getDefaultOStream())
30{
31 isInitialized_ = false;
32 setParameterList(pList_);
33 *out_ << "\n\nDamping coeff c = " << c_ << "\n";
34 *out_ << "Forcing coeff f = " << f_ << "\n";
35 *out_ << "x coeff k = " << k_ << "\n";
36 *out_ << "Mass coeff m = " << m_ << "\n";
37 //Divide all coefficients by m_
38 k_ /= m_;
39 f_ /= m_;
40 c_ /= m_;
41 m_ = 1.0;
42 //Set up space and initial guess for solution vector
43 vecLength_ = 1;
44 x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(vecLength_);
45 x_vec_ = createMember(x_space_);
46 Thyra::put_scalar(0.0, x_vec_.ptr());
47 x_dot_vec_ = createMember(x_space_);
48 Thyra::put_scalar(1.0, x_dot_vec_.ptr());
49 x_dot_dot_vec_ = createMember(x_space_);
50 //The following is the initial condition for the acceleration
51 //Commenting this out to check that IC for acceleration
52 //is computed correctly using displacement and velocity ICs
53 //inside 2nd order steppers.
54 //Thyra::put_scalar(f_-c_, x_dot_dot_vec_.ptr());
55 if (use_accel_IC == true) {
56 Thyra::put_scalar(-2.0, x_dot_dot_vec_.ptr());
57 }
58 else {
59 //Instead of real IC, putting arbitrary, incorrect IC to check correctness
60 //in stepper involving calculation of a IC.
61 Thyra::put_scalar(7.0, x_dot_dot_vec_.ptr());
62 }
63
64 //Set up responses
65 numResponses_ = 1;
66 g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(numResponses_);
67
69}
70
71template<class Scalar>
72Thyra::ModelEvaluatorBase::InArgs<Scalar>
74getExactSolution(double t) const
75{
77 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
78 "Error, setupInOutArgs_ must be called first!\n");
79 Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs = inArgs_;
80 double exact_t = t;
81 inArgs.set_t(exact_t);
82 Teuchos::RCP<VectorBase<Scalar> > exact_x = createMember(x_space_);
83 { // scope to delete DetachedVectorView
84 Thyra::DetachedVectorView<Scalar> exact_x_view(*exact_x);
85 if (k_ == 0) {
86 if (c_ == 0)
87 exact_x_view[0] = t*(1.0+0.5*f_*t);
88 else
89 exact_x_view[0] = (c_-f_)/(c_*c_)*(1.0-exp(-c_*t))
90 + f_*t/c_;
91 }
92 else {
93 exact_x_view[0] = 1.0/sqrt(k_)*sin(sqrt(k_)*t) + f_/k_*(1.0-cos(sqrt(k_)*t));
94 }
95 }
96 inArgs.set_x(exact_x);
97 Teuchos::RCP<VectorBase<Scalar> > exact_x_dot = createMember(x_space_);
98 { // scope to delete DetachedVectorView
99 Thyra::DetachedVectorView<Scalar> exact_x_dot_view(*exact_x_dot);
100 if (k_ == 0) {
101 if (c_ == 0)
102 exact_x_dot_view[0] = 1.0+f_*t;
103 else
104 exact_x_dot_view[0] = (c_-f_)/c_*exp(-c_*t)+f_/c_;
105 }
106 else {
107 exact_x_dot_view[0] = cos(sqrt(k_)*t) + f_/sqrt(k_)*sin(sqrt(k_)*t);
108 }
109 }
110 inArgs.set_x_dot(exact_x_dot);
111 Teuchos::RCP<VectorBase<Scalar> > exact_x_dot_dot = createMember(x_space_);
112 { // scope to delete DetachedVectorView
113 Thyra::DetachedVectorView<Scalar> exact_x_dot_dot_view(*exact_x_dot_dot);
114 if (k_ == 0) {
115 if (c_ == 0)
116 exact_x_dot_dot_view[0] = f_;
117 else
118 exact_x_dot_dot_view[0] = (f_-c_)*exp(-c_*t);
119 }
120 else {
121 exact_x_dot_dot_view[0] = f_*cos(sqrt(k_)*t) - sqrt(k_)*sin(sqrt(k_)*t);
122 }
123 }
124 inArgs.set_x_dot_dot(exact_x_dot_dot);
125 return(inArgs);
126}
127
128template<class Scalar>
129Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
131get_x_space() const
132{
133 return x_space_;
134}
135
136
137template<class Scalar>
138Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
140get_f_space() const
141{
142 return x_space_;
143}
144
145
146template<class Scalar>
147Thyra::ModelEvaluatorBase::InArgs<Scalar>
149getNominalValues() const
150{
151 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
152 "Error, setupInOutArgs_ must be called first!\n");
153 return nominalValues_;
154}
155
156
157template<class Scalar>
158Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
160create_W() const
161{
162 Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = this->get_W_factory();
163 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
164 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > matrix_mv = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,true);
165 Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix_mv );
166 //IKT: is it necessary for W to be non-singular when initialized?
167 matrix_view(0,0) = 1.0;
168 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > W =
169 Thyra::linearOpWithSolve<Scalar>(*W_factory, matrix );
170 return W;
171}
172
173
174template<class Scalar>
175Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
177create_W_op() const
178{
179 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > matrix = Thyra::createMembers(x_space_, vecLength_);
180 return(matrix);
181}
182
183
184template<class Scalar>
185Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
187get_W_factory() const
188{
189 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory =
190 Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
191 return W_factory;
192}
193
194
195template<class Scalar>
196Thyra::ModelEvaluatorBase::InArgs<Scalar>
198createInArgs() const
199{
200 setupInOutArgs_();
201 return inArgs_;
202}
203
204
205// Private functions overridden from ModelEvaluatorDefaultBase
206
207
208template<class Scalar>
209Thyra::ModelEvaluatorBase::OutArgs<Scalar>
211createOutArgsImpl() const
212{
213 setupInOutArgs_();
214 return outArgs_;
215}
216
217
218template<class Scalar>
219void
222 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
223 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
224 ) const
225{
226 using Teuchos::RCP;
227 using Thyra::VectorBase;
228 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
229 "Error, setupInOutArgs_ must be called first!\n");
230
231 RCP<const VectorBase<Scalar> > x_in = inArgs.get_x();
232 double beta = inArgs.get_beta();
233 if (!x_in.get()) {
234 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
235 "\n ERROR: HarmonicOscillatorModel requires x as InArgs.\n");
236 }
237 Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
238 //IKT, FIXME: check that subDim() is the write routine to get local length of a Thyra::ConstDetachedVectorView
239 auto myVecLength = x_in_view.subDim();
240
241 RCP<const VectorBase<Scalar> > x_dot_in = inArgs.get_x_dot();
242 double alpha = inArgs.get_alpha();
243
244 RCP<const VectorBase<Scalar> > x_dotdot_in = inArgs.get_x_dot_dot();
245 double omega = inArgs.get_W_x_dot_dot_coeff();
246
247 //Parse OutArgs
248 RCP<VectorBase<Scalar> > f_out = outArgs.get_f();
249 RCP<VectorBase<Scalar> > g_out = outArgs.get_g(0);
250 const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
251
252 Scalar neg_sign = 1.0;
253 //Explicit ODE
254 if (inArgs.get_x_dot_dot().is_null()) neg_sign = -1.0;
255
256 //Populate residual and Jacobian
257 if (f_out != Teuchos::null) {
258 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
259 for (int i=0; i<myVecLength; i++) {
260 f_out_view[i] = f_;
261 }
262 if (x_dotdot_in != Teuchos::null) {
263 Thyra::ConstDetachedVectorView<Scalar> x_dotdot_in_view( *x_dotdot_in);
264 for (int i=0; i<myVecLength; i++) {
265 f_out_view[i] = x_dotdot_in_view[i] - f_out_view[i];
266 }
267 }
268 if (x_dot_in != Teuchos::null) {
269 Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in);
270 for (int i=0; i<myVecLength; i++) {
271 f_out_view[i] += neg_sign*c_*x_dot_in_view[i];
272 }
273 }
274 if (x_in != Teuchos::null) {
275 for (int i=0; i<myVecLength; i++) {
276 f_out_view[i] += neg_sign*k_*x_in_view[i];
277 }
278 }
279 }
280
281 // Note: W = alpha*df/dxdot + beta*df/dx + omega*df/dxdotdot
282 if (W_out != Teuchos::null) {
283 RCP<Thyra::MultiVectorBase<Scalar> > matrix = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
284 Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
285 if (omega == 0.0) {
286 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
287 "\n ERROR: omega = 0 in HarmonicOscillatorModel!\n");
288 }
289 matrix_view(0,0) = omega;
290 if (x_dot_in != Teuchos::null) {
291 matrix_view(0,0) += neg_sign*c_*alpha;
292 }
293 if (x_in != Teuchos::null) {
294 matrix_view(0,0) += neg_sign*k_*beta;
295 }
296 }
297
298 //Calculated response(s) g
299 //g = mean value of x
300 if (g_out != Teuchos::null) {
301 Thyra::DetachedVectorView<Scalar> g_out_view(*g_out);
302 g_out_view[0] = Thyra::sum(*x_in)/vecLength_;
303 }
304}
305
306template<class Scalar>
307Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
309get_p_space(int l) const
310{
311 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
312 "\n Error! HarmonicOscillatorModel::get_p_space() is not supported!\n");
313 return Teuchos::null;
314}
315
316template<class Scalar>
317Teuchos::RCP<const Teuchos::Array<std::string> >
319get_p_names(int l) const
320{
321 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
322 "\n Error! HarmonicOscillatorModel::get_p_names() is not supported!\n");
323 return Teuchos::null;
324}
325
326template<class Scalar>
327Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
329get_g_space(int j) const
330{
331 TEUCHOS_TEST_FOR_EXCEPTION(j != 0, std::logic_error,
332 "\n Error! HarmonicOscillatorModel::get_g_space() only " <<
333 " supports 1 parameter vector. Supplied index l = " <<
334 j << "\n");
335 return g_space_;
336}
337
338// private
339
340template<class Scalar>
341void
343setupInOutArgs_() const
344{
345 if (isInitialized_) return;
346
347 //Set up InArgs
348 Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
349 inArgs.setModelEvalDescription(this->description());
350 inArgs.set_Np(0);
351 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x);
352 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot);
353 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot_dot);
354 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_t);
355 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_W_x_dot_dot_coeff);
356 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_alpha);
357 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_beta);
358 inArgs_ = inArgs;
359
360 //Set up OutArgs
361 Thyra::ModelEvaluatorBase::OutArgsSetup<Scalar> outArgs;
362 outArgs.setModelEvalDescription(this->description());
363 outArgs.set_Np_Ng(0, numResponses_);
364
365 outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f);
366 outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W_op);
367 //outArgs.setSupports(OUT_ARG_W,true);
368 //IKT, what is the following supposed to do??
369 //outArgs.set_W_properties( DerivativeProperties(
370 // DERIV_LINEARITY_UNKNOWN, DERIV_RANK_FULL, true));
371 outArgs_ = outArgs;
372
373 // Set up nominal values
374 nominalValues_ = inArgs_;
375 nominalValues_.set_t(0.0);
376 nominalValues_.set_x(x_vec_);
377 nominalValues_.set_x_dot(x_dot_vec_);
378 nominalValues_.set_x_dot_dot(x_dot_dot_vec_);
379
380 isInitialized_ = true;
381
382}
383
384template<class Scalar>
385void
387setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
388{
389 using Teuchos::get;
390 using Teuchos::RCP;
391 using Teuchos::ParameterList;
392 RCP<ParameterList> tmpPL = Teuchos::rcp(new ParameterList("HarmonicOscillatorModel"));
393 if (paramList != Teuchos::null) tmpPL = paramList;
394 tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
395 this->setMyParamList(tmpPL);
396 RCP<ParameterList> pl = this->getMyNonconstParamList();
397 c_ = get<Scalar>(*pl,"Damping coeff c");
398 f_ = get<Scalar>(*pl,"Forcing coeff f");
399 k_ = get<Scalar>(*pl,"x coeff k");
400 m_ = get<Scalar>(*pl,"Mass coeff m");
401 if (m_ <= 0.0) {
402 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
403 "Error: invalid value of Mass coeff m = " << m_ <<"! Mass coeff m must be > 0.\n");
404 }
405 if (k_ < 0.0) {
406 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
407 "Error: invalid value of x coeff k = " << k_ <<"! x coeff k must be >= 0.\n");
408 }
409 if ((k_ > 0.0) && (c_ != 0.0)) {
410 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
411 "Error: HarmonicOscillator model only supports x coeff k > 0 when Damping coeff c = 0. You have "
412 << "specified x coeff k = " << k_ << " and Damping coeff c = " << c_ << ".\n");
413 }
414}
415
416template<class Scalar>
417Teuchos::RCP<const Teuchos::ParameterList>
419getValidParameters() const
420{
421 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
422 if (is_null(validPL)) {
423 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
424 validPL = pl;
425 Teuchos::setDoubleParameter(
426 "Damping coeff c", 0.0, "Damping coefficient in model", &*pl);
427 Teuchos::setDoubleParameter(
428 "Forcing coeff f", -1.0, "Forcing coefficient in model", &*pl);
429 Teuchos::setDoubleParameter(
430 "x coeff k", 0.0, "x coefficient in model", &*pl);
431 Teuchos::setDoubleParameter(
432 "Mass coeff m", 1.0, "Mass coefficient in model", &*pl);
433 }
434 return validPL;
435}
436
437} // namespace Tempus_Test
438#endif // TEMPUS_TEST_HARMONIC_OSCILLATOR_MODEL_IMPL_HPP
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > g_space_
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Teuchos::RCP< Thyra::VectorBase< Scalar > > x_dot_vec_
Teuchos::RCP< Thyra::VectorBase< Scalar > > x_vec_
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getExactSolution(double t) const
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
Teuchos::RCP< Teuchos::FancyOStream > out_
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const
Teuchos::RCP< Thyra::VectorBase< Scalar > > x_dot_dot_vec_
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
HarmonicOscillatorModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null, const bool use_accel_IC=false)
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > x_space_