Page MenuHomec4science
Files F90977472

AbstractSolver.cpp
No OneTemporary
Actions

File Metadata

Created
Wed, Nov 6, 14:22

AbstractSolver.cpp
View Options

//
// Created by lionel on 25.11.19.
//
#include <cassert>
#include <sstream>
#include <iostream>
#include "AbstractSolver.h"
AbstractSolver::AbstractSolver(){
stepSize=0;
initialTime=0;
finalTime=0;
dimension=0;
nSteps=0;
};
AbstractSolver::AbstractSolver(double t0,double tf,int n,Eigen::VectorXd& y0):AbstractSolver(){
setInitialTime(t0);
setFinalTime(tf);
setNumberOfSteps(n);
setInitialValue(y0);
}
AbstractSolver::~AbstractSolver(){};
void AbstractSolver::setInitialTime(const double t0) {
initialTime = t0;
}
void AbstractSolver::setFinalTime(const double tf) {
finalTime = tf;
}
void AbstractSolver::setStepSize(const double h) {
checkStepSize(h);
if (stepSize!=0){ // need to warn user that the stepsize was changed to avoid confusion
std::cerr << "[WARNING] : Step size " << stepSize << " was already assigned, new step size is now " << h << std::endl;
}
stepSize = h;
}
void AbstractSolver::checkStepSize(double h) {
if (h <= 0) { // negative step size will make our solver fail
throw std::domain_error("Step size h must be strictly positive");
}
}
void AbstractSolver::checkTimes() {
if (finalTime <= initialTime){ // final time must be larger for the solve method to work
throw std::domain_error("Initial time is not smaller than the final time.");
}
}
void AbstractSolver::setNumberOfSteps(int n){
if (n < 1){
throw std::domain_error("Number of step must be strictly positive ");
} else {
double h = (finalTime - initialTime)/ (double) n; // linearly spaced steps
setStepSize(h);
}
nSteps=n;
}
void AbstractSolver::checkSpaceDimension(int dim) {
if (dim < 1) {
throw std::domain_error("Dimension must be greater than 0");
} else{
if (dimension == 0) { // the dimension was never specified if it is 0, so we set it to dim
dimension = dim;
} else if (dimension != dim){ // return an error with the conflicting values of dimension
std::stringstream errorMessage;
errorMessage << "Invalid dimension : " << dim << ", it was already assigned " << dimension << ".";
throw std::length_error(errorMessage.str());
} else { // do nothing if the values of dimension match
;
}
}
}
void AbstractSolver::setInitialValue(Eigen::VectorXd & y0) {
// check it is a valid input
checkSpaceDimension(y0.size());
initialValue = y0;
}
void AbstractSolver::setRightHandSide(Eigen::MatrixXd & rhsMatrix, Eigen::VectorXd (*g_)(double t)){
// check the right hand side matrix is valid before storing it
try {
checkSpaceDimension(rhsMatrix.rows());
checkSpaceDimension(rhsMatrix.cols());
} catch (const std::length_error& e){
std::stringstream errorMessage; // creates a more complex message to help the user correcting his error
errorMessage << "The linear Matrix has dimension ("<< rhsMatrix.rows() <<","<< rhsMatrix.cols() <<
") which is invalid : "<< e.what();
throw std::length_error(errorMessage.str());
}
A = rhsMatrix;
// check if the function is callable and that it returns object of desired dimension
try {
// checks if no exception happens for final and initial time, hope no exception will happen for in between times
auto testVector = g_(initialTime);
auto testVectorCopy = g_(finalTime);
} catch(const std::exception& e) {
std::stringstream errorMessage;
errorMessage << "Function g(t) could not be evaluated, an exception was caught with message : " << e.what();
throw std::runtime_error(errorMessage.str());
} catch (...) { // any other type of error will make our solver fail
throw std::runtime_error("Function g(t) could not be evaluated. No exception was caught.");
}
// checks the output of g is in the required format
try {
checkSpaceDimension(g_(initialTime).size());
} catch (const std::length_error& e){
std::stringstream errorMessage;
errorMessage << "Function g(t) has invalid dimension in output : " << e.what();
throw std::length_error(errorMessage.str());
}
g = g_;
}
double AbstractSolver::getInitialTime() const{
return initialTime;
}
double AbstractSolver::getFinalTime() const{
return finalTime;
}
double AbstractSolver::getStepSize() const{
return stepSize;
}
int AbstractSolver::getSpaceDimension() const {
return dimension;
}
void AbstractSolver::solve(std::ostream &stream) {
auto y = initialValue;
double t = initialTime;
assert(initialTime < finalTime && stepSize > 0); // ensures next while loop will not be infinite
stream << t << " " << y.transpose() << std::endl; //write the initial time
while(t< finalTime){
step(y,t); // step updates the y values
t += stepSize;
stream << t << " " << y.transpose() << std::endl;
}
}
void AbstractSolver::setRightHandSide(Eigen::MatrixXd &rhsMatrix, std::function<Eigen::VectorXd(double)> g_) {
// check the right hand side matrix is valid before storing it
try {
checkSpaceDimension(rhsMatrix.rows());
checkSpaceDimension(rhsMatrix.cols());
} catch (const std::length_error& e){
std::stringstream errorMessage; // creates a more complex message to help the user correcting his error
errorMessage << "The linear Matrix has dimension ("<< rhsMatrix.rows() <<","<< rhsMatrix.cols() <<
") which is invalid : "<< e.what();
throw std::length_error(errorMessage.str());
}
A = rhsMatrix;
// check if the function is callable and that it returns object of desired dimension
try {
// checks if no exception happens for final and initial time, hope no exception will happen for in between times
auto testVector = g_(initialTime);
auto testVectorCopy = g_(finalTime);
} catch(const std::exception& e) {
std::stringstream errorMessage;
errorMessage << "Function g(t) could not be evaluated, an exception was caught with message : " << e.what();
throw std::runtime_error(errorMessage.str());
} catch (...) { // any other type of error will make our solver fail
throw std::runtime_error("Function g(t) could not be evaluated. No exception was caught.");
}
// checks the output of g is in the required format
try {
checkSpaceDimension(g_(initialTime).size());
} catch (const std::length_error& e){
std::stringstream errorMessage;
errorMessage << "Function g(t) has invalid dimension in output : " << e.what();
throw std::length_error(errorMessage.str());
}
g = g_;
}

Event Timeline

c4science · Help