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_IntegratorForwardSensitivity.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 const bool use_combined_method,
45 const bool use_dfdp_as_tangent,
46 Teuchos::FancyOStream &out,
bool &success)
48 std::vector<std::string> RKMethods;
49 RKMethods.push_back(
"RK Backward Euler");
50 RKMethods.push_back(
"IRK 1 Stage Theta Method");
51 RKMethods.push_back(
"SDIRK 1 Stage 1st order");
52 RKMethods.push_back(
"SDIRK 2 Stage 2nd order");
53 RKMethods.push_back(
"SDIRK 2 Stage 3rd order");
54 RKMethods.push_back(
"EDIRK 2 Stage 3rd order");
55 RKMethods.push_back(
"EDIRK 2 Stage Theta Method");
56 RKMethods.push_back(
"SDIRK 3 Stage 4th order");
57 RKMethods.push_back(
"SDIRK 5 Stage 4th order");
58 RKMethods.push_back(
"SDIRK 5 Stage 5th order");
61 if (method_name !=
"") {
62 auto it = std::find(RKMethods.begin(), RKMethods.end(),
method_name);
63 TEUCHOS_TEST_FOR_EXCEPTION(it == RKMethods.end(), std::logic_error,
67 std::vector<double> RKMethodErrors;
68 if (use_combined_method) {
69 RKMethodErrors.push_back(0.0428449);
70 RKMethodErrors.push_back(0.000297933);
71 RKMethodErrors.push_back(0.0428449);
72 RKMethodErrors.push_back(0.000144507);
73 RKMethodErrors.push_back(8.65434e-06);
74 RKMethodErrors.push_back(1.3468e-06);
75 RKMethodErrors.push_back(0.000297933);
76 RKMethodErrors.push_back(5.44037e-07);
77 RKMethodErrors.push_back(2.77342e-09);
78 RKMethodErrors.push_back(1.21689e-10);
81 RKMethodErrors.push_back(0.0428449);
82 RKMethodErrors.push_back(0.000221049);
83 RKMethodErrors.push_back(0.0428449);
84 RKMethodErrors.push_back(0.000125232);
85 RKMethodErrors.push_back(4.79475e-06);
86 RKMethodErrors.push_back(9.63899e-07);
87 RKMethodErrors.push_back(0.0141323);
88 RKMethodErrors.push_back(2.9362e-07);
89 RKMethodErrors.push_back(9.20081e-08);
90 RKMethodErrors.push_back(9.16252e-08);
93 Teuchos::RCP<const Teuchos::Comm<int> > comm =
94 Teuchos::DefaultComm<int>::getComm();
95 Teuchos::RCP<Teuchos::FancyOStream> my_out =
96 Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
97 my_out->setProcRankAndSize(comm->getRank(), comm->getSize());
98 my_out->setOutputToRootOnly(0);
100 for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
103 if (method_name !=
"" && RKMethods[m] != method_name)
106 std::string RKMethod_ = RKMethods[m];
107 std::replace(RKMethod_.begin(), RKMethod_.end(),
' ',
'_');
108 std::replace(RKMethod_.begin(), RKMethod_.end(),
'/',
'.');
109 std::vector<double> StepSize;
110 std::vector<double> ErrorNorm;
111 const int nTimeStepSizes = 3;
114 for (
int n=0; n<nTimeStepSizes; n++) {
117 RCP<ParameterList> pList =
118 getParametersFromXmlFile(
"Tempus_DIRK_SinCos.xml");
121 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
122 scm_pl->set(
"Use DfDp as Tangent", use_dfdp_as_tangent);
123 RCP<SinCosModel<double> > model =
127 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
128 pl->sublist(
"Default Stepper").set(
"Stepper Type", RKMethods[m]);
129 if (RKMethods[m] ==
"SDIRK 2 Stage 2nd order") {
130 pl->sublist(
"Default Stepper").set(
"gamma", 0.2928932190);
131 }
else if (RKMethods[m] ==
"SDIRK 2 Stage 3rd order") {
132 pl->sublist(
"Default Stepper").set(
"3rd Order A-stable",
true);
133 pl->sublist(
"Default Stepper").set(
"2nd Order L-stable",
false);
134 pl->sublist(
"Default Stepper").set(
"gamma", 0.7886751347);
140 ParameterList& sens_pl = pl->sublist(
"Sensitivities");
141 if (use_combined_method)
142 sens_pl.set(
"Sensitivity Method",
"Combined");
144 sens_pl.set(
"Sensitivity Method",
"Staggered");
145 sens_pl.set(
"Use DfDp as Tangent", use_dfdp_as_tangent);
146 ParameterList& interp_pl =
147 pl->sublist(
"Default Integrator").sublist(
"Solution History").sublist(
"Interpolator");
148 interp_pl.set(
"Interpolator Type",
"Lagrange");
149 interp_pl.set(
"Order", 4);
152 pl->sublist(
"Default Integrator")
153 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
154 RCP<Tempus::IntegratorForwardSensitivity<double> > integrator =
155 Tempus::integratorForwardSensitivity<double>(pl, model);
156 order = integrator->getStepper()->getOrder();
159 if (use_combined_method ==
false &&
160 RKMethods[m] ==
"EDIRK 2 Stage Theta Method") order = 1.0;
166 RCP<Thyra::VectorBase<double> > x0 =
167 model->getNominalValues().get_x()->clone_v();
168 const int num_param = model->get_p_space(0)->dim();
169 RCP<Thyra::MultiVectorBase<double> > DxDp0 =
170 Thyra::createMembers(model->get_x_space(), num_param);
171 for (
int i=0; i<num_param; ++i)
172 Thyra::assign(DxDp0->col(i).ptr(),
173 *(model->getExactSensSolution(i, 0.0).get_x()));
174 integrator->initializeSolutionHistory(0.0, x0, Teuchos::null, Teuchos::null,
175 DxDp0, Teuchos::null, Teuchos::null);
178 bool integratorStatus = integrator->advanceTime();
179 TEST_ASSERT(integratorStatus)
182 double time = integrator->getTime();
183 double timeFinal = pl->sublist(
"Default Integrator")
184 .sublist(
"Time Step Control").get<
double>(
"Final Time");
185 double tol = 100.0 * std::numeric_limits<double>::epsilon();
186 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
189 RCP<const Thyra::VectorBase<double> > x = integrator->getX();
190 RCP<const Thyra::MultiVectorBase<double> > DxDp = integrator->getDxDp();
191 RCP<const Thyra::VectorBase<double> > x_exact =
192 model->getExactSolution(time).get_x();
193 RCP<Thyra::MultiVectorBase<double> > DxDp_exact =
194 Thyra::createMembers(model->get_x_space(), num_param);
195 for (
int i=0; i<num_param; ++i)
196 Thyra::assign(DxDp_exact->col(i).ptr(),
197 *(model->getExactSensSolution(i, time).get_x()));
200 if (comm->getRank() == 0 && n == nTimeStepSizes-1) {
201 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
203 std::ofstream ftmp(
"Tempus_"+RKMethod_+
"_SinCos_Sens.dat");
205 integrator->getSolutionHistory();
206 RCP< Thyra::MultiVectorBase<double> > DxDp_exact_plot =
207 Thyra::createMembers(model->get_x_space(), num_param);
208 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
209 RCP<const SolutionState<double> > solutionState =
210 (*solutionHistory)[i];
211 double time_i = solutionState->getTime();
212 RCP<const DMVPV> x_prod_plot =
213 Teuchos::rcp_dynamic_cast<
const DMVPV>(solutionState->getX());
214 RCP<const Thyra::VectorBase<double> > x_plot =
215 x_prod_plot->getMultiVector()->col(0);
216 RCP<const Thyra::MultiVectorBase<double> > DxDp_plot =
217 x_prod_plot->getMultiVector()->subView(Teuchos::Range1D(1,num_param));
218 RCP<const Thyra::VectorBase<double> > x_exact_plot =
219 model->getExactSolution(time_i).get_x();
220 for (
int j=0; j<num_param; ++j)
221 Thyra::assign(DxDp_exact_plot->col(j).ptr(),
222 *(model->getExactSensSolution(j, time_i).get_x()));
223 ftmp << std::fixed << std::setprecision(7)
225 << std::setw(11) << get_ele(*(x_plot), 0)
226 << std::setw(11) << get_ele(*(x_plot), 1);
227 for (
int j=0; j<num_param; ++j)
228 ftmp << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 0)
229 << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 1);
230 ftmp << std::setw(11) << get_ele(*(x_exact_plot), 0)
231 << std::setw(11) << get_ele(*(x_exact_plot), 1);
232 for (
int j=0; j<num_param; ++j)
233 ftmp << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 0)
234 << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 1);
241 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
242 RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
243 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
244 Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
245 StepSize.push_back(dt);
246 double L2norm = Thyra::norm_2(*xdiff);
248 Teuchos::Array<double> L2norm_DxDp(num_param);
249 Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
250 for (
int i=0; i<num_param; ++i)
251 L2norm += L2norm_DxDp[i]*L2norm_DxDp[i];
252 L2norm = std::sqrt(L2norm);
253 ErrorNorm.push_back(L2norm);
255 *my_out <<
" n = " << n <<
" dt = " << dt <<
" error = " << L2norm
259 if (comm->getRank() == 0) {
260 std::ofstream ftmp(
"Tempus_"+RKMethod_+
"_SinCos_Sens-Error.dat");
261 double error0 = 0.8*ErrorNorm[0];
262 for (
int n=0; n<(int)StepSize.size(); n++) {
263 ftmp << StepSize[n] <<
" " << ErrorNorm[n] <<
" "
264 << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
278 double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
279 *my_out <<
" Stepper = " << RKMethods[m] << std::endl;
280 *my_out <<
" =========================" << std::endl;
281 *my_out <<
" Expected order: " << order << std::endl;
282 *my_out <<
" Observed order: " << slope << std::endl;
283 *my_out <<
" =========================" << std::endl;
286 double order_expected = use_combined_method ? order : std::min(order,4.0);
287 TEST_FLOATING_EQUALITY( slope, order_expected, 0.03 );
288 TEST_FLOATING_EQUALITY( ErrorNorm[0], RKMethodErrors[m], 5.0e-4 );
291 Teuchos::TimeMonitor::summarize();
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
void test_sincos_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)
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...