Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_Test_BackwardEuler_OptInterface.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 #include "Teuchos_UnitTestHarness.hpp"
11 #include "Teuchos_XMLParameterListHelpers.hpp"
12 #include "Teuchos_TimeMonitor.hpp"
13 #include "Teuchos_DefaultComm.hpp"
14 
15 #include "Thyra_MultiVectorStdOps.hpp"
16 
17 #include "Tempus_config.hpp"
18 #include "Tempus_IntegratorBasic.hpp"
19 #include "Tempus_StepperBackwardEuler.hpp"
20 
21 #include "../TestModels/SinCosModel.hpp"
22 #include "../TestModels/VanDerPolModel.hpp"
23 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
24 
25 #include <vector>
26 #include <fstream>
27 #include <sstream>
28 #include <limits>
29 
30 namespace Tempus_Test {
31 
32 using Teuchos::getParametersFromXmlFile;
33 using Teuchos::ParameterList;
34 using Teuchos::RCP;
35 using Teuchos::rcp;
36 using Teuchos::rcp_const_cast;
37 using Teuchos::sublist;
38 
39 using Tempus::IntegratorBasic;
40 using Tempus::SolutionHistory;
41 using Tempus::SolutionState;
42 
43 // ************************************************************
44 // ************************************************************
45 TEUCHOS_UNIT_TEST(BackwardEuler, OptInterface)
46 {
47  // Read params from .xml file
48  RCP<ParameterList> pList =
49  getParametersFromXmlFile("Tempus_BackwardEuler_SinCos.xml");
50 
51  // Setup the SinCosModel
52  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
53  auto model = rcp(new SinCosModel<double>(scm_pl));
54 
55  // Setup the Integrator
56  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
57  RCP<Tempus::IntegratorBasic<double> > integrator =
58  Tempus::createIntegratorBasic<double>(pl, model);
59 
60  // Integrate to timeMax
61  bool integratorStatus = integrator->advanceTime();
62  TEST_ASSERT(integratorStatus);
63 
64  // Get solution history
65  RCP<const SolutionHistory<double> > solutionHistory =
66  integrator->getSolutionHistory();
67 
68  // Get the stepper and cast to optimization interface
69  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
70  RCP<Tempus::StepperOptimizationInterface<double> > opt_stepper =
71  Teuchos::rcp_dynamic_cast<Tempus::StepperOptimizationInterface<double> >(
72  stepper, true);
73 
74  // Check stencil length
75  TEST_EQUALITY(opt_stepper->stencilLength(), 2);
76 
77  // Create needed vectors/multivectors
78  Teuchos::Array<RCP<const Thyra::VectorBase<double> > > x(2);
79  Teuchos::Array<double> t(2);
80  RCP<const Thyra::VectorBase<double> > p = model->getNominalValues().get_p(0);
81  RCP<Thyra::VectorBase<double> > x_dot =
82  Thyra::createMember(model->get_x_space());
83  RCP<Thyra::VectorBase<double> > f = Thyra::createMember(model->get_f_space());
84  RCP<Thyra::VectorBase<double> > f2 =
85  Thyra::createMember(model->get_f_space());
86  RCP<Thyra::LinearOpBase<double> > dfdx = model->create_W_op();
87  RCP<Thyra::LinearOpBase<double> > dfdx2 = model->create_W_op();
88  RCP<Thyra::MultiVectorBase<double> > dfdx_mv =
89  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<double> >(dfdx, true);
90  RCP<Thyra::MultiVectorBase<double> > dfdx_mv2 =
91  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<double> >(dfdx2, true);
92  const int num_p = p->range()->dim();
93  RCP<Thyra::MultiVectorBase<double> > dfdp =
94  Thyra::createMembers(model->get_f_space(), num_p);
95  RCP<Thyra::MultiVectorBase<double> > dfdp2 =
96  Thyra::createMembers(model->get_f_space(), num_p);
97  RCP<Thyra::LinearOpWithSolveBase<double> > W = model->create_W();
98  RCP<Thyra::LinearOpWithSolveBase<double> > W2 = model->create_W();
99  RCP<Thyra::MultiVectorBase<double> > tmp =
100  Thyra::createMembers(model->get_x_space(), num_p);
101  RCP<Thyra::MultiVectorBase<double> > tmp2 =
102  Thyra::createMembers(model->get_x_space(), num_p);
103  std::vector<double> nrms(num_p);
104  double err;
105 
106  // Loop over states, checking residuals and derivatives
107  const int n = solutionHistory->getNumStates();
108  for (int i = 1; i < n; ++i) {
109  RCP<const SolutionState<double> > state = (*solutionHistory)[i];
110  RCP<const SolutionState<double> > prev_state = (*solutionHistory)[i - 1];
111 
112  // Fill x, t stencils
113  x[0] = state->getX();
114  x[1] = prev_state->getX();
115  t[0] = state->getTime();
116  t[1] = prev_state->getTime();
117 
118  // Compute x_dot
119  const double dt = t[0] - t[1];
120  Thyra::V_StVpStV(x_dot.ptr(), 1.0 / dt, *(x[0]), -1.0 / dt, *(x[1]));
121 
122  // Create model inargs
123  typedef Thyra::ModelEvaluatorBase MEB;
124  MEB::InArgs<double> in_args = model->createInArgs();
125  MEB::OutArgs<double> out_args = model->createOutArgs();
126  in_args.set_x(x[0]);
127  in_args.set_x_dot(x_dot);
128  in_args.set_t(t[0]);
129  in_args.set_p(0, p);
130 
131  const double tol = 1.0e-14;
132 
133  // Check residual
134  opt_stepper->computeStepResidual(*f, x, t, *p, 0);
135  out_args.set_f(f2);
136  model->evalModel(in_args, out_args);
137  out_args.set_f(Teuchos::null);
138  Thyra::V_VmV(f.ptr(), *f, *f2);
139  err = Thyra::norm(*f);
140  TEST_FLOATING_EQUALITY(err, 0.0, tol);
141 
142  // Check df/dx_n
143  // df/dx_n = df/dx_dot * dx_dot/dx_n + df/dx_n = 1/dt*df/dx_dot + df/dx_n
144  opt_stepper->computeStepJacobian(*dfdx, x, t, *p, 0, 0);
145  out_args.set_W_op(dfdx2);
146  in_args.set_alpha(1.0 / dt);
147  in_args.set_beta(1.0);
148  model->evalModel(in_args, out_args);
149  out_args.set_W_op(Teuchos::null);
150  Thyra::V_VmV(dfdx_mv.ptr(), *dfdx_mv, *dfdx_mv2);
151  Thyra::norms(*dfdx_mv, Teuchos::arrayViewFromVector(nrms));
152  err = 0.0;
153  for (auto nrm : nrms) err += nrm;
154  TEST_FLOATING_EQUALITY(err, 0.0, tol);
155 
156  // Check df/dx_{n-1}
157  // df/dx_{n-1} = df/dx_dot * dx_dot/dx_{n-1} = -1/dt*df/dx_dot
158  opt_stepper->computeStepJacobian(*dfdx, x, t, *p, 0, 1);
159  out_args.set_W_op(dfdx2);
160  in_args.set_alpha(-1.0 / dt);
161  in_args.set_beta(0.0);
162  model->evalModel(in_args, out_args);
163  out_args.set_W_op(Teuchos::null);
164  Thyra::V_VmV(dfdx_mv.ptr(), *dfdx_mv, *dfdx_mv2);
165  Thyra::norms(*dfdx_mv, Teuchos::arrayViewFromVector(nrms));
166  err = 0.0;
167  for (auto nrm : nrms) err += nrm;
168  TEST_FLOATING_EQUALITY(err, 0.0, tol);
169 
170  // Check df/dp
171  opt_stepper->computeStepParamDeriv(*dfdp, x, t, *p, 0);
172  out_args.set_DfDp(
173  0, MEB::Derivative<double>(dfdp2, MEB::DERIV_MV_JACOBIAN_FORM));
174  model->evalModel(in_args, out_args);
175  out_args.set_DfDp(0, MEB::Derivative<double>());
176  Thyra::V_VmV(dfdp.ptr(), *dfdp, *dfdp2);
177  Thyra::norms(*dfdp, Teuchos::arrayViewFromVector(nrms));
178  err = 0.0;
179  for (auto nrm : nrms) err += nrm;
180  TEST_FLOATING_EQUALITY(err, 0.0, tol);
181 
182  // Check W
183  opt_stepper->computeStepSolver(*W, x, t, *p, 0);
184  out_args.set_W(W2);
185  in_args.set_alpha(1.0 / dt);
186  in_args.set_beta(1.0);
187  model->evalModel(in_args, out_args);
188  out_args.set_W(Teuchos::null);
189  // note: dfdp overwritten above so dfdp != dfdp2
190  Thyra::solve(*W, Thyra::NOTRANS, *dfdp2, tmp.ptr());
191  Thyra::solve(*W2, Thyra::NOTRANS, *dfdp2, tmp2.ptr());
192  Thyra::V_VmV(tmp.ptr(), *tmp, *tmp2);
193  Thyra::norms(*tmp, Teuchos::arrayViewFromVector(nrms));
194  err = 0.0;
195  for (auto nrm : nrms) err += nrm;
196  TEST_FLOATING_EQUALITY(err, 0.0, tol);
197  }
198 
199  Teuchos::TimeMonitor::summarize();
200 }
201 
202 } // namespace Tempus_Test
NOTRANS
TEST_FLOATING_EQUALITY
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
TEST_ASSERT
#define TEST_ASSERT(v1)
Tempus_Test::SinCosModel
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
Definition: SinCosModel_decl.hpp:109
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition: Tempus_BackwardEuler_ASA.cpp:45
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::TimeMonitor::summarize
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
Teuchos::RCP::ptr
Ptr< T > ptr() const
Teuchos_UnitTestHarness.hpp
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:28
Teuchos_XMLParameterListHelpers.hpp
Tempus::StepperOptimizationInterface
Stepper interface to support full-space optimization.
Definition: Tempus_StepperOptimizationInterface.hpp:41
TEST_EQUALITY
#define TEST_EQUALITY(v1, v2)
Teuchos_TimeMonitor.hpp
Tempus::SolutionState
Solution state for integrators and steppers.
Definition: Tempus_SolutionState_decl.hpp:137
Generated on Thu Dec 26 2024 09:25:16 for Tempus by   doxygen 1.8.5