9 #include "Teuchos_UnitTestHarness.hpp"
10 #include "Teuchos_XMLParameterListHelpers.hpp"
11 #include "Teuchos_TimeMonitor.hpp"
12 #include "Teuchos_DefaultComm.hpp"
14 #include "Thyra_VectorStdOps.hpp"
15 #include "Thyra_MultiVectorStdOps.hpp"
17 #include "Tempus_IntegratorBasic.hpp"
18 #include "Tempus_IntegratorAdjointSensitivity.hpp"
20 #include "Thyra_DefaultMultiVectorProductVector.hpp"
21 #include "Thyra_DefaultProductVector.hpp"
23 #include "../TestModels/SinCosModel.hpp"
24 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
29 namespace Tempus_Test {
32 using Teuchos::ParameterList;
33 using Teuchos::sublist;
34 using Teuchos::getParametersFromXmlFile;
44 std::vector<std::string> RKMethods;
45 RKMethods.push_back(
"RK Backward Euler");
46 RKMethods.push_back(
"IRK 1 Stage Theta Method");
47 RKMethods.push_back(
"SDIRK 1 Stage 1st order");
48 RKMethods.push_back(
"SDIRK 2 Stage 2nd order");
49 RKMethods.push_back(
"SDIRK 2 Stage 3rd order");
50 RKMethods.push_back(
"EDIRK 2 Stage 3rd order");
51 RKMethods.push_back(
"EDIRK 2 Stage Theta Method");
52 RKMethods.push_back(
"SDIRK 3 Stage 4th order");
53 RKMethods.push_back(
"SDIRK 5 Stage 4th order");
54 RKMethods.push_back(
"SDIRK 5 Stage 5th order");
56 std::vector<double> RKMethodErrors;
57 RKMethodErrors.push_back(0.0383339);
58 RKMethodErrors.push_back(0.000221028);
59 RKMethodErrors.push_back(0.0383339);
60 RKMethodErrors.push_back(8.48235e-05);
61 RKMethodErrors.push_back(4.87848e-06);
62 RKMethodErrors.push_back(7.30827e-07);
63 RKMethodErrors.push_back(0.0144662);
64 RKMethodErrors.push_back(3.10132e-07);
65 RKMethodErrors.push_back(7.56838e-10);
66 RKMethodErrors.push_back(1.32374e-10);
68 Teuchos::RCP<const Teuchos::Comm<int> > comm =
69 Teuchos::DefaultComm<int>::getComm();
70 Teuchos::RCP<Teuchos::FancyOStream> my_out =
71 Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
72 my_out->setProcRankAndSize(comm->getRank(), comm->getSize());
73 my_out->setOutputToRootOnly(0);
75 for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
77 std::string RKMethod_ = RKMethods[m];
78 std::replace(RKMethod_.begin(), RKMethod_.end(),
' ',
'_');
79 std::replace(RKMethod_.begin(), RKMethod_.end(),
'/',
'.');
80 std::vector<double> StepSize;
81 std::vector<double> ErrorNorm;
82 const int nTimeStepSizes = 3;
85 for (
int n=0; n<nTimeStepSizes; n++) {
88 RCP<ParameterList> pList =
89 getParametersFromXmlFile(
"Tempus_DIRK_SinCos.xml");
92 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
93 RCP<SinCosModel<double> > model =
97 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
98 pl->sublist(
"Default Stepper").set(
"Stepper Type", RKMethods[m]);
99 if (RKMethods[m] ==
"IRK 1 Stage Theta Method" ||
100 RKMethods[m] ==
"EDIRK 2 Stage Theta Method") {
101 pl->sublist(
"Default Stepper").set<
double>(
"theta", 0.5);
102 }
else if (RKMethods[m] ==
"SDIRK 2 Stage 2nd order") {
103 pl->sublist(
"Default Stepper").set(
"gamma", 0.2928932188134524);
104 }
else if (RKMethods[m] ==
"SDIRK 2 Stage 3rd order") {
105 pl->sublist(
"Default Stepper").set(
"3rd Order A-stable",
true);
106 pl->sublist(
"Default Stepper").set(
"2nd Order L-stable",
false);
107 pl->sublist(
"Default Stepper").set(
"gamma", 0.7886751345948128);
113 ParameterList& sens_pl = pl->sublist(
"Sensitivities");
114 sens_pl.set(
"Mass Matrix Is Identity",
true);
115 ParameterList& interp_pl =
116 pl->sublist(
"Default Integrator").sublist(
"Solution History").sublist(
"Interpolator");
117 interp_pl.set(
"Interpolator Type",
"Lagrange");
118 interp_pl.set(
"Order", 4);
121 pl->sublist(
"Default Integrator")
122 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
123 RCP<Tempus::IntegratorAdjointSensitivity<double> > integrator =
124 Tempus::integratorAdjointSensitivity<double>(pl, model);
125 order = integrator->getStepper()->getOrder();
128 if (RKMethods[m] ==
"EDIRK 2 Stage Theta Method") order = 1.0;
134 RCP<Thyra::VectorBase<double> > x0 =
135 model->getNominalValues().get_x()->clone_v();
136 const int num_param = model->get_p_space(0)->dim();
137 RCP<Thyra::MultiVectorBase<double> > DxDp0 =
138 Thyra::createMembers(model->get_x_space(), num_param);
139 for (
int i=0; i<num_param; ++i)
140 Thyra::assign(DxDp0->col(i).ptr(),
141 *(model->getExactSensSolution(i, 0.0).get_x()));
142 integrator->initializeSolutionHistory(0.0, x0, Teuchos::null, Teuchos::null,
143 DxDp0, Teuchos::null, Teuchos::null);
146 bool integratorStatus = integrator->advanceTime();
147 TEST_ASSERT(integratorStatus)
150 double time = integrator->getTime();
151 double timeFinal = pl->sublist(
"Default Integrator")
152 .sublist(
"Time Step Control").get<
double>(
"Final Time");
153 double tol = 100.0 * std::numeric_limits<double>::epsilon();
154 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
159 RCP<const Thyra::VectorBase<double> > x = integrator->getX();
160 RCP<const Thyra::MultiVectorBase<double> > DgDp = integrator->getDgDp();
161 RCP<Thyra::MultiVectorBase<double> > DxDp =
162 Thyra::createMembers(model->get_x_space(), num_param);
164 Thyra::ConstDetachedMultiVectorView<double> dgdp_view(*DgDp);
165 Thyra::DetachedMultiVectorView<double> dxdp_view(*DxDp);
166 const int num_g = DgDp->domain()->dim();
167 for (
int i=0; i<num_g; ++i)
168 for (
int j=0; j<num_param; ++j)
169 dxdp_view(i,j) = dgdp_view(j,i);
171 RCP<const Thyra::VectorBase<double> > x_exact =
172 model->getExactSolution(time).get_x();
173 RCP<Thyra::MultiVectorBase<double> > DxDp_exact =
174 Thyra::createMembers(model->get_x_space(), num_param);
175 for (
int i=0; i<num_param; ++i)
176 Thyra::assign(DxDp_exact->col(i).ptr(),
177 *(model->getExactSensSolution(i, time).get_x()));
180 if (comm->getRank() == 0 && n == nTimeStepSizes-1) {
181 typedef Thyra::DefaultProductVector<double> DPV;
182 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
184 std::ofstream ftmp(
"Tempus_"+RKMethod_+
"_SinCos_AdjSens.dat");
186 integrator->getSolutionHistory();
187 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
188 RCP<const SolutionState<double> > solutionState =
189 (*solutionHistory)[i];
190 const double time_i = solutionState->getTime();
191 RCP<const DPV> x_prod_plot =
192 Teuchos::rcp_dynamic_cast<
const DPV>(solutionState->getX());
193 RCP<const Thyra::VectorBase<double> > x_plot =
194 x_prod_plot->getVectorBlock(0);
195 RCP<const DMVPV > adjoint_prod_plot =
196 Teuchos::rcp_dynamic_cast<
const DMVPV>(x_prod_plot->getVectorBlock(1));
197 RCP<const Thyra::MultiVectorBase<double> > adjoint_plot =
198 adjoint_prod_plot->getMultiVector();
199 RCP<const Thyra::VectorBase<double> > x_exact_plot =
200 model->getExactSolution(time_i).get_x();
201 ftmp << std::fixed << std::setprecision(7)
203 << std::setw(11) << get_ele(*(x_plot), 0)
204 << std::setw(11) << get_ele(*(x_plot), 1)
205 << std::setw(11) << get_ele(*(adjoint_plot->col(0)), 0)
206 << std::setw(11) << get_ele(*(adjoint_plot->col(0)), 1)
207 << std::setw(11) << get_ele(*(adjoint_plot->col(1)), 0)
208 << std::setw(11) << get_ele(*(adjoint_plot->col(1)), 1)
209 << std::setw(11) << get_ele(*(x_exact_plot), 0)
210 << std::setw(11) << get_ele(*(x_exact_plot), 1)
217 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
218 RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
219 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
220 Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
221 StepSize.push_back(dt);
222 double L2norm = Thyra::norm_2(*xdiff);
224 Teuchos::Array<double> L2norm_DxDp(num_param);
225 Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
226 for (
int i=0; i<num_param; ++i)
227 L2norm += L2norm_DxDp[i]*L2norm_DxDp[i];
228 L2norm = std::sqrt(L2norm);
229 ErrorNorm.push_back(L2norm);
231 *my_out <<
" n = " << n <<
" dt = " << dt <<
" error = " << L2norm
235 if (comm->getRank() == 0) {
236 std::ofstream ftmp(
"Tempus_"+RKMethod_+
"_SinCos_AdjSens-Error.dat");
237 double error0 = 0.8*ErrorNorm[0];
238 for (
int n=0; n<(int)StepSize.size(); n++) {
239 ftmp << StepSize[n] <<
" " << ErrorNorm[n] <<
" "
240 << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
254 double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
255 *my_out <<
" Stepper = " << RKMethods[m] << std::endl;
256 *my_out <<
" =========================" << std::endl;
257 *my_out <<
" Expected order: " << order << std::endl;
258 *my_out <<
" Observed order: " << slope << std::endl;
259 *my_out <<
" =========================" << std::endl;
260 TEST_FLOATING_EQUALITY( slope, order, 0.03 );
261 TEST_FLOATING_EQUALITY( ErrorNorm[0], RKMethodErrors[m], 5.0e-4 );
264 Teuchos::TimeMonitor::summarize();
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
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...