Tempus  Version of the Day
Time Integration
Tempus_ExplicitRK_FSA.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 #include "Teuchos_UnitTestHarness.hpp"
10 #include "Teuchos_XMLParameterListHelpers.hpp"
11 #include "Teuchos_TimeMonitor.hpp"
12 #include "Teuchos_DefaultComm.hpp"
13 
14 #include "Thyra_VectorStdOps.hpp"
15 #include "Thyra_MultiVectorStdOps.hpp"
16 
17 #include "Tempus_IntegratorBasic.hpp"
18 #include "Tempus_IntegratorForwardSensitivity.hpp"
19 
20 #include "Thyra_DefaultMultiVectorProductVector.hpp"
21 #include "Thyra_DefaultProductVector.hpp"
22 
23 #include "../TestModels/SinCosModel.hpp"
24 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
25 
26 #include <fstream>
27 #include <vector>
28 
29 namespace Tempus_Test {
30 
31 using Teuchos::RCP;
32 using Teuchos::ParameterList;
33 using Teuchos::sublist;
34 using Teuchos::getParametersFromXmlFile;
35 
36 using Tempus::IntegratorBasic;
37 using Tempus::SolutionHistory;
38 using Tempus::SolutionState;
39 
40 // ************************************************************
41 // ************************************************************
42 void test_sincos_fsa(const bool use_combined_method,
43  const bool use_dfdp_as_tangent,
44  Teuchos::FancyOStream &out, bool &success)
45 {
46  std::vector<std::string> RKMethods;
47  RKMethods.push_back("RK Forward Euler");
48  RKMethods.push_back("RK Explicit 4 Stage");
49  RKMethods.push_back("RK Explicit 3/8 Rule");
50  RKMethods.push_back("RK Explicit 4 Stage 3rd order by Runge");
51  RKMethods.push_back("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
52  RKMethods.push_back("RK Explicit 3 Stage 3rd order");
53  RKMethods.push_back("RK Explicit 3 Stage 3rd order TVD");
54  RKMethods.push_back("RK Explicit 3 Stage 3rd order by Heun");
55  RKMethods.push_back("RK Explicit 2 Stage 2nd order by Runge");
56  RKMethods.push_back("RK Explicit Trapezoidal");
57  RKMethods.push_back("General ERK");
58  std::vector<double> RKMethodErrors;
59  if (use_combined_method) {
60  RKMethodErrors.push_back(0.183799);
61  RKMethodErrors.push_back(6.88637e-06);
62  RKMethodErrors.push_back(6.88637e-06);
63  RKMethodErrors.push_back(0.000264154);
64  RKMethodErrors.push_back(5.22798e-05);
65  RKMethodErrors.push_back(0.000261896);
66  RKMethodErrors.push_back(0.000261896);
67  RKMethodErrors.push_back(0.000261896);
68  RKMethodErrors.push_back(0.00934377);
69  RKMethodErrors.push_back(0.00934377);
70  RKMethodErrors.push_back(6.88637e-06);
71  }
72  else {
73  RKMethodErrors.push_back(0.183799);
74  RKMethodErrors.push_back(2.1915e-05);
75  RKMethodErrors.push_back(2.23367e-05);
76  RKMethodErrors.push_back(0.000205051);
77  RKMethodErrors.push_back(2.85141e-05);
78  RKMethodErrors.push_back(0.000126478);
79  RKMethodErrors.push_back(9.64964e-05);
80  RKMethodErrors.push_back(0.000144616);
81  RKMethodErrors.push_back(0.00826159);
82  RKMethodErrors.push_back(0.00710492);
83  RKMethodErrors.push_back(2.1915e-05);
84  }
85  Teuchos::RCP<const Teuchos::Comm<int> > comm =
86  Teuchos::DefaultComm<int>::getComm();
87  Teuchos::RCP<Teuchos::FancyOStream> my_out =
88  Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
89  my_out->setProcRankAndSize(comm->getRank(), comm->getSize());
90  my_out->setOutputToRootOnly(0);
91 
92  for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
93 
94  std::string RKMethod_ = RKMethods[m];
95  std::replace(RKMethod_.begin(), RKMethod_.end(), ' ', '_');
96  std::replace(RKMethod_.begin(), RKMethod_.end(), '/', '.');
97  std::vector<double> StepSize;
98  std::vector<double> ErrorNorm;
99  const int nTimeStepSizes = 7;
100  double dt = 0.2;
101  double order = 0.0;
102  for (int n=0; n<nTimeStepSizes; n++) {
103 
104  // Read params from .xml file
105  RCP<ParameterList> pList =
106  getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
107 
108  // Setup the SinCosModel
109  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
110  scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent);
111  RCP<SinCosModel<double> > model =
112  Teuchos::rcp(new SinCosModel<double>(scm_pl));
113 
114  // Set the Stepper
115  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
116  if (RKMethods[m] == "General ERK") {
117  pl->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper 2");
118  } else {
119  pl->sublist("Demo Stepper").set("Stepper Type", RKMethods[m]);
120  }
121 
122 
123  dt /= 2;
124 
125  // Setup sensitivities
126  ParameterList& sens_pl = pl->sublist("Sensitivities");
127  if (use_combined_method)
128  sens_pl.set("Sensitivity Method", "Combined");
129  else
130  sens_pl.set("Sensitivity Method", "Staggered");
131  sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
132  ParameterList& interp_pl =
133  pl->sublist("Demo Integrator").sublist("Solution History").sublist("Interpolator");
134  interp_pl.set("Interpolator Type", "Lagrange");
135  interp_pl.set("Order", 3); // All RK methods here are at most 4th order
136 
137  // Setup the Integrator and reset initial time step
138  pl->sublist("Demo Integrator")
139  .sublist("Time Step Control").set("Initial Time Step", dt);
140  RCP<Tempus::IntegratorForwardSensitivity<double> > integrator =
141  Tempus::integratorForwardSensitivity<double>(pl, model);
142  order = integrator->getStepper()->getOrder();
143 
144  // Initial Conditions
145  double t0 = pl->sublist("Demo Integrator")
146  .sublist("Time Step Control").get<double>("Initial Time");
147  // RCP<const Thyra::VectorBase<double> > x0 =
148  // model->getExactSolution(t0).get_x()->clone_v();
149  RCP<Thyra::VectorBase<double> > x0 =
150  model->getNominalValues().get_x()->clone_v();
151  const int num_param = model->get_p_space(0)->dim();
152  RCP<Thyra::MultiVectorBase<double> > DxDp0 =
153  Thyra::createMembers(model->get_x_space(), num_param);
154  for (int i=0; i<num_param; ++i)
155  Thyra::assign(DxDp0->col(i).ptr(),
156  *(model->getExactSensSolution(i, t0).get_x()));
157  integrator->setInitialState(t0, x0, Teuchos::null, Teuchos::null,
158  DxDp0, Teuchos::null, Teuchos::null);
159 
160  // Integrate to timeMax
161  bool integratorStatus = integrator->advanceTime();
162  TEST_ASSERT(integratorStatus)
163 
164  // Test if at 'Final Time'
165  double time = integrator->getTime();
166  double timeFinal = pl->sublist("Demo Integrator")
167  .sublist("Time Step Control").get<double>("Final Time");
168  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
169 
170  // Time-integrated solution and the exact solution
171  RCP<const Thyra::VectorBase<double> > x = integrator->getX();
172  RCP<const Thyra::MultiVectorBase<double> > DxDp = integrator->getDxDp();
173  RCP<const Thyra::VectorBase<double> > x_exact =
174  model->getExactSolution(time).get_x();
175  RCP<Thyra::MultiVectorBase<double> > DxDp_exact =
176  Thyra::createMembers(model->get_x_space(), num_param);
177  for (int i=0; i<num_param; ++i)
178  Thyra::assign(DxDp_exact->col(i).ptr(),
179  *(model->getExactSensSolution(i, time).get_x()));
180 
181  // Plot sample solution and exact solution
182  if (comm->getRank() == 0 && n == nTimeStepSizes-1) {
183  typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
184 
185  std::ofstream ftmp("Tempus_"+RKMethod_+"_SinCos_Sens.dat");
186  RCP<const SolutionHistory<double> > solutionHistory =
187  integrator->getSolutionHistory();
188  RCP< Thyra::MultiVectorBase<double> > DxDp_exact_plot =
189  Thyra::createMembers(model->get_x_space(), num_param);
190  for (int i=0; i<solutionHistory->getNumStates(); i++) {
191  RCP<const SolutionState<double> > solutionState =
192  (*solutionHistory)[i];
193  double time = solutionState->getTime();
194  RCP<const DMVPV> x_prod_plot =
195  Teuchos::rcp_dynamic_cast<const DMVPV>(solutionState->getX());
196  RCP<const Thyra::VectorBase<double> > x_plot =
197  x_prod_plot->getMultiVector()->col(0);
198  RCP<const Thyra::MultiVectorBase<double> > DxDp_plot =
199  x_prod_plot->getMultiVector()->subView(Teuchos::Range1D(1,num_param));
200  RCP<const Thyra::VectorBase<double> > x_exact_plot =
201  model->getExactSolution(time).get_x();
202  for (int j=0; j<num_param; ++j)
203  Thyra::assign(DxDp_exact_plot->col(j).ptr(),
204  *(model->getExactSensSolution(j, time).get_x()));
205  ftmp << std::fixed << std::setprecision(7)
206  << time
207  << std::setw(11) << get_ele(*(x_plot), 0)
208  << std::setw(11) << get_ele(*(x_plot), 1);
209  for (int j=0; j<num_param; ++j)
210  ftmp << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 0)
211  << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 1);
212  ftmp << std::setw(11) << get_ele(*(x_exact_plot), 0)
213  << std::setw(11) << get_ele(*(x_exact_plot), 1);
214  for (int j=0; j<num_param; ++j)
215  ftmp << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 0)
216  << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 1);
217  ftmp << std::endl;
218  }
219  ftmp.close();
220  }
221 
222  // Calculate the error
223  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
224  RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
225  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
226  Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
227  StepSize.push_back(dt);
228  double L2norm = Thyra::norm_2(*xdiff);
229  L2norm *= L2norm;
230  Teuchos::Array<double> L2norm_DxDp(num_param);
231  Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
232  for (int i=0; i<num_param; ++i)
233  L2norm += L2norm_DxDp[i]*L2norm_DxDp[i];
234  L2norm = std::sqrt(L2norm);
235  ErrorNorm.push_back(L2norm);
236 
237  *my_out << " n = " << n << " dt = " << dt << " error = " << L2norm
238  << std::endl;
239  }
240 
241  // Check the order and intercept
242  double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
243  *my_out << " Stepper = " << RKMethods[m] << std::endl;
244  *my_out << " =========================" << std::endl;
245  *my_out << " Expected order: " << order << std::endl;
246  *my_out << " Observed order: " << slope << std::endl;
247  *my_out << " =========================" << std::endl;
248  TEST_FLOATING_EQUALITY( slope, order, 0.04 );
249  TEST_FLOATING_EQUALITY( ErrorNorm[0], RKMethodErrors[m], 1.0e-4 );
250 
251  if (comm->getRank() == 0) {
252  std::ofstream ftmp("Tempus_"+RKMethod_+"_SinCos_Sens-Error.dat");
253  double error0 = 0.8*ErrorNorm[0];
254  for (int n=0; n<nTimeStepSizes; n++) {
255  ftmp << StepSize[n] << " " << ErrorNorm[n] << " "
256  << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
257  }
258  ftmp.close();
259  }
260  }
261 
262  Teuchos::TimeMonitor::summarize();
263 }
264 
265 } // namespace Tempus_Test
Tempus_Test
Definition: Tempus_BackwardEuler_ASA.cpp:34
Tempus::solutionHistory
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
Definition: Tempus_SolutionHistory_impl.hpp:504
Tempus::IntegratorBasic
Basic time integrator.
Definition: Tempus_IntegratorBasic_decl.hpp:35
Tempus::SolutionState
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Definition: Tempus_SolutionState_decl.hpp:56
Tempus_Test::test_sincos_fsa
void test_sincos_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)
Definition: Tempus_BackwardEuler_FSA.hpp:48
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:25