12 #include "Teuchos_DefaultComm.hpp"
14 #include "Tempus_config.hpp"
15 #include "Tempus_IntegratorBasic.hpp"
16 #include "Tempus_IntegratorForwardSensitivity.hpp"
17 #include "Tempus_IntegratorPseudoTransientForwardSensitivity.hpp"
19 #include "Thyra_VectorStdOps.hpp"
20 #include "Thyra_MultiVectorStdOps.hpp"
22 #include "../TestModels/SinCosModel.hpp"
23 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
25 #include "Thyra_DefaultMultiVectorProductVector.hpp"
32 namespace Tempus_Test {
34 using Teuchos::getParametersFromXmlFile;
37 using Teuchos::sublist;
49 std::vector<double> StepSize;
50 std::vector<double> ErrorNorm;
51 const int nTimeStepSizes = 7;
56 for (
int n = 0; n < nTimeStepSizes; n++) {
59 getParametersFromXmlFile(
"Tempus_BackwardEuler_SinCos.xml");
63 scm_pl->
set(
"Use DfDp as Tangent", use_dfdp_as_tangent);
72 if (use_combined_method)
73 sens_pl.
set(
"Sensitivity Method",
"Combined");
75 sens_pl.
set(
"Sensitivity Method",
"Staggered");
76 sens_pl.
set(
"Reuse State Linear Solver",
true);
78 sens_pl.
set(
"Use DfDp as Tangent", use_dfdp_as_tangent);
81 pl->
sublist(
"Default Integrator")
83 .
set(
"Initial Time Step", dt);
85 Tempus::createIntegratorForwardSensitivity<double>(pl, model);
86 order = integrator->getStepper()->getOrder();
89 double t0 = pl->sublist(
"Default Integrator")
90 .sublist(
"Time Step Control")
91 .
get<
double>(
"Initial Time");
93 model->getExactSolution(t0).get_x();
94 const int num_param = model->get_p_space(0)->dim();
96 Thyra::createMembers(model->get_x_space(), num_param);
97 for (
int i = 0; i < num_param; ++i)
98 Thyra::assign(DxDp0->col(i).
ptr(),
99 *(model->getExactSensSolution(i, t0).get_x()));
100 integrator->initializeSolutionHistory(t0, x0, Teuchos::null, Teuchos::null,
101 DxDp0, Teuchos::null, Teuchos::null);
104 bool integratorStatus = integrator->advanceTime();
108 double time = integrator->getTime();
109 double timeFinal = pl->sublist(
"Default Integrator")
110 .sublist(
"Time Step Control")
111 .
get<
double>(
"Final Time");
118 model->getExactSolution(time).get_x();
120 Thyra::createMembers(model->get_x_space(), num_param);
121 for (
int i = 0; i < num_param; ++i)
122 Thyra::assign(DxDp_exact->col(i).
ptr(),
123 *(model->getExactSensSolution(i, time).get_x()));
126 if (comm->getRank() == 0 && n == nTimeStepSizes - 1) {
129 std::ofstream ftmp(
"Tempus_BackwardEuler_SinCos_Sens.dat");
131 integrator->getSolutionHistory();
133 Thyra::createMembers(model->get_x_space(), num_param);
134 for (
int i = 0; i < solutionHistory->getNumStates(); i++) {
136 double time_i = solutionState->getTime();
138 Teuchos::rcp_dynamic_cast<
const DMVPV>(solutionState->getX());
140 x_prod_plot->getMultiVector()->col(0);
142 x_prod_plot->getMultiVector()->subView(
145 model->getExactSolution(time_i).get_x();
146 for (
int j = 0; j < num_param; ++j)
147 Thyra::assign(DxDp_exact_plot->col(j).
ptr(),
148 *(model->getExactSensSolution(j, time_i).get_x()));
149 ftmp << std::fixed << std::setprecision(7) << time_i << std::setw(11)
150 << get_ele(*(x_plot), 0) << std::setw(11) << get_ele(*(x_plot), 1);
151 for (
int j = 0; j < num_param; ++j)
152 ftmp << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 0)
153 << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 1);
154 ftmp << std::setw(11) << get_ele(*(x_exact_plot), 0) << std::setw(11)
155 << get_ele(*(x_exact_plot), 1);
156 for (
int j = 0; j < num_param; ++j)
157 ftmp << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 0)
158 << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 1);
167 Thyra::V_StVpStV(xdiff.
ptr(), 1.0, *x_exact, -1.0, *(x));
168 Thyra::V_VmV(DxDpdiff.
ptr(), *DxDp_exact, *DxDp);
169 StepSize.push_back(dt);
170 double L2norm = Thyra::norm_2(*xdiff);
173 Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
174 for (
int i = 0; i < num_param; ++i)
175 L2norm += L2norm_DxDp[i] * L2norm_DxDp[i];
176 L2norm = std::sqrt(L2norm);
177 ErrorNorm.push_back(L2norm);
179 out <<
" n = " << n <<
" dt = " << dt <<
" error = " << L2norm << std::endl;
183 double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
184 out <<
" Stepper = BackwardEuler" << std::endl;
185 out <<
" =========================" << std::endl;
186 out <<
" Expected order: " << order << std::endl;
187 out <<
" Observed order: " << slope << std::endl;
188 out <<
" =========================" << std::endl;
192 if (comm->getRank() == 0) {
193 std::ofstream ftmp(
"Tempus_BackwardEuler_SinCos_Sens-Error.dat");
194 double error0 = 0.8 * ErrorNorm[0];
195 for (
int n = 0; n < nTimeStepSizes; n++) {
196 ftmp << StepSize[n] <<
" " << ErrorNorm[n] <<
" "
197 << error0 * (StepSize[n] / StepSize[0]) << std::endl;
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
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...
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Solution state for integrators and steppers.