Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_BDF2_ASA.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_IntegratorAdjointSensitivity.hpp"
17
18#include "Thyra_VectorStdOps.hpp"
19#include "Thyra_MultiVectorStdOps.hpp"
20
21#include "../TestModels/SinCosModel.hpp"
22#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
23
24#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
25#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
26#include "Thyra_DefaultMultiVectorProductVector.hpp"
27
28#include <fstream>
29#include <limits>
30#include <sstream>
31#include <vector>
32
33namespace Tempus_Test {
34
35using Teuchos::RCP;
36using Teuchos::ParameterList;
37using Teuchos::sublist;
38using Teuchos::getParametersFromXmlFile;
39
43
44// ************************************************************
45// ************************************************************
46TEUCHOS_UNIT_TEST(BDF2, SinCos_ASA)
47{
48 std::vector<double> StepSize;
49 std::vector<double> ErrorNorm;
50 const int nTimeStepSizes = 7;
51 double dt = 0.2;
52 double order = 0.0;
53 Teuchos::RCP<const Teuchos::Comm<int> > comm =
54 Teuchos::DefaultComm<int>::getComm();
55 Teuchos::RCP<Teuchos::FancyOStream> my_out =
56 Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
57 my_out->setProcRankAndSize(comm->getRank(), comm->getSize());
58 my_out->setOutputToRootOnly(0);
59 for (int n=0; n<nTimeStepSizes; n++) {
60
61 // Read params from .xml file
62 RCP<ParameterList> pList =
63 getParametersFromXmlFile("Tempus_BDF2_SinCos_SA.xml");
64
65 // Setup the SinCosModel
66 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
67 RCP<SinCosModel<double> > model =
68 Teuchos::rcp(new SinCosModel<double>(scm_pl));
69
70 dt /= 2;
71
72 // Setup sensitivities
73 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
74 //ParameterList& sens_pl = pl->sublist("Sensitivities");
75 ParameterList& interp_pl =
76 pl->sublist("Default Integrator").sublist("Solution History")
77 .sublist("Interpolator");
78 interp_pl.set("Interpolator Type", "Lagrange");
79 interp_pl.set("Order", 1);
80
81 // Setup the Integrator and reset initial time step
82 pl->sublist("Default Integrator")
83 .sublist("Time Step Control").set("Initial Time Step", dt);
84 RCP<Tempus::IntegratorAdjointSensitivity<double> > integrator =
85 Tempus::createIntegratorAdjointSensitivity<double>(pl, model);
86 order = integrator->getStepper()->getOrder();
87
88 // Initial Conditions
89 double t0 = pl->sublist("Default Integrator")
90 .sublist("Time Step Control").get<double>("Initial Time");
91 RCP<const Thyra::VectorBase<double> > x0 =
92 model->getExactSolution(t0).get_x();
93 const int num_param = model->get_p_space(0)->dim();
94 RCP<Thyra::MultiVectorBase<double> > DxDp0 =
95 Thyra::createMembers(model->get_x_space(), num_param);
96 for (int i=0; i<num_param; ++i)
97 Thyra::assign(DxDp0->col(i).ptr(),
98 *(model->getExactSensSolution(i, t0).get_x()));
99 integrator->initializeSolutionHistory(t0, x0, Teuchos::null, Teuchos::null,
100 DxDp0, Teuchos::null, Teuchos::null);
101
102 // Integrate to timeMax
103 bool integratorStatus = integrator->advanceTime();
104 TEST_ASSERT(integratorStatus)
105
106 // Test if at 'Final Time'
107 double time = integrator->getTime();
108 double timeFinal =pl->sublist("Default Integrator")
109 .sublist("Time Step Control").get<double>("Final Time");
110 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
111
112 // Time-integrated solution and the exact solution along with
113 // sensitivities (relying on response g(x) = x). Note we must transpose
114 // dg/dp since the integrator returns it in gradient form.
115 RCP<const Thyra::VectorBase<double> > x = integrator->getX();
116 RCP<const Thyra::MultiVectorBase<double> > DgDp = integrator->getDgDp();
117 RCP<Thyra::MultiVectorBase<double> > DxDp =
118 Thyra::createMembers(model->get_x_space(), num_param);
119 {
120 Thyra::ConstDetachedMultiVectorView<double> dgdp_view(*DgDp);
121 Thyra::DetachedMultiVectorView<double> dxdp_view(*DxDp);
122 const int num_g = DgDp->domain()->dim();
123 for (int i=0; i<num_g; ++i)
124 for (int j=0; j<num_param; ++j)
125 dxdp_view(i,j) = dgdp_view(j,i);
126 }
127 RCP<const Thyra::VectorBase<double> > x_exact =
128 model->getExactSolution(time).get_x();
129 RCP<Thyra::MultiVectorBase<double> > DxDp_exact =
130 Thyra::createMembers(model->get_x_space(), num_param);
131 for (int i=0; i<num_param; ++i)
132 Thyra::assign(DxDp_exact->col(i).ptr(),
133 *(model->getExactSensSolution(i, time).get_x()));
134
135 // Plot sample solution and exact solution
136 if (comm->getRank() == 0 && n == nTimeStepSizes-1) {
137 typedef Thyra::DefaultProductVector<double> DPV;
138 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
139
140 std::ofstream ftmp("Tempus_BDF2_SinCos_AdjSens.dat");
141 RCP<const SolutionHistory<double> > solutionHistory =
142 integrator->getSolutionHistory();
143 for (int i=0; i<solutionHistory->getNumStates(); i++) {
144 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
145 const double time_i = solutionState->getTime();
146 RCP<const DPV> x_prod_plot =
147 Teuchos::rcp_dynamic_cast<const DPV>(solutionState->getX());
148 RCP<const Thyra::VectorBase<double> > x_plot =
149 x_prod_plot->getVectorBlock(0);
150 RCP<const DMVPV > adjoint_prod_plot =
151 Teuchos::rcp_dynamic_cast<const DMVPV>(x_prod_plot->getVectorBlock(1));
152 RCP<const Thyra::MultiVectorBase<double> > adjoint_plot =
153 adjoint_prod_plot->getMultiVector();
154 RCP<const Thyra::VectorBase<double> > x_exact_plot =
155 model->getExactSolution(time_i).get_x();
156 ftmp << std::fixed << std::setprecision(7)
157 << time_i
158 << std::setw(11) << get_ele(*(x_plot), 0)
159 << std::setw(11) << get_ele(*(x_plot), 1)
160 << std::setw(11) << get_ele(*(adjoint_plot->col(0)), 0)
161 << std::setw(11) << get_ele(*(adjoint_plot->col(0)), 1)
162 << std::setw(11) << get_ele(*(adjoint_plot->col(1)), 0)
163 << std::setw(11) << get_ele(*(adjoint_plot->col(1)), 1)
164 << std::setw(11) << get_ele(*(x_exact_plot), 0)
165 << std::setw(11) << get_ele(*(x_exact_plot), 1)
166 << std::endl;
167 }
168 ftmp.close();
169 }
170
171 // Calculate the error
172 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
173 RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
174 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
175 Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
176 StepSize.push_back(dt);
177 double L2norm = Thyra::norm_2(*xdiff);
178 L2norm *= L2norm;
179 Teuchos::Array<double> L2norm_DxDp(num_param);
180 Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
181 for (int i=0; i<num_param; ++i)
182 L2norm += L2norm_DxDp[i]*L2norm_DxDp[i];
183 L2norm = std::sqrt(L2norm);
184 ErrorNorm.push_back(L2norm);
185
186 //*my_out << " n = " << n << " dt = " << dt << " error = " << L2norm
187 // << std::endl;
188 }
189
190 // Check the order and intercept
191 double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
192 *my_out << " Stepper = BDF2" << std::endl;
193 *my_out << " =========================" << std::endl;
194 *my_out << " Expected order: " << order << std::endl;
195 *my_out << " Observed order: " << slope << std::endl;
196 *my_out << " =========================" << std::endl;
197 TEST_FLOATING_EQUALITY( slope, 1.95006, 0.015 );
198 TEST_FLOATING_EQUALITY( ErrorNorm[0], 0.0137394, 1.0e-4 );
199
200 if (comm->getRank() == 0) {
201 std::ofstream ftmp("Tempus_BDF2_SinCos_AdjSens-Error.dat");
202 double error0 = 0.8*ErrorNorm[0];
203 for (int n=0; n<nTimeStepSizes; n++) {
204 ftmp << StepSize[n] << " " << ErrorNorm[n] << " "
205 << error0*(StepSize[n]/StepSize[0]) << std::endl;
206 }
207 ftmp.close();
208 }
209
210}
211
212} // namespace Tempus_Test
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)