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