12 #include "Teuchos_DefaultComm.hpp"
14 #include "Tempus_config.hpp"
15 #include "Tempus_IntegratorBasic.hpp"
16 #include "Tempus_StepperBDF2.hpp"
18 #ifdef TEMPUS_ENABLE_EPETRA_STACK
19 #include "../TestModels/CDR_Model.hpp"
20 #ifdef Tempus_ENABLE_MPI
21 #include "Epetra_MpiComm.h"
23 #include "Epetra_SerialComm.h"
26 #ifdef TEMPUS_ENABLE_TPETRA_STACK
27 #include "../TestModels/CDR_Model_Tpetra.hpp"
28 #include "Tpetra_Core.hpp"
30 #include "../TestModels/VanDerPolModel.hpp"
31 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
33 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
40 namespace Tempus_Test {
42 using Teuchos::getParametersFromXmlFile;
46 using Teuchos::rcp_const_cast;
47 using Teuchos::sublist;
55 template <
typename SC,
typename Model,
typename Comm>
59 RCP<Tempus::IntegratorBasic<double>> integrator;
60 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
61 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
62 std::vector<double> StepSize;
65 auto pList = getParametersFromXmlFile(
"Tempus_BDF2_CDR.xml");
69 auto pl = sublist(pList,
"Tempus",
true);
70 auto dt = pl->sublist(
"Demo Integrator")
71 .sublist(
"Time Step Control")
72 .get<
double>(
"Initial Time Step");
74 auto model_pl = sublist(pList,
"CDR Model",
true);
76 const auto nTimeStepSizes =
77 model_pl->get<
int>(
"Number of Time Step Sizes", 5);
79 for (
int n = 0; n < nTimeStepSizes; n++) {
81 const int num_elements = model_pl->get<
int>(
"num elements");
82 const auto left_end = model_pl->get<SC>(
"left end");
83 const auto right_end = model_pl->get<SC>(
"right end");
84 const auto a_convection = model_pl->get<SC>(
"a (convection)");
85 const auto k_source = model_pl->get<SC>(
"k (source)");
87 auto model =
rcp(
new Model(comm, num_elements, left_end, right_end,
88 a_convection, k_source));
91 ::Stratimikos::DefaultLinearSolverBuilder builder;
93 auto p =
rcp(
new ParameterList);
94 p->set(
"Linear Solver Type",
"Belos");
95 p->set(
"Preconditioner Type",
"None");
96 builder.setParameterList(p);
98 auto lowsFactory = builder.createLinearSolveStrategy(
"");
100 model->set_W_factory(lowsFactory);
106 pl->sublist(
"Demo Integrator")
107 .sublist(
"Time Step Control")
108 .set(
"Initial Time Step", dt);
109 integrator = Tempus::createIntegratorBasic<double>(pl, model);
112 bool integratorStatus = integrator->advanceTime();
116 double time = integrator->getTime();
117 double timeFinal = pl->sublist("Demo Integrator")
118 .sublist("Time Step Control")
119 .get<
double>("Final Time");
120 double tol = 100.0 * std::numeric_limits<
double>::epsilon();
121 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
124 StepSize.push_back(dt);
125 auto solution = Thyra::createMember(model->get_x_space());
126 Thyra::copy(*(integrator->getX()), solution.ptr());
127 solutions.push_back(solution);
128 auto solutionDot = Thyra::createMember(model->get_x_space());
129 Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
130 solutionsDot.push_back(solutionDot);
134 if ((n == nTimeStepSizes - 1) && (commSize == 1)) {
135 std::ofstream ftmp(
"Tempus_BDF2_CDR.dat");
136 ftmp <<
"TITLE=\"BDF2 Solution to CDR\"\n"
137 <<
"VARIABLES=\"z\",\"T\"\n";
139 std::fabs(left_end - right_end) /
static_cast<double>(num_elements);
140 auto solutionHistory = integrator->getSolutionHistory();
141 int nStates = solutionHistory->getNumStates();
142 for (
int i = 0; i < nStates; i++) {
143 auto solutionState = (*solutionHistory)[i];
144 auto x = solutionState->getX();
145 auto ttime = solutionState->getTime();
146 ftmp <<
"ZONE T=\"Time=" << ttime <<
"\", I=" << num_elements + 1
148 for (
int j = 0; j < num_elements + 1; j++) {
149 const auto x_coord = left_end +
static_cast<double>(j) * dx;
150 ftmp << x_coord <<
" ";
153 for (
int j = 0; j < num_elements + 1; j++)
154 ftmp << get_ele(*x, j) <<
" ";
162 if (nTimeStepSizes > 2) {
164 double xDotSlope = 0.0;
165 std::vector<double> xErrorNorm;
166 std::vector<double> xDotErrorNorm;
167 auto stepper = integrator->getStepper();
168 auto order = stepper->getOrder();
169 writeOrderError(
"Tempus_BDF2_CDR-Error.dat", stepper, StepSize, solutions,
170 xErrorNorm, xSlope, solutionsDot, xDotErrorNorm, xDotSlope,
184 auto pListCDR = getParametersFromXmlFile(
"Tempus_BDF2_CDR.xml");
185 auto model_pl_CDR = sublist(pListCDR,
"CDR Model",
true);
186 const auto num_elements = model_pl_CDR->get<
int>(
"num elements");
187 const auto left_end = model_pl_CDR->get<
double>(
"left end");
188 const auto right_end = model_pl_CDR->get<
double>(
"right end");
190 const auto& x = *(solutions[solutions.size() - 1]);
192 std::ofstream ftmp(
"Tempus_BDF2_CDR-Solution.dat");
193 for (
int n = 0; n < num_elements + 1; n++) {
195 std::fabs(left_end - right_end) /
static_cast<double>(num_elements);
196 const auto x_coord = left_end +
static_cast<double>(n) * dx;
197 ftmp << x_coord <<
" " << Thyra::get_ele(x, n) << std::endl;
205 #ifdef TEMPUS_ENABLE_EPETRA_STACK
211 RCP<Epetra_Comm> comm;
212 #ifdef Tempus_ENABLE_MPI
213 comm =
rcp(
new Epetra_MpiComm(MPI_COMM_WORLD));
215 comm =
rcp(
new Epetra_SerialComm);
218 CDR_Test<double, Tempus_Test::CDR_Model<double>>(comm, comm->NumProc(), out,
223 #ifdef TEMPUS_ENABLE_TPETRA_STACK
229 using SC = Tpetra::Vector<>::scalar_type;
230 using LO = Tpetra::Vector<>::local_ordinal_type;
231 using GO = Tpetra::Vector<>::global_ordinal_type;
232 using Node = Tpetra::Vector<>::node_type;
234 auto comm = Tpetra::getDefaultComm();
236 CDR_Test<SC, Tempus_Test::CDR_Model_Tpetra<SC, LO, GO, Node>>(
237 comm, comm->getSize(), out, success);
#define TEST_COMPARE(v1, comp, v2)
void writeOrderError(const std::string filename, Teuchos::RCP< Tempus::Stepper< Scalar >> stepper, std::vector< Scalar > &StepSize, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutions, std::vector< Scalar > &xErrorNorm, Scalar &xSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutionsDot, std::vector< Scalar > &xDotErrorNorm, Scalar &xDotSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutionsDotDot, std::vector< Scalar > &xDotDotErrorNorm, Scalar &xDotDotSlope, Teuchos::FancyOStream &out)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
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)
void CDR_Test(const Comm &comm, const int commSize, Teuchos::FancyOStream &out, bool &success)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers.