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 {
33 using Teuchos::sublist;
34 using Teuchos::getParametersFromXmlFile;
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 Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
79 for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
81 std::string RKMethod_ = RKMethods[m];
82 std::replace(RKMethod_.begin(), RKMethod_.end(),
' ',
'_');
83 std::replace(RKMethod_.begin(), RKMethod_.end(),
'/',
'.');
84 std::vector<double> StepSize;
85 std::vector<double> ErrorNorm;
86 const int nTimeStepSizes = 6;
89 for (
int n=0; n<nTimeStepSizes; n++) {
93 getParametersFromXmlFile(
"Tempus_ExplicitRK_SinCos.xml");
102 if (RKMethods[m] ==
"General ERK") {
103 pl->
sublist(
"Demo Integrator").
set(
"Stepper Name",
"Demo Stepper 2");
104 pl->
sublist(
"Demo Stepper 2").
set(
"Initial Condition Consistency",
106 pl->
sublist(
"Demo Stepper 2").
set(
"Initial Condition Consistency Check",
109 pl->
sublist(
"Demo Stepper").
set(
"Stepper Type", RKMethods[m]);
110 pl->
sublist(
"Demo Stepper").
set(
"Initial Condition Consistency",
112 pl->
sublist(
"Demo Stepper").
set(
"Initial Condition Consistency Check",
121 sens_pl.
set(
"Mass Matrix Is Identity",
true);
124 interp_pl.
set(
"Interpolator Type",
"Lagrange");
125 interp_pl.
set(
"Order", 3);
129 .
sublist(
"Time Step Control").
set(
"Initial Time Step", dt);
131 Tempus::integratorAdjointSensitivity<double>(pl, model);
132 order = integrator->getStepper()->getOrder();
135 double t0 = pl->sublist(
"Demo Integrator")
136 .sublist(
"Time Step Control").
get<
double>(
"Initial Time");
140 model->getNominalValues().get_x()->clone_v();
141 const int num_param = model->get_p_space(0)->dim();
143 Thyra::createMembers(model->get_x_space(), num_param);
144 for (
int i=0; i<num_param; ++i)
145 Thyra::assign(DxDp0->col(i).
ptr(),
146 *(model->getExactSensSolution(i, t0).get_x()));
147 integrator->initializeSolutionHistory(t0, x0, Teuchos::null, Teuchos::null,
148 DxDp0, Teuchos::null, Teuchos::null);
151 bool integratorStatus = integrator->advanceTime();
155 double time = integrator->getTime();
156 double timeFinal = pl->sublist(
"Demo Integrator")
157 .sublist(
"Time Step Control").
get<
double>(
"Final Time");
166 Thyra::createMembers(model->get_x_space(), num_param);
170 const int num_g = DgDp->domain()->dim();
171 for (
int i=0; i<num_g; ++i)
172 for (
int j=0; j<num_param; ++j)
173 dxdp_view(i,j) = dgdp_view(j,i);
176 model->getExactSolution(time).get_x();
178 Thyra::createMembers(model->get_x_space(), num_param);
179 for (
int i=0; i<num_param; ++i)
180 Thyra::assign(DxDp_exact->col(i).
ptr(),
181 *(model->getExactSensSolution(i, time).get_x()));
184 if (comm->getRank() == 0 && n == nTimeStepSizes-1) {
188 std::ofstream ftmp(
"Tempus_"+RKMethod_+
"_SinCos_AdjSens.dat");
190 integrator->getSolutionHistory();
191 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
193 (*solutionHistory)[i];
194 const double time_i = solutionState->getTime();
196 Teuchos::rcp_dynamic_cast<
const DPV>(solutionState->getX());
198 x_prod_plot->getVectorBlock(0);
200 Teuchos::rcp_dynamic_cast<
const DMVPV>(x_prod_plot->getVectorBlock(1));
202 adjoint_prod_plot->getMultiVector();
204 model->getExactSolution(time_i).get_x();
205 ftmp << std::fixed << std::setprecision(7)
207 << std::setw(11) << get_ele(*(x_plot), 0)
208 << std::setw(11) << get_ele(*(x_plot), 1)
209 << std::setw(11) << get_ele(*(adjoint_plot->col(0)), 0)
210 << std::setw(11) << get_ele(*(adjoint_plot->col(0)), 1)
211 << std::setw(11) << get_ele(*(adjoint_plot->col(1)), 0)
212 << std::setw(11) << get_ele(*(adjoint_plot->col(1)), 1)
213 << std::setw(11) << get_ele(*(x_exact_plot), 0)
214 << std::setw(11) << get_ele(*(x_exact_plot), 1)
223 Thyra::V_StVpStV(xdiff.
ptr(), 1.0, *x_exact, -1.0, *(x));
224 Thyra::V_VmV(DxDpdiff.
ptr(), *DxDp_exact, *DxDp);
225 StepSize.push_back(dt);
226 double L2norm = Thyra::norm_2(*xdiff);
229 Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
230 for (
int i=0; i<num_param; ++i)
231 L2norm += L2norm_DxDp[i]*L2norm_DxDp[i];
232 L2norm = std::sqrt(L2norm);
233 ErrorNorm.push_back(L2norm);
240 double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
241 *my_out <<
" Stepper = " << RKMethods[m] << std::endl;
242 *my_out <<
" =========================" << std::endl;
243 *my_out <<
" Expected order: " << order << std::endl;
244 *my_out <<
" Observed order: " << slope << std::endl;
245 *my_out <<
" =========================" << std::endl;
249 if (comm->getRank() == 0) {
250 std::ofstream ftmp(
"Tempus_"+RKMethod_+
"_SinCos_AdjSens-Error.dat");
251 double error0 = 0.8*ErrorNorm[0];
252 for (
int n=0; n<nTimeStepSizes; n++) {
253 ftmp << StepSize[n] <<
" " << ErrorNorm[n] <<
" "
254 << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
basic_FancyOStream & setProcRankAndSize(const int procRank, const int numProcs)
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...
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...