12 #include "Teuchos_DefaultComm.hpp"
14 #include "Tempus_config.hpp"
15 #include "Tempus_IntegratorBasic.hpp"
16 #include "Tempus_IntegratorForwardSensitivity.hpp"
18 #include "Thyra_VectorStdOps.hpp"
19 #include "Thyra_MultiVectorStdOps.hpp"
21 #include "../TestModels/SinCosModel.hpp"
22 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
24 #include "Thyra_DefaultMultiVectorProductVector.hpp"
31 namespace Tempus_Test {
33 using Teuchos::getParametersFromXmlFile;
36 using Teuchos::sublist;
48 std::vector<double> StepSize;
49 std::vector<double> ErrorNorm;
50 const int nTimeStepSizes = 7;
55 for (
int n = 0; n < nTimeStepSizes; n++) {
57 RCP<ParameterList> pList =
58 getParametersFromXmlFile(
"Tempus_BDF2_SinCos_SA.xml");
61 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
62 scm_pl->set(
"Use DfDp as Tangent", use_dfdp_as_tangent);
63 RCP<SinCosModel<double> > model =
69 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
70 ParameterList& sens_pl = pl->sublist(
"Sensitivities");
71 if (use_combined_method)
72 sens_pl.set(
"Sensitivity Method",
"Combined");
74 sens_pl.set(
"Sensitivity Method",
"Staggered");
77 sens_pl.set(
"Use DfDp as Tangent", use_dfdp_as_tangent);
78 ParameterList& interp_pl = pl->sublist(
"Default Integrator")
79 .sublist(
"Solution History")
80 .sublist(
"Interpolator");
81 interp_pl.set(
"Interpolator Type",
"Lagrange");
82 interp_pl.set(
"Order", 1);
85 pl->sublist(
"Default Integrator")
86 .sublist(
"Time Step Control")
87 .set(
"Initial Time Step", dt);
88 RCP<Tempus::IntegratorForwardSensitivity<double> > integrator =
89 Tempus::createIntegratorForwardSensitivity<double>(pl, model);
90 order = integrator->getStepper()->getOrder();
93 double t0 = pl->sublist(
"Default Integrator")
94 .sublist(
"Time Step Control")
95 .get<
double>(
"Initial Time");
96 RCP<const Thyra::VectorBase<double> > x0 =
97 model->getExactSolution(t0).get_x();
98 const int num_param = model->get_p_space(0)->dim();
99 RCP<Thyra::MultiVectorBase<double> > DxDp0 =
100 Thyra::createMembers(model->get_x_space(), num_param);
101 for (
int i = 0; i < num_param; ++i)
102 Thyra::assign(DxDp0->col(i).ptr(),
103 *(model->getExactSensSolution(i, t0).get_x()));
104 integrator->initializeSolutionHistory(t0, x0, Teuchos::null, Teuchos::null,
105 DxDp0, Teuchos::null, Teuchos::null);
108 bool integratorStatus = integrator->advanceTime();
112 double time = integrator->getTime();
113 double timeFinal = pl->sublist("Default Integrator")
114 .sublist("Time Step Control")
115 .get<
double>("Final Time");
116 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
119 RCP<const Thyra::VectorBase<
double> > x = integrator->getX();
120 RCP<const Thyra::MultiVectorBase<
double> > DxDp = integrator->getDxDp();
121 RCP<const Thyra::VectorBase<
double> > x_exact =
122 model->getExactSolution(time).get_x();
123 RCP<Thyra::MultiVectorBase<
double> > DxDp_exact =
124 Thyra::createMembers(model->get_x_space(), num_param);
125 for (
int i = 0; i < num_param; ++i)
126 Thyra::assign(DxDp_exact->col(i).ptr(),
127 *(model->getExactSensSolution(i, time).get_x()));
130 if (comm->getRank() == 0 && n == nTimeStepSizes - 1) {
133 std::ofstream ftmp(
"Tempus_BDF2_SinCos_Sens.dat");
134 RCP<const SolutionHistory<double> > solutionHistory =
135 integrator->getSolutionHistory();
136 RCP<Thyra::MultiVectorBase<double> > DxDp_exact_plot =
137 Thyra::createMembers(model->get_x_space(), num_param);
138 for (
int i = 0; i < solutionHistory->getNumStates(); i++) {
139 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
140 double time_i = solutionState->getTime();
141 RCP<const DMVPV> x_prod_plot =
142 Teuchos::rcp_dynamic_cast<
const DMVPV>(solutionState->getX());
143 RCP<const Thyra::VectorBase<double> > x_plot =
144 x_prod_plot->getMultiVector()->col(0);
145 RCP<const Thyra::MultiVectorBase<double> > DxDp_plot =
146 x_prod_plot->getMultiVector()->subView(
148 RCP<const Thyra::VectorBase<double> > x_exact_plot =
149 model->getExactSolution(time_i).get_x();
150 for (
int j = 0; j < num_param; ++j)
151 Thyra::assign(DxDp_exact_plot->col(j).ptr(),
152 *(model->getExactSensSolution(j, time_i).get_x()));
153 ftmp << std::fixed << std::setprecision(7) << time_i << std::setw(11)
154 << get_ele(*(x_plot), 0) << std::setw(11) << get_ele(*(x_plot), 1);
155 for (
int j = 0; j < num_param; ++j)
156 ftmp << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 0)
157 << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 1);
158 ftmp << std::setw(11) << get_ele(*(x_exact_plot), 0) << std::setw(11)
159 << get_ele(*(x_exact_plot), 1);
160 for (
int j = 0; j < num_param; ++j)
161 ftmp << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 0)
162 << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 1);
169 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
170 RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
171 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
172 Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
173 StepSize.push_back(dt);
174 double L2norm = Thyra::norm_2(*xdiff);
177 Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
178 for (
int i = 0; i < num_param; ++i)
179 L2norm += L2norm_DxDp[i] * L2norm_DxDp[i];
180 L2norm = std::sqrt(L2norm);
181 ErrorNorm.push_back(L2norm);
183 out <<
" n = " << n <<
" dt = " << dt <<
" error = " << L2norm << std::endl;
187 double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
188 out <<
" Stepper = BDF2" << std::endl;
189 out <<
" =========================" << std::endl;
190 out <<
" Expected order: " << order << std::endl;
191 out <<
" Observed order: " << slope << std::endl;
192 out <<
" =========================" << std::endl;
196 if (comm->getRank() == 0) {
197 std::ofstream ftmp(
"Tempus_BDF2_SinCos_Sens-Error.dat");
198 double error0 = 0.8 * ErrorNorm[0];
199 for (
int n = 0; n < nTimeStepSizes; n++) {
200 ftmp << StepSize[n] <<
" " << ErrorNorm[n] <<
" "
201 << error0 * (StepSize[n] / StepSize[0]) << std::endl;
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void test_sincos_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers.