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 {
31 using Teuchos::getParametersFromXmlFile;
34 using Teuchos::sublist;
44 std::vector<std::string> RKMethods;
45 RKMethods.push_back(
"RK Forward Euler");
46 RKMethods.push_back(
"RK Explicit 4 Stage");
47 RKMethods.push_back(
"RK Explicit 3/8 Rule");
48 RKMethods.push_back(
"RK Explicit 4 Stage 3rd order by Runge");
49 RKMethods.push_back(
"RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
50 RKMethods.push_back(
"RK Explicit 3 Stage 3rd order");
51 RKMethods.push_back(
"RK Explicit 3 Stage 3rd order TVD");
52 RKMethods.push_back(
"RK Explicit 3 Stage 3rd order by Heun");
53 RKMethods.push_back(
"RK Explicit Midpoint");
54 RKMethods.push_back(
"RK Explicit Trapezoidal");
55 RKMethods.push_back(
"Heuns Method");
56 RKMethods.push_back(
"General ERK");
58 std::vector<double> RKMethodErrors;
59 RKMethodErrors.push_back(0.154904);
60 RKMethodErrors.push_back(4.55982e-06);
61 RKMethodErrors.push_back(4.79132e-06);
62 RKMethodErrors.push_back(0.000113603);
63 RKMethodErrors.push_back(4.98796e-05);
64 RKMethodErrors.push_back(0.00014564);
65 RKMethodErrors.push_back(0.000121968);
66 RKMethodErrors.push_back(0.000109495);
67 RKMethodErrors.push_back(0.00559871);
68 RKMethodErrors.push_back(0.00710492);
69 RKMethodErrors.push_back(0.00710492);
70 RKMethodErrors.push_back(4.55982e-06);
75 for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
76 std::string RKMethod_ = RKMethods[m];
77 std::replace(RKMethod_.begin(), RKMethod_.end(),
' ',
'_');
78 std::replace(RKMethod_.begin(), RKMethod_.end(),
'/',
'.');
79 std::vector<double> StepSize;
80 std::vector<double> ErrorNorm;
81 const int nTimeStepSizes = 6;
84 for (
int n = 0; n < nTimeStepSizes; n++) {
87 getParametersFromXmlFile(
"Tempus_ExplicitRK_SinCos.xml");
96 if (RKMethods[m] ==
"General ERK") {
97 pl->
sublist(
"Demo Integrator").
set(
"Stepper Name",
"Demo Stepper 2");
99 .
set(
"Initial Condition Consistency",
"None");
101 .
set(
"Initial Condition Consistency Check",
false);
104 pl->
sublist(
"Demo Stepper").
set(
"Stepper Type", RKMethods[m]);
106 .
set(
"Initial Condition Consistency",
"None");
108 .
set(
"Initial Condition Consistency Check",
false);
115 sens_pl.
set(
"Mass Matrix Is Identity",
true);
119 interp_pl.
set(
"Interpolator Type",
"Lagrange");
120 interp_pl.
set(
"Order", 3);
125 .
set(
"Initial Time Step", dt);
127 Tempus::createIntegratorAdjointSensitivity<double>(pl, model);
128 order = integrator->getStepper()->getOrder();
131 double t0 = pl->sublist(
"Demo Integrator")
132 .sublist(
"Time Step Control")
133 .
get<
double>(
"Initial Time");
137 model->getNominalValues().get_x()->clone_v();
138 const int num_param = model->get_p_space(0)->dim();
140 Thyra::createMembers(model->get_x_space(), num_param);
141 for (
int i = 0; i < num_param; ++i)
142 Thyra::assign(DxDp0->col(i).
ptr(),
143 *(model->getExactSensSolution(i, t0).get_x()));
144 integrator->initializeSolutionHistory(t0, x0, Teuchos::null,
145 Teuchos::null, DxDp0, Teuchos::null,
149 bool integratorStatus = integrator->advanceTime();
153 double time = integrator->getTime();
154 double timeFinal = pl->sublist(
"Demo Integrator")
155 .sublist(
"Time Step Control")
156 .
get<
double>(
"Final Time");
165 Thyra::createMembers(model->get_x_space(), num_param);
169 const int num_g = DgDp->domain()->dim();
170 for (
int i = 0; i < num_g; ++i)
171 for (
int j = 0; j < num_param; ++j) dxdp_view(i, j) = dgdp_view(j, i);
174 model->getExactSolution(time).get_x();
176 Thyra::createMembers(model->get_x_space(), num_param);
177 for (
int i = 0; i < num_param; ++i)
178 Thyra::assign(DxDp_exact->col(i).
ptr(),
179 *(model->getExactSensSolution(i, time).get_x()));
182 if (comm->getRank() == 0 && n == nTimeStepSizes - 1) {
186 std::ofstream ftmp(
"Tempus_" + RKMethod_ +
"_SinCos_AdjSens.dat");
188 integrator->getSolutionHistory();
189 for (
int i = 0; i < solutionHistory->getNumStates(); i++) {
191 (*solutionHistory)[i];
192 const double time_i = solutionState->getTime();
194 Teuchos::rcp_dynamic_cast<
const DPV>(solutionState->getX());
196 x_prod_plot->getVectorBlock(0);
198 Teuchos::rcp_dynamic_cast<
const DMVPV>(
199 x_prod_plot->getVectorBlock(1));
201 adjoint_prod_plot->getMultiVector();
203 model->getExactSolution(time_i).get_x();
204 ftmp << std::fixed << std::setprecision(7) << time_i << std::setw(11)
205 << get_ele(*(x_plot), 0) << std::setw(11)
206 << get_ele(*(x_plot), 1) << std::setw(11)
207 << get_ele(*(adjoint_plot->col(0)), 0) << std::setw(11)
208 << get_ele(*(adjoint_plot->col(0)), 1) << std::setw(11)
209 << get_ele(*(adjoint_plot->col(1)), 0) << std::setw(11)
210 << get_ele(*(adjoint_plot->col(1)), 1) << std::setw(11)
211 << get_ele(*(x_exact_plot), 0) << std::setw(11)
212 << get_ele(*(x_exact_plot), 1) << std::endl;
220 Thyra::V_StVpStV(xdiff.
ptr(), 1.0, *x_exact, -1.0, *(x));
221 Thyra::V_VmV(DxDpdiff.
ptr(), *DxDp_exact, *DxDp);
222 StepSize.push_back(dt);
223 double L2norm = Thyra::norm_2(*xdiff);
226 Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
227 for (
int i = 0; i < num_param; ++i)
228 L2norm += L2norm_DxDp[i] * L2norm_DxDp[i];
229 L2norm = std::sqrt(L2norm);
230 ErrorNorm.push_back(L2norm);
237 double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
238 out <<
" Stepper = " << RKMethods[m] << std::endl;
239 out <<
" =========================" << std::endl;
240 out <<
" Expected order: " << order << std::endl;
241 out <<
" Observed order: " << slope << std::endl;
242 out <<
" =========================" << std::endl;
246 if (comm->getRank() == 0) {
247 std::ofstream ftmp(
"Tempus_" + RKMethod_ +
"_SinCos_AdjSens-Error.dat");
248 double error0 = 0.8 * ErrorNorm[0];
249 for (
int n = 0; n < nTimeStepSizes; n++) {
250 ftmp << StepSize[n] <<
" " << ErrorNorm[n] <<
" "
251 << error0 * (pow(StepSize[n] / StepSize[0], order)) << 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_UNIT_TEST(BackwardEuler, SinCos_ASA)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
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.