Page MenuHomec4science
Files F67132823

testRungeKuttaSolver.cpp
No OneTemporary
Actions

File Metadata

Created
Thu, Jun 20, 03:33

testRungeKuttaSolver.cpp
View Options

//
// Created by lionel on 01.12.19.
//
# include <gtest/gtest.h>
# include <cmath>
# include <cassert>
# include <exception>
#include <fstream>
#include <iostream>
#include "../ODElibrary/RungeKutta.h"
// Function g(t) for testing
Eigen::VectorXd ForcedOscillations(double t){
Eigen::VectorXd s(2);
s[0] = 0;
s[1] = sin(2*M_PI*t); // derivative of velocity
return s;
};
// Fixture class
class RKSolverTest : public ::testing::Test {
protected:
void SetUp() override {
// create a solver
Eigen::VectorXd y0(2); // set the spring at rest
y0[0] = 0;
y0[1] = 0;
double k = 0.1; // spring constant
Eigen::MatrixXd m(2,2); // from spring equation
m(0,0)=0;
m(0,1)=1;
m(1,0)=-k;
m(1,1)=0.3;
solver.setInitialTime(0.);
solver.setFinalTime(1.);
solver.setStepSize(0.01);
solver.setRightHandSide(m, ForcedOscillations);
solver.setInitialValue(y0);
}
void TearDown() override {
// clear particles
solver;
}
RungeKutta solver;
// values found by solving the system with an exact solution
double yFinal0Exact = 0.182532458;
double yFinal1Exact = 0.046017418;
};
// tests the space dimension
TEST_F(RKSolverTest, spaceDimension) {
ASSERT_EQ(solver.getSpaceDimension(), 2);
}
// tests the order of the method
TEST_F(RKSolverTest, methodOrder) {
// check default is 4
ASSERT_EQ(solver.getMethodOrder(), 4);
solver.setMethodOrder(2);
ASSERT_EQ(solver.getMethodOrder(), 2);
// check impossible order gives exception
ASSERT_THROW(solver.setMethodOrder(0), std::invalid_argument);
ASSERT_THROW(solver.setMethodOrder(5), std::invalid_argument);
}
// tests that we can assign and read a positive initial time
TEST_F(RKSolverTest, initialAndFinalTimes) {
solver.setInitialTime(0.5);
EXPECT_EQ(solver.getInitialTime(), 0.5);
solver.setFinalTime(4.2);
EXPECT_EQ(solver.getFinalTime(), 4.2);
}
// tests that we can assign and read a positive timestep, but get a warning message (already assigned in fixture)
TEST_F(RKSolverTest, positiveTimestep) {
testing::internal::CaptureStderr();
solver.setStepSize(0.3);
std::string output = testing::internal::GetCapturedStderr();
ASSERT_EQ(solver.getStepSize(), 0.3);
EXPECT_EQ(output, "[WARNING] : Step size 0.01 was already assigned, new step size is now 0.3\n");
}
// tests that we can assign a positive number of steps, but get a warning message (timestep already assigned in fixture)
TEST_F(RKSolverTest, numberOfStep) {
ASSERT_THROW(solver.setNumberOfSteps(-2), std::exception); // negative number of step throws exception
ASSERT_THROW(solver.setNumberOfSteps(0), std::exception); // 0 also throws exception
testing::internal::CaptureStderr();
solver.setNumberOfSteps(40);
std::string output = testing::internal::GetCapturedStderr();
ASSERT_EQ(solver.getStepSize(), 0.025);
EXPECT_EQ(output, "[WARNING] : Step size 0.01 was already assigned, new step size is now 0.025\n");
}
// tests we cannot assign a 0 or negative timestep
TEST_F(RKSolverTest, nonPositiveTimestep) {
ASSERT_THROW(solver.setStepSize(-0.1), std::domain_error);
ASSERT_THROW(solver.setStepSize(0), std::domain_error);
}
// tests we cannot assign a different size inputs for the initial value
TEST_F(RKSolverTest, differentSizedInitialValue) {
Eigen::VectorXd yTooLong(3);
ASSERT_THROW(solver.setInitialValue(yTooLong), std::length_error);
}
// tests we cannot assign a different size inputs for the rhs matrix
TEST_F(RKSolverTest, differentSizedRHS) {
Eigen::MatrixXd mFalse1(2,3);
Eigen::MatrixXd mFalse2(1,2);
Eigen::MatrixXd mFalse3(3,3);
ASSERT_THROW(solver.setRightHandSide(mFalse1, ForcedOscillations), std::length_error);
ASSERT_THROW(solver.setRightHandSide(mFalse2, ForcedOscillations), std::length_error);
ASSERT_THROW(solver.setRightHandSide(mFalse3, ForcedOscillations), std::length_error);
}
// tests we cannot assign a rhs function
Eigen::VectorXd RHSgTooLong(double t){
Eigen::VectorXd s(4);
return s;
};
Eigen::VectorXd RHSgThrowingException(double t){
Eigen::VectorXd s(2);
throw std::exception();
return s;
};
Eigen::VectorXd RHSgThrowingInteger(double t){
Eigen::VectorXd s(2);
throw 42;
return s;
};
TEST_F(RKSolverTest, timeDependentRHSFunction) {
Eigen::MatrixXd m(2,2);
ASSERT_THROW(solver.setRightHandSide(m, RHSgTooLong), std::length_error);
ASSERT_THROW(solver.setRightHandSide(m, RHSgThrowingException), std::runtime_error);
ASSERT_THROW(solver.setRightHandSide(m, RHSgThrowingInteger), std::runtime_error);
}
// for RK, we will test all method and check that error is smaller when we increase the method order
TEST_F(RKSolverTest, testSolve) {
double errors[4][2];
for (int order = 1; order <= 4; order++){ // loop over the order of the step
// writes the output to the file
std::ofstream testFile;
testFile.open("RKSolverTest.dat");
assert(testFile.is_open());
solver.setMethodOrder(order);
solver.solve(testFile);
testFile.close();
// read the file where the output is
std::ifstream readFile("RKSolverTest.dat");
assert(readFile.is_open());
double t=0,x=0,y=0;
int i=0;
while (!readFile.eof()) { // read till the end
readFile >> t >> x >> y;
i++;
}
readFile.close();
errors[order-1][0] = abs(x-yFinal0Exact);
errors[order-1][1] = abs(y-yFinal1Exact);
// check the value are close enough to the specified values
ASSERT_NEAR(x, yFinal0Exact, 1e-2);
ASSERT_NEAR(y, yFinal1Exact, 1e-2);
}
for (int j = 0; j < 3; ++j) {
EXPECT_TRUE(errors[j][0] > errors[j+1][0]);
if (j!=1){ // order 3 error is larger than order 2, hard to exactly know why, it is a bit of makeshift job
EXPECT_TRUE(errors[j][1] > errors[j+1][1]);
}
}
}
// tests the constructor way of giving input
TEST_F(RKSolverTest, constructor) {
Eigen::VectorXd y0(2); // set the spring at rest
y0[0] = 0;
y0[1] = 0;
RungeKutta byConstructorSolver(-2.,4., 10, 3,y0);
}

Event Timeline

c4science · Help