Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_BackwardEuler_ASA_AlgebraicEquations.cpp
Go to the documentation of this file.
1//@HEADER
2// *****************************************************************************
3// Tempus: Time Integration and Sensitivity Analysis Package
4//
5// Copyright 2017 NTESS and the Tempus contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8//@HEADER
9
10// Adjoint sensitivities of a semi-explicit index-1 DAE ("Algebraic Equations
11// Mask" option of IntegratorAdjointSensitivity), for the "Final Time", "Sum
12// Over Times" and "Time Integral" responses.
13
14#include <cmath>
15#include <string>
16#include <vector>
17
18#include "Teuchos_UnitTestHarness.hpp"
19#include "Teuchos_XMLParameterListHelpers.hpp"
20
21#include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
22#include "Thyra_DefaultSpmdVectorSpace.hpp"
23#include "Thyra_DetachedMultiVectorView.hpp"
24#include "Thyra_DetachedVectorView.hpp"
25#include "Thyra_MultiVectorStdOps.hpp"
26#include "Thyra_StateFuncModelEvaluatorBase.hpp"
27#include "Thyra_VectorStdOps.hpp"
28
29#include "Tempus_DAEStructure.hpp"
30#include "Tempus_IntegratorAdjointSensitivity.hpp"
31
32#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
33
34namespace Tempus_Test {
35
36using Teuchos::ParameterList;
37using Teuchos::RCP;
38using Teuchos::rcp;
39
53class AlgebraicDAEModel : public Thyra::StateFuncModelEvaluatorBase<double> {
54 public:
55 typedef Thyra::ModelEvaluatorBase MEB;
56
57 AlgebraicDAEModel(const double kappa = 0.5, const double h0 = 1.0,
58 const double p0 = 0.3, const double q0 = 0.1)
59 : kappa_(kappa), h0_(h0), p0_(p0), q0_(q0)
60 {
61 space_ = Thyra::defaultSpmdVectorSpace<double>(2);
62 p_names_.push_back("p");
63 p_names_.push_back("q");
64 g_names_.push_back("g = x");
65
66 // The descriptions must be this->description(), which is checked by
67 // Thyra::ModelEvaluatorDefaultBase in debug builds
68 MEB::InArgsSetup<double> inArgs;
69 inArgs.setModelEvalDescription(this->description());
70 inArgs.setSupports(MEB::IN_ARG_t);
71 inArgs.setSupports(MEB::IN_ARG_x);
72 inArgs.setSupports(MEB::IN_ARG_x_dot);
73 inArgs.setSupports(MEB::IN_ARG_alpha);
74 inArgs.setSupports(MEB::IN_ARG_beta);
75 inArgs.set_Np(1);
76 inArgs_ = inArgs;
77
78 MEB::OutArgsSetup<double> outArgs;
79 outArgs.setModelEvalDescription(this->description());
80 outArgs.setSupports(MEB::OUT_ARG_f);
81 outArgs.setSupports(MEB::OUT_ARG_W_op);
82 outArgs.set_Np_Ng(1, 1);
83 outArgs.setSupports(MEB::OUT_ARG_DfDp, 0, MEB::DERIV_MV_JACOBIAN_FORM);
84 outArgs.setSupports(MEB::OUT_ARG_DgDx, 0, MEB::DERIV_MV_GRADIENT_FORM);
85 outArgs.setSupports(MEB::OUT_ARG_DgDp, 0, 0, MEB::DERIV_MV_JACOBIAN_FORM);
86 outArgs_ = outArgs;
87
88 nominal_ = inArgs_;
89 nominal_.set_t(0.0);
90 RCP<Thyra::VectorBase<double> > x = Thyra::createMember(space_);
91 RCP<Thyra::VectorBase<double> > x_dot = Thyra::createMember(space_);
92 RCP<Thyra::VectorBase<double> > p = Thyra::createMember(space_);
93 {
94 Thyra::DetachedVectorView<double> xv(*x), xdv(*x_dot), pv(*p);
95 xv[0] = kappa_ * h0_ + p0_;
96 xv[1] = h0_;
97 xdv[0] = 0.0;
98 xdv[1] = a() * h0_ + p0_ + 2.0 * q0_;
99 pv[0] = p0_;
100 pv[1] = q0_;
101 }
102 nominal_.set_x(x);
103 nominal_.set_x_dot(x_dot);
104 nominal_.set_p(0, p);
105 }
106
107 double a() const { return kappa_ - 1.0; }
108 double kappa() const { return kappa_; }
109
111 RCP<const Thyra::VectorBase<double> > getInitialCondition() const
112 {
113 return nominal_.get_x();
114 }
115
117 RCP<Thyra::MultiVectorBase<double> > getInitialSensitivity() const
118 {
119 RCP<Thyra::MultiVectorBase<double> > dxdp = Thyra::createMembers(space_, 2);
120 Thyra::assign(dxdp.ptr(), 0.0);
121 Thyra::DetachedMultiVectorView<double> v(*dxdp);
122 v(0, 0) = 1.0;
123 return dxdp;
124 }
125
127 void exactSensitivity(const double t, double s[2][2]) const
128 {
129 const double st = (std::exp(a() * t) - 1.0) / a();
130 s[0][0] = kappa_ * st + 1.0;
131 s[0][1] = 2.0 * kappa_ * st;
132 s[1][0] = st;
133 s[1][1] = 2.0 * st;
134 }
135
137 void exactIntegratedSensitivity(const double t, double s[2][2]) const
138 {
139 const double St = ((std::exp(a() * t) - 1.0) / a() - t) / a();
140 s[0][0] = kappa_ * St + t;
141 s[0][1] = 2.0 * kappa_ * St;
142 s[1][0] = St;
143 s[1][1] = 2.0 * St;
144 }
145
146 RCP<const Thyra::VectorSpaceBase<double> > get_x_space() const
147 {
148 return space_;
149 }
150 RCP<const Thyra::VectorSpaceBase<double> > get_f_space() const
151 {
152 return space_;
153 }
154 RCP<const Thyra::VectorSpaceBase<double> > get_p_space(int l) const
155 {
156 TEUCHOS_ASSERT_EQUALITY(l, 0);
157 return space_;
158 }
159 RCP<const Teuchos::Array<std::string> > get_p_names(int l) const
160 {
161 TEUCHOS_ASSERT_EQUALITY(l, 0);
162 return Teuchos::rcpFromRef(p_names_);
163 }
164 RCP<const Thyra::VectorSpaceBase<double> > get_g_space(int j) const
165 {
166 TEUCHOS_ASSERT_EQUALITY(j, 0);
167 return space_;
168 }
169 Teuchos::ArrayView<const std::string> get_g_names(int j) const
170 {
171 TEUCHOS_ASSERT_EQUALITY(j, 0);
172 return g_names_();
173 }
174 MEB::InArgs<double> getNominalValues() const { return nominal_; }
175 MEB::InArgs<double> createInArgs() const { return inArgs_; }
176 RCP<Thyra::LinearOpBase<double> > create_W_op() const
177 {
178 return Thyra::createMembers(space_, 2);
179 }
180 RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > get_W_factory()
181 const
182 {
183 return Thyra::defaultSerialDenseLinearOpWithSolveFactory<double>();
184 }
185
186 private:
187 MEB::OutArgs<double> createOutArgsImpl() const { return outArgs_; }
188
189 void evalModelImpl(const MEB::InArgs<double>& inArgs,
190 const MEB::OutArgs<double>& outArgs) const
191 {
192 Thyra::ConstDetachedVectorView<double> x(*inArgs.get_x().assert_not_null());
193 RCP<const Thyra::VectorBase<double> > p_in = inArgs.get_p(0);
194 if (p_in == Teuchos::null) p_in = nominal_.get_p(0);
195 Thyra::ConstDetachedVectorView<double> p(*p_in);
196 double h_dot = 0.0;
197 if (inArgs.get_x_dot() != Teuchos::null)
198 h_dot = Thyra::get_ele(*inArgs.get_x_dot(), 1);
199
200 RCP<Thyra::VectorBase<double> > f = outArgs.get_f();
201 if (f != Teuchos::null) {
202 Thyra::DetachedVectorView<double> fv(*f);
203 fv[0] = x[0] - kappa_ * x[1] - p[0];
204 fv[1] = h_dot + x[1] - x[0] - 2.0 * p[1];
205 }
206
207 // W = alpha df/dx_dot + beta df/dx, with df/dx_dot = [0 0; 0 1]
208 RCP<Thyra::LinearOpBase<double> > W_op = outArgs.get_W_op();
209 if (W_op != Teuchos::null) {
210 const double alpha = inArgs.get_alpha();
211 const double beta = inArgs.get_beta();
212 RCP<Thyra::MultiVectorBase<double> > W =
213 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<double> >(W_op, true);
214 Thyra::DetachedMultiVectorView<double> Wv(*W);
215 Wv(0, 0) = beta;
216 Wv(0, 1) = -kappa_ * beta;
217 Wv(1, 0) = -beta;
218 Wv(1, 1) = alpha + beta;
219 }
220
221 RCP<Thyra::MultiVectorBase<double> > dfdp = outArgs.get_DfDp(0).getMultiVector();
222 if (dfdp != Teuchos::null) {
223 Thyra::DetachedMultiVectorView<double> v(*dfdp);
224 v(0, 0) = -1.0;
225 v(0, 1) = 0.0;
226 v(1, 0) = 0.0;
227 v(1, 1) = -2.0;
228 }
229
230 RCP<Thyra::VectorBase<double> > g = outArgs.get_g(0);
231 if (g != Teuchos::null) Thyra::assign(g.ptr(), *inArgs.get_x());
232
233 RCP<Thyra::MultiVectorBase<double> > dgdx = outArgs.get_DgDx(0).getMultiVector();
234 if (dgdx != Teuchos::null) {
235 Thyra::DetachedMultiVectorView<double> v(*dgdx);
236 v(0, 0) = 1.0;
237 v(0, 1) = 0.0;
238 v(1, 0) = 0.0;
239 v(1, 1) = 1.0;
240 }
241
242 RCP<Thyra::MultiVectorBase<double> > dgdp = outArgs.get_DgDp(0, 0).getMultiVector();
243 if (dgdp != Teuchos::null) Thyra::assign(dgdp.ptr(), 0.0);
244 }
245
246 double kappa_, h0_, p0_, q0_;
247 RCP<const Thyra::VectorSpaceBase<double> > space_;
248 Teuchos::Array<std::string> p_names_;
249 Teuchos::Array<std::string> g_names_;
250 MEB::InArgs<double> inArgs_;
251 MEB::OutArgs<double> outArgs_;
252 MEB::InArgs<double> nominal_;
253};
254
256RCP<ParameterList> getTempusParameters(const double dt)
257{
258 RCP<ParameterList> pl = Teuchos::getParametersFromXmlString(
259 "<ParameterList name=\"Tempus\">"
260 " <Parameter name=\"Integrator Name\" type=\"string\""
261 " value=\"Default Integrator\"/>"
262 " <ParameterList name=\"Default Integrator\">"
263 " <Parameter name=\"Integrator Type\" type=\"string\""
264 " value=\"Integrator Basic\"/>"
265 " <Parameter name=\"Screen Output Index Interval\" type=\"int\""
266 " value=\"100000\"/>"
267 " <Parameter name=\"Stepper Name\" type=\"string\""
268 " value=\"Default Stepper\"/>"
269 " <ParameterList name=\"Solution History\">"
270 " <Parameter name=\"Storage Type\" type=\"string\""
271 " value=\"Unlimited\"/>"
272 " <Parameter name=\"Storage Limit\" type=\"int\" value=\"20\"/>"
273 " <ParameterList name=\"Interpolator\">"
274 " <Parameter name=\"Interpolator Type\" type=\"string\""
275 " value=\"Lagrange\"/>"
276 " <Parameter name=\"Order\" type=\"int\" value=\"0\"/>"
277 " </ParameterList>"
278 " </ParameterList>"
279 " <ParameterList name=\"Time Step Control\">"
280 " <Parameter name=\"Initial Time\" type=\"double\" value=\"0.0\"/>"
281 " <Parameter name=\"Final Time\" type=\"double\" value=\"1.0\"/>"
282 " <Parameter name=\"Final Time Index\" type=\"int\""
283 " value=\"100000\"/>"
284 " <Parameter name=\"Minimum Time Step\" type=\"double\""
285 " value=\"1.0e-6\"/>"
286 " <Parameter name=\"Maximum Time Step\" type=\"double\""
287 " value=\"1.0\"/>"
288 " </ParameterList>"
289 " </ParameterList>"
290 " <ParameterList name=\"Default Stepper\">"
291 " <Parameter name=\"Stepper Type\" type=\"string\""
292 " value=\"Backward Euler\"/>"
293 " <Parameter name=\"Use FSAL\" type=\"bool\" value=\"false\"/>"
294 " <Parameter name=\"Initial Condition Consistency\" type=\"string\""
295 " value=\"None\"/>"
296 " <Parameter name=\"Initial Condition Consistency Check\""
297 " type=\"bool\" value=\"false\"/>"
298 " <Parameter name=\"Solver Name\" type=\"string\""
299 " value=\"Default Solver\"/>"
300 " <Parameter name=\"Zero Initial Guess\" type=\"bool\" value=\"0\"/>"
301 " <ParameterList name=\"Default Solver\">"
302 " <ParameterList name=\"NOX\">"
303 " <ParameterList name=\"Direction\">"
304 " <Parameter name=\"Method\" type=\"string\" value=\"Newton\"/>"
305 " </ParameterList>"
306 " <ParameterList name=\"Line Search\">"
307 " <Parameter name=\"Method\" type=\"string\""
308 " value=\"Full Step\"/>"
309 " </ParameterList>"
310 " <Parameter name=\"Nonlinear Solver\" type=\"string\""
311 " value=\"Line Search Based\"/>"
312 " <ParameterList name=\"Status Tests\">"
313 " <Parameter name=\"Test Type\" type=\"string\" value=\"Combo\"/>"
314 " <Parameter name=\"Combo Type\" type=\"string\" value=\"OR\"/>"
315 " <Parameter name=\"Number of Tests\" type=\"int\" value=\"2\"/>"
316 " <ParameterList name=\"Test 0\">"
317 " <Parameter name=\"Test Type\" type=\"string\""
318 " value=\"NormF\"/>"
319 " <Parameter name=\"Tolerance\" type=\"double\""
320 " value=\"1.0e-12\"/>"
321 " </ParameterList>"
322 " <ParameterList name=\"Test 1\">"
323 " <Parameter name=\"Test Type\" type=\"string\""
324 " value=\"MaxIters\"/>"
325 " <Parameter name=\"Maximum Iterations\" type=\"int\""
326 " value=\"10\"/>"
327 " </ParameterList>"
328 " </ParameterList>"
329 " </ParameterList>"
330 " </ParameterList>"
331 " </ParameterList>"
332 "</ParameterList>");
333 pl->sublist("Default Integrator")
334 .sublist("Time Step Control")
335 .set("Initial Time Step", dt);
336 return pl;
337}
338
340RCP<const Thyra::VectorBase<double> > createMask(
341 const AlgebraicDAEModel& model, const double m0, const double m1)
342{
343 RCP<Thyra::VectorBase<double> > mask = Thyra::createMember(model.get_f_space());
344 Thyra::DetachedVectorView<double> v(*mask);
345 v[0] = m0;
346 v[1] = m1;
347 return mask;
348}
349
351RCP<Thyra::MultiVectorBase<double> > runAdjoint(
352 const RCP<AlgebraicDAEModel>& model, const std::string& responseType,
353 const std::string& responseTimes, const double dt,
354 const RCP<const Thyra::VectorBase<double> >& mask, bool& status)
355{
356 RCP<ParameterList> pl = getTempusParameters(dt);
357 ParameterList& sens_pl = pl->sublist("Sensitivities");
358 sens_pl.set("Mass Matrix Is Identity", false);
359 sens_pl.set("Response Type", responseType);
360 if (!responseTimes.empty()) sens_pl.set("Response Time List", responseTimes);
361 sens_pl.set(Tempus::DAEStructure<double>::algebraicMaskName(), mask);
362 sens_pl.sublist("Algebraic Jump Solver").set("Relative Tolerance", 1.0e-12);
363
364 RCP<Tempus::IntegratorAdjointSensitivity<double> > integrator =
365 Tempus::createIntegratorAdjointSensitivity<double>(pl, model);
366 integrator->initializeSolutionHistory(
367 0.0, model->getInitialCondition(), Teuchos::null, Teuchos::null,
368 model->getInitialSensitivity(), Teuchos::null, Teuchos::null);
369 status = integrator->advanceTime();
370
371 // getDgDp() is in gradient form: (p, g)
372 RCP<const Thyra::MultiVectorBase<double> > DgDp = integrator->getDgDp();
373 RCP<Thyra::MultiVectorBase<double> > dgdp = Thyra::createMembers(model->get_g_space(0), 2);
374 Thyra::ConstDetachedMultiVectorView<double> grad(*DgDp);
375 Thyra::DetachedMultiVectorView<double> jac(*dgdp);
376 for (int i = 0; i < 2; ++i)
377 for (int j = 0; j < 2; ++j) jac(i, j) = grad(j, i);
378 return dgdp;
379}
380
382void checkConvergence(const std::string& responseType,
383 Teuchos::FancyOStream& out, bool& success)
384{
385 RCP<AlgebraicDAEModel> model = rcp(new AlgebraicDAEModel());
386 const std::vector<double> times = {0.25, 0.5, 1.0};
387 const std::string timeList = "0.25, 0.5, 1.0";
388
389 // Exact sensitivity of the response
390 double exact[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
391 if (responseType == "Final Time")
392 model->exactSensitivity(1.0, exact);
393 else if (responseType == "Time Integral")
394 model->exactIntegratedSensitivity(1.0, exact);
395 else {
396 for (const double t : times) {
397 double s[2][2];
398 model->exactSensitivity(t, s);
399 for (int i = 0; i < 2; ++i)
400 for (int j = 0; j < 2; ++j) exact[i][j] += s[i][j];
401 }
402 }
403
404 std::vector<double> stepSize, errorNorm;
405 for (int n = 0; n < 5; ++n) {
406 const double dt = 1.0 / (16 << n); // exact in binary
407 bool status = false;
408 RCP<Thyra::MultiVectorBase<double> > dgdp =
409 runAdjoint(model, responseType,
410 responseType == "Sum Over Times" ? timeList : "", dt,
411 createMask(*model, 1.0, 0.0), status);
412 TEST_ASSERT(status);
413
414 // Backward Euler sensitivity s_k of h at t_k = k dt (exact for the
415 // discrete solution, for comparison only)
416 double s_BE = 0.0, sum_BE = 0.0;
417 const int N = 16 << n;
418 for (int k = 1; k <= N; ++k) {
419 s_BE = (s_BE + dt) / (1.0 - model->a() * dt);
420 for (const double t : times)
421 if (std::abs(k * dt - t) < 1.0e-12) sum_BE += s_BE;
422 }
423 const double s_disc = responseType == "Final Time" ? s_BE : sum_BE;
424 const double ns = responseType == "Final Time" ? 1.0 : 3.0;
425
426 Thyra::ConstDetachedMultiVectorView<double> v(*dgdp);
427 double err2 = 0.0;
428 for (int i = 0; i < 2; ++i)
429 for (int j = 0; j < 2; ++j)
430 err2 += (v(i, j) - exact[i][j]) * (v(i, j) - exact[i][j]);
431 stepSize.push_back(dt);
432 errorNorm.push_back(std::sqrt(err2));
433 out << " " << responseType << ", dt = " << dt << ": du/dp = " << v(0, 0)
434 << ", du/dq = " << v(0, 1) << ", dh/dp = " << v(1, 0)
435 << ", dh/dq = " << v(1, 1) << ", |error| = " << std::sqrt(err2);
436 if (responseType != "Time Integral")
437 out << ", du/dp - (discrete BE) = "
438 << v(0, 0) - (model->kappa() * s_disc + ns);
439 out << std::endl;
440 }
441
442 const double slope = computeLinearRegressionLogLog<double>(stepSize, errorNorm);
443 out << " " << responseType << ": order = " << slope << std::endl;
444 TEST_FLOATING_EQUALITY(slope, 1.0, 0.05);
445 TEST_COMPARE(errorNorm.back(), <, 1.0e-2);
446}
447
448// ************************************************************
449TEUCHOS_UNIT_TEST(BackwardEuler, ASA_AlgebraicEquations_FinalTime)
450{
451 checkConvergence("Final Time", out, success);
452}
453
454// ************************************************************
455TEUCHOS_UNIT_TEST(BackwardEuler, ASA_AlgebraicEquations_SumOverTimes)
456{
457 checkConvergence("Sum Over Times", out, success);
458}
459
460// ************************************************************
461TEUCHOS_UNIT_TEST(BackwardEuler, ASA_AlgebraicEquations_TimeIntegral)
462{
463 checkConvergence("Time Integral", out, success);
464}
465
466// ************************************************************
467// The algebraic component of the jump: for g = u at the final time, the
468// sensitivity to p jumps by 1 (the ODE formula would need df/dx_dot^{-T},
469// which does not exist).
470TEUCHOS_UNIT_TEST(BackwardEuler, ASA_AlgebraicEquations_AlgebraicJump)
471{
472 RCP<AlgebraicDAEModel> model = rcp(new AlgebraicDAEModel());
473 bool status = false;
474 const double dt = 1.0 / 64;
475 RCP<Thyra::MultiVectorBase<double> > dgdp = runAdjoint(
476 model, "Final Time", "", dt, createMask(*model, 1.0, 0.0), status);
477 TEST_ASSERT(status);
478 Thyra::ConstDetachedMultiVectorView<double> v(*dgdp);
479 // du/dp = kappa dh/dp + 1 and du/dq = kappa dh/dq hold exactly for the
480 // discrete solution
481 TEST_FLOATING_EQUALITY(v(0, 0), model->kappa() * v(1, 0) + 1.0, 1.0e-10);
482 TEST_FLOATING_EQUALITY(v(0, 1), model->kappa() * v(1, 1), 1.0e-10);
483}
484
485// ************************************************************
486TEUCHOS_UNIT_TEST(BackwardEuler, ASA_AlgebraicEquations_InvalidMask)
487{
488 RCP<AlgebraicDAEModel> model = rcp(new AlgebraicDAEModel());
489 bool status = false;
490
491 // Entries must be 0 or 1
492 TEST_THROW(runAdjoint(model, "Final Time", "", 0.1,
493 createMask(*model, 0.5, 0.0), status),
494 std::invalid_argument);
495
496 // The differential equation depends on x_dot and cannot be flagged
497 bool massMatrixError = false;
498 try {
499 runAdjoint(model, "Final Time", "", 0.1, createMask(*model, 1.0, 1.0),
500 status);
501 }
502 catch (const std::logic_error& e) {
503 massMatrixError = std::string(e.what()).find(
504 "the rows of the mass matrix df/dx_dot") !=
505 std::string::npos;
506 if (!massMatrixError) out << "Unexpected exception: " << e.what() << std::endl;
507 }
508 TEST_ASSERT(massMatrixError);
509}
510
511// ************************************************************
512TEUCHOS_UNIT_TEST(BackwardEuler, ASA_AlgebraicEquations_DAEStructure)
513{
514 RCP<AlgebraicDAEModel> model = rcp(new AlgebraicDAEModel());
515 Tempus::DAEStructure<double> dae(createMask(*model, 1.0, 0.0));
516 TEST_EQUALITY(Thyra::get_ele(*dae.getDifferentialMask(), 0), 0.0);
517 TEST_EQUALITY(Thyra::get_ele(*dae.getDifferentialMask(), 1), 1.0);
518
519 RCP<Thyra::MultiVectorBase<double> > v = Thyra::createMembers(model->get_f_space(), 2);
520 Thyra::assign(v.ptr(), 3.0);
521 dae.zeroAlgebraicComponents(*v);
522 {
523 Thyra::ConstDetachedMultiVectorView<double> w(*v);
524 TEST_EQUALITY(w(0, 0), 0.0);
525 TEST_EQUALITY(w(0, 1), 0.0);
526 TEST_EQUALITY(w(1, 0), 3.0);
527 TEST_EQUALITY(w(1, 1), 3.0);
528 }
529 Thyra::assign(v.ptr(), 3.0);
530 dae.zeroDifferentialComponents(*v);
531 {
532 Thyra::ConstDetachedMultiVectorView<double> w(*v);
533 TEST_EQUALITY(w(0, 0), 3.0);
534 TEST_EQUALITY(w(1, 1), 0.0);
535 }
536
537 // Absent or null mask: no DAE structure
538 ParameterList pl;
539 TEST_ASSERT(Tempus::getDAEStructure<double>(pl) == Teuchos::null);
540 pl.set(Tempus::DAEStructure<double>::algebraicMaskName(),
541 createMask(*model, 1.0, 0.0));
542 TEST_ASSERT(Tempus::getDAEStructure<double>(pl) != Teuchos::null);
543}
544
545} // namespace Tempus_Test
Tempus_ConvergenceTestUtils.hpp
Tempus_DAEStructure.hpp
Tempus::DAEStructure
Structure of an index-1 differential-algebraic system.
Definition Tempus_DAEStructure.hpp:55
Tempus::DAEStructure::zeroDifferentialComponents
void zeroDifferentialComponents(MultiVector &v) const
v <- P_a v: zero the differential components of each column of v
Definition Tempus_DAEStructure.hpp:114
Tempus::DAEStructure::getDifferentialMask
Teuchos::RCP< const Vector > getDifferentialMask() const
Mask of the other equations (1 - mask of the algebraic equations).
Definition Tempus_DAEStructure.hpp:101
Tempus::DAEStructure::algebraicMaskName
static const std::string & algebraicMaskName()
Name of the parameter giving the mask of the algebraic equations.
Definition Tempus_DAEStructure.hpp:62
Tempus::DAEStructure::zeroAlgebraicComponents
void zeroAlgebraicComponents(MultiVector &v) const
v <- P_d v: zero the algebraic components of each column of v
Definition Tempus_DAEStructure.hpp:107
Tempus_Test
Definition Tempus_BackwardEuler_ASA.cpp:34
Tempus_Test::getTempusParameters
RCP< ParameterList > getTempusParameters(const double dt)
Tempus parameters: Backward Euler with a constant time step dt on [0, 1].
Definition Tempus_BackwardEuler_ASA_AlgebraicEquations.cpp:256
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition Tempus_BackwardEuler_ASA.cpp:47
Tempus_Test::createMask
RCP< const Thyra::VectorBase< double > > createMask(const AlgebraicDAEModel &model, const double m0, const double m1)
Mask of the algebraic equations.
Definition Tempus_BackwardEuler_ASA_AlgebraicEquations.cpp:340
Tempus_Test::runAdjoint
RCP< Thyra::MultiVectorBase< double > > runAdjoint(const RCP< AlgebraicDAEModel > &model, const std::string &responseType, const std::string &responseTimes, const double dt, const RCP< const Thyra::VectorBase< double > > &mask, bool &status)
Run the adjoint sensitivity analysis; returns dg/dp (Jacobian form).
Definition Tempus_BackwardEuler_ASA_AlgebraicEquations.cpp:351
Tempus_Test::checkConvergence
void checkConvergence(const std::string &responseType, Teuchos::FancyOStream &out, bool &success)
Convergence of the adjoint sensitivities to the exact ones (order 1).
Definition Tempus_BackwardEuler_ASA_AlgebraicEquations.cpp:382
Tempus_Test::computeLinearRegressionLogLog
Scalar computeLinearRegressionLogLog(std::vector< Scalar > &x, std::vector< Scalar > &y)
Definition Tempus_ConvergenceTestUtils.hpp:159
Tempus::getDAEStructure
Teuchos::RCP< const DAEStructure< Scalar > > getDAEStructure(const Teuchos::ParameterList &pl)
The DAE structure given by the "Algebraic Equations Mask" parameter of pl, or null if the parameter i...
Definition Tempus_DAEStructure.hpp:171
Tempus::createIntegratorAdjointSensitivity
Teuchos::RCP< IntegratorAdjointSensitivity< Scalar > > createIntegratorAdjointSensitivity()
Nonmember constructor.
Definition Tempus_IntegratorAdjointSensitivity_impl.hpp:900
Generated by doxygen 1.15.0