Skip to content
Snippets Groups Projects
Commit d67c78dd authored by Christoph Lehmann's avatar Christoph Lehmann
Browse files

[T] compare Newton to Picard

parent 2ddbf99c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Tests/NumLib/TestODEInt.cpp
+ 88
− 91
View file @ d67c78dd
...@@ -24,6 +24,12 @@ ...@@ -24,6 +24,12 @@
using GMatrix = GlobalMatrix; using GMatrix = GlobalMatrix;
using GVector = GlobalVector; using GVector = GlobalVector;
template <typename Vector>
struct Solution
{
std::vector<double> ts;
std::vector<Vector> solutions;
};
template<typename Matrix, typename Vector, NumLib::NonlinearSolverTag NLTag> template<typename Matrix, typename Vector, NumLib::NonlinearSolverTag NLTag>
class TestOutput class TestOutput
...@@ -36,18 +42,15 @@ public: ...@@ -36,18 +42,15 @@ public:
: _file_name_part(name) : _file_name_part(name)
{} {}
template<typename ODE>
void run_test(ODE& ode, TimeDisc& timeDisc)
{
run_test<ODE>(ode, timeDisc, 10); // by default make 10 timesteps
}
template <class ODE> template <class ODE>
void run_test(ODE& ode, TimeDisc& timeDisc, const unsigned num_timesteps) Solution<Vector> run_test(ODE& ode, TimeDisc& timeDisc,
const unsigned num_timesteps)
{ {
using ODE_ = ODE; using ODE_ = ODE;
using ODET = ODETraits<ODE>; using ODET = ODETraits<ODE>;
Solution<Vector> sol;
NumLib::TimeDiscretizedODESystem<ODE_::ODETag, NLTag> NumLib::TimeDiscretizedODESystem<ODE_::ODETag, NLTag>
ode_sys(ode, timeDisc); ode_sys(ode, timeDisc);
...@@ -70,23 +73,30 @@ public: ...@@ -70,23 +73,30 @@ public:
INFO("Running test %s with %u timesteps of size %g s.", INFO("Running test %s with %u timesteps of size %g s.",
_file_name_part.c_str(), num_timesteps, delta_t); _file_name_part.c_str(), num_timesteps, delta_t);
init_file(delta_t); // init_file(delta_t);
// initial condition // initial condition
Vector x0(ode.getMatrixSpecifications().nrows); Vector x0(ode.getMatrixSpecifications().nrows);
ODET::setIC(x0); ODET::setIC(x0);
write(t0, x0, x0); // write(t0, x0, x0);
sol.ts.push_back(t0);
sol.solutions.push_back(x0);
auto cb = [this](const double t, Vector const& x) { auto cb = [this, &sol](const double t, Vector const& x) {
loopCallback<ODE>(t, x); // loopCallback<ODE>(t, x);
sol.ts.push_back(t);
sol.solutions.push_back(x);
}; };
if (num_timesteps > 0) if (num_timesteps > 0)
EXPECT_TRUE(loop.loop(t0, x0, t_end, delta_t, cb)); EXPECT_TRUE(loop.loop(t0, x0, t_end, delta_t, cb));
return sol;
} }
private: private:
/*
void init_file(const double delta_t) void init_file(const double delta_t)
{ {
std::string path(BaseLib::BuildInfo::tests_tmp_path + "ODEInt_"); std::string path(BaseLib::BuildInfo::tests_tmp_path + "ODEInt_");
...@@ -113,6 +123,7 @@ private: ...@@ -113,6 +123,7 @@ private:
{ {
write(t, x, ODETraits<Ode>::solution(t)); write(t, x, ODETraits<Ode>::solution(t));
} }
*/
const std::string _file_name_part; const std::string _file_name_part;
std::unique_ptr<std::ofstream> _file; std::unique_ptr<std::ofstream> _file;
...@@ -121,72 +132,71 @@ private: ...@@ -121,72 +132,71 @@ private:
const unsigned _maxiter = 20; const unsigned _maxiter = 20;
}; };
template <typename Matrix, typename Vector, typename TimeDisc, typename ODE,
template<typename Matrix, typename Vector, typename TimeDisc, typename ODE, NumLib::NonlinearSolverTag NLTag>
NumLib::NonlinearSolverTag NLTag> typename std::enable_if<std::is_same<TimeDisc, NumLib::BackwardEuler>::value,
typename std::enable_if<std::is_same<TimeDisc, NumLib::BackwardEuler >::value>::type Solution<Vector>>::type
run_test_case(const unsigned num_timesteps, const char* name) run_test_case(const unsigned num_timesteps, const char* name)
{ {
ODE ode; ODE ode;
TimeDisc timeDisc; TimeDisc timeDisc;
TestOutput<Matrix, Vector, NLTag> test(name); TestOutput<Matrix, Vector, NLTag> test(name);
test.run_test(ode, timeDisc, num_timesteps); return test.run_test(ode, timeDisc, num_timesteps);
} }
template<typename Matrix, typename Vector, typename TimeDisc, typename ODE, template <typename Matrix, typename Vector, typename TimeDisc, typename ODE,
NumLib::NonlinearSolverTag NLTag> NumLib::NonlinearSolverTag NLTag>
typename std::enable_if<std::is_same<TimeDisc, NumLib::ForwardEuler >::value>::type typename std::enable_if<std::is_same<TimeDisc, NumLib::ForwardEuler>::value,
Solution<Vector>>::type
run_test_case(const unsigned num_timesteps, const char* name) run_test_case(const unsigned num_timesteps, const char* name)
{ {
ODE ode; ODE ode;
TimeDisc timeDisc; TimeDisc timeDisc;
TestOutput<Matrix, Vector, NLTag> test(name); TestOutput<Matrix, Vector, NLTag> test(name);
test.run_test(ode, timeDisc, num_timesteps); return test.run_test(ode, timeDisc, num_timesteps);
} }
template<typename Matrix, typename Vector, typename TimeDisc, typename ODE, template <typename Matrix, typename Vector, typename TimeDisc, typename ODE,
NumLib::NonlinearSolverTag NLTag> NumLib::NonlinearSolverTag NLTag>
typename std::enable_if<std::is_same<TimeDisc, NumLib::CrankNicolson >::value>::type typename std::enable_if<std::is_same<TimeDisc, NumLib::CrankNicolson>::value,
Solution<Vector>>::type
run_test_case(const unsigned num_timesteps, const char* name) run_test_case(const unsigned num_timesteps, const char* name)
{ {
ODE ode; ODE ode;
TimeDisc timeDisc(0.5); TimeDisc timeDisc(0.5);
TestOutput<Matrix, Vector, NLTag> test(name); TestOutput<Matrix, Vector, NLTag> test(name);
test.run_test(ode, timeDisc, num_timesteps); return test.run_test(ode, timeDisc, num_timesteps);
} }
template<typename Matrix, typename Vector, typename TimeDisc, typename ODE, template <typename Matrix, typename Vector, typename TimeDisc, typename ODE,
NumLib::NonlinearSolverTag NLTag> NumLib::NonlinearSolverTag NLTag>
typename std::enable_if< typename std::enable_if<
std::is_same<TimeDisc, NumLib::BackwardDifferentiationFormula >::value>::type std::is_same<TimeDisc, NumLib::BackwardDifferentiationFormula>::value,
Solution<Vector>>::type
run_test_case(const unsigned num_timesteps, const char* name) run_test_case(const unsigned num_timesteps, const char* name)
{ {
ODE ode; ODE ode;
TimeDisc timeDisc(3); TimeDisc timeDisc(3);
TestOutput<Matrix, Vector, NLTag> test(name); TestOutput<Matrix, Vector, NLTag> test(name);
test.run_test(ode, timeDisc, num_timesteps); return test.run_test(ode, timeDisc, num_timesteps);
} }
// This class is only here s.t. I don't have to put the members into // This class is only here s.t. I don't have to put the members into
// the definition of the macro TCLITEM below. // the definition of the macro TCLITEM below.
template <typename Matrix_, typename Vector_, class ODE_, class TimeDisc_, template <typename Matrix_, typename Vector_, class ODE_, class TimeDisc_>
NumLib::NonlinearSolverTag NLTag_>
struct TestCaseBase struct TestCaseBase
{ {
using Matrix = Matrix_; using Matrix = Matrix_;
using Vector = Vector_; using Vector = Vector_;
using ODE = ODE_; using ODE = ODE_;
using TimeDisc = TimeDisc_; using TimeDisc = TimeDisc_;
static const NumLib::NonlinearSolverTag NLTag = NLTag_;
}; };
template <typename Matrix_, typename Vector_, class ODE_, class TimeDisc_, template <typename Matrix_, typename Vector_, class ODE_, class TimeDisc_>
NumLib::NonlinearSolverTag NLTag_>
struct TestCase; struct TestCase;
...@@ -197,47 +207,29 @@ struct TestCase; ...@@ -197,47 +207,29 @@ struct TestCase;
// ///////////////////////////////////// // /////////////////////////////////////
#define TESTCASESLIST \ #define TESTCASESLIST \
/* Global sparse matrix */ \ /* Global sparse matrix */ \
TCLITEM(GMatrix, GVector, ODE1, BackwardEuler, Newton) TCLSEP \ TCLITEM(GMatrix, GVector, ODE1, BackwardEuler ) TCLSEP \
TCLITEM(GMatrix, GVector, ODE1, ForwardEuler, Newton) TCLSEP \ TCLITEM(GMatrix, GVector, ODE1, ForwardEuler ) TCLSEP \
TCLITEM(GMatrix, GVector, ODE1, CrankNicolson, Newton) TCLSEP \ TCLITEM(GMatrix, GVector, ODE1, CrankNicolson ) TCLSEP \
TCLITEM(GMatrix, GVector, ODE1, BackwardDifferentiationFormula, Newton) TCLSEP \ TCLITEM(GMatrix, GVector, ODE1, BackwardDifferentiationFormula) TCLSEP \
\
TCLITEM(GMatrix, GVector, ODE1, BackwardEuler, Picard) TCLSEP \
TCLITEM(GMatrix, GVector, ODE1, ForwardEuler, Picard) TCLSEP \
TCLITEM(GMatrix, GVector, ODE1, CrankNicolson, Picard) TCLSEP \
TCLITEM(GMatrix, GVector, ODE1, BackwardDifferentiationFormula, Picard) TCLSEP \
\ \
TCLITEM(GMatrix, GVector, ODE2, BackwardEuler, Newton) TCLSEP \ TCLITEM(GMatrix, GVector, ODE2, BackwardEuler ) TCLSEP \
TCLITEM(GMatrix, GVector, ODE2, ForwardEuler, Newton) TCLSEP \ TCLITEM(GMatrix, GVector, ODE2, ForwardEuler ) TCLSEP \
TCLITEM(GMatrix, GVector, ODE2, CrankNicolson, Newton) TCLSEP \ TCLITEM(GMatrix, GVector, ODE2, CrankNicolson ) TCLSEP \
TCLITEM(GMatrix, GVector, ODE2, BackwardDifferentiationFormula, Newton) TCLSEP \ TCLITEM(GMatrix, GVector, ODE2, BackwardDifferentiationFormula) TCLSEP \
\ \
TCLITEM(GMatrix, GVector, ODE2, BackwardEuler, Picard) TCLSEP \ TCLITEM(GMatrix, GVector, ODE3, BackwardEuler ) TCLSEP \
TCLITEM(GMatrix, GVector, ODE2, ForwardEuler, Picard) TCLSEP \ TCLITEM(GMatrix, GVector, ODE3, ForwardEuler ) TCLSEP \
TCLITEM(GMatrix, GVector, ODE2, CrankNicolson, Picard) TCLSEP \ TCLITEM(GMatrix, GVector, ODE3, CrankNicolson ) TCLSEP \
TCLITEM(GMatrix, GVector, ODE2, BackwardDifferentiationFormula, Picard) TCLSEP \ TCLITEM(GMatrix, GVector, ODE3, BackwardDifferentiationFormula)
\
TCLITEM(GMatrix, GVector, ODE3, BackwardEuler, Newton) TCLSEP \ #define TCLITEM(MAT, VEC, ODE, TIMEDISC) \
TCLITEM(GMatrix, GVector, ODE3, ForwardEuler, Newton) TCLSEP \ template <> \
TCLITEM(GMatrix, GVector, ODE3, CrankNicolson, Newton) TCLSEP \ struct TestCase<MAT, VEC, ODE, NumLib::TIMEDISC> \
TCLITEM(GMatrix, GVector, ODE3, BackwardDifferentiationFormula, Newton) TCLSEP \ : TestCaseBase<MAT, VEC, ODE, NumLib::TIMEDISC> { \
\ static const char name[]; \
TCLITEM(GMatrix, GVector, ODE3, BackwardEuler, Picard) TCLSEP \ }; \
TCLITEM(GMatrix, GVector, ODE3, ForwardEuler, Picard) TCLSEP \ const char TestCase<MAT, VEC, ODE, NumLib::TIMEDISC>::name[] = \
TCLITEM(GMatrix, GVector, ODE3, CrankNicolson, Picard) TCLSEP \ #MAT "_" #VEC "_" #ODE "_" #TIMEDISC;
TCLITEM(GMatrix, GVector, ODE3, BackwardDifferentiationFormula, Picard)
#define TCLITEM(MAT, VEC, ODE, TIMEDISC, NLTAG) \
template<> \
struct TestCase<MAT, VEC, ODE, NumLib::TIMEDISC, NumLib::NonlinearSolverTag::NLTAG> \
: TestCaseBase<MAT, VEC, ODE, NumLib::TIMEDISC, NumLib::NonlinearSolverTag::NLTAG> \
{ \
static const char name[]; \
}; \
const char TestCase<MAT, VEC, ODE, NumLib::TIMEDISC, \
NumLib::NonlinearSolverTag::NLTAG>::name[] \
= #MAT "_" #VEC "_" #ODE "_" #TIMEDISC "_" #NLTAG;
#define TCLSEP #define TCLSEP
TESTCASESLIST TESTCASESLIST
...@@ -245,8 +237,8 @@ TESTCASESLIST ...@@ -245,8 +237,8 @@ TESTCASESLIST
#undef TCLITEM #undef TCLITEM
#undef TCLSEP #undef TCLSEP
#define TCLITEM(MAT, VEC, ODE, TIMEDISC, NLTAG) \ #define TCLITEM(MAT, VEC, ODE, TIMEDISC) \
TestCase<MAT, VEC, ODE, NumLib::TIMEDISC, NumLib::NonlinearSolverTag::NLTAG> TestCase<MAT, VEC, ODE, NumLib::TIMEDISC>
#define TCLSEP , #define TCLSEP ,
typedef ::testing::Types<TESTCASESLIST> TestCases; typedef ::testing::Types<TESTCASESLIST> TestCases;
...@@ -269,15 +261,31 @@ public: ...@@ -269,15 +261,31 @@ public:
static void test() static void test()
{ {
for (auto num_timesteps : { 100u }) const unsigned num_timesteps = 100;
{
run_test_case<Matrix, Vector, TimeDisc, ODE, NLTag>( auto const sol_picard =
num_timesteps, TestParams::name); run_test_case<Matrix, Vector, TimeDisc, ODE,
NumLib::NonlinearSolverTag::Picard>(num_timesteps,
TestParams::name);
auto const sol_newton =
run_test_case<Matrix, Vector, TimeDisc, ODE,
NumLib::NonlinearSolverTag::Newton>(num_timesteps,
TestParams::name);
const double tol_picard_newton = 2e-9;
ASSERT_EQ(sol_picard.ts.size(), sol_newton.ts.size());
for (std::size_t i = 0; i < sol_picard.ts.size(); ++i) {
ASSERT_EQ(sol_picard.ts[i], sol_newton.ts[i]);
for (std::size_t comp = 0; comp < sol_picard.solutions[i].size(); ++comp) {
EXPECT_NEAR(sol_picard.solutions[i][comp],
sol_newton.solutions[i][comp], tol_picard_newton);
}
} }
} }
}; };
TYPED_TEST_CASE(NumLibODEIntTyped, TestCases); TYPED_TEST_CASE(NumLibODEIntTyped, TestCases);
TYPED_TEST(NumLibODEIntTyped, T1) TYPED_TEST(NumLibODEIntTyped, T1)
...@@ -286,17 +294,6 @@ TYPED_TEST(NumLibODEIntTyped, T1) ...@@ -286,17 +294,6 @@ TYPED_TEST(NumLibODEIntTyped, T1)
} }
TEST(NumLibODEInt, ODE3)
{
const char* name = "dummy";
{
run_test_case<GMatrix, GVector, NumLib::BackwardEuler, ODE3,
NumLib::NonlinearSolverTag::Newton>(0u, name);
}
}
/* TODO Other possible test cases: /* TODO Other possible test cases:
* *
* * check that results are within a specified tolerance * * check that results are within a specified tolerance
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment