ExternalVendorCode/trick
1
0
Fork 0
mirror of https://github.com/nasa/trick.git synced 2024-12-18 20:57:55 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Improve unittest code coverage for SAIntegrator library. #1081

Browse Source
This commit is contained in:
Penn, John M 047828115 2020-11-17 13:37:06 -06:00
parent f04dcd7567
commit 7b34af2e54
8 changed files with 668 additions and 349 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
trick_source/trick_utils/SAIntegrator/examples/CannonBall/CannonBall.cpp
View File

@ -27,15 +27,9 @@ void calc_derivs( double t, double state[], double derivs[], void* udata) {
double impact( double t, double state[], RootFinder* root_finder, void* udata) {
double root_error = root_finder->find_roots(t, state[1]);
if (root_error == 0.0) {
printf("---------------------------------------------------------------\n");
printf("Impact at t = %5.10f x = %5.10f y = %5.10f.\n", t, state[0], state[1]);
printf("---------------------------------------------------------------\n");
root_finder->init();
state[1] = 0.0000000001; // pos_y
state[2] = 0.0; // vel_x
state[3] = 0.0; // vel_y
state[4] = 0.0; // acc_x
state[5] = 0.0; // acc_y
state[2] = 0.9 * state[2];
state[3] = -0.9 * state[3];
}
return (root_error);
}
@ -55,17 +49,14 @@ int main ( int argc, char* argv[]) {
double t = 0.0;
SA::RK2Integrator integ(dt, 6, state_var_p, state_var_p, calc_derivs, NULL);
// Uncomment the following two lines to have the RootFinder detect impact.
// RootFinder root_finder(0.00000000001, Negative);
// integ.add_Rootfinder(&root_finder, &impact);
RootFinder root_finder(0.00000000001, Negative);
integ.add_Rootfinder(&root_finder, &impact);
init_state(cannon);
print_header();
print_state( t, cannon);
while (t < 5.3) {
integ.load();
integ.step();
integ.unload();
while (t < 20.0) {
integ.integrate();
t = integ.getIndyVar();
print_state( t, cannon);
}

8
trick_source/trick_utils/SAIntegrator/examples/DoubleIntegral/DoubleIntegral.cpp
View File

@ -23,9 +23,7 @@ void deriv_x( double x, double state_x[], double derivs_x[], void* udata) {
integ_y->set_out_vars(state_y);
integ_y->setIndyVar(context->start);
while (integ_y->getIndyVar() <= context->end) {
integ_y->load();
integ_y->step();
integ_y->unload();
integ_y->integrate();
}
derivs_x[0] = area; // volume = ∫ area dx
}
@ -40,9 +38,7 @@ double doubleIntegral( double x_start, double x_end, double y_start, double y_en
SA::RK2Integrator integ_x (dx, 1, state_x, state_x, deriv_x, &y_integ_context);
integ_x.setIndyVar(x_start);
while (integ_x.getIndyVar()<=x_end) {
integ_x.load();
integ_x.step();
integ_x.unload();
integ_x.integrate();
}
return volume;
}

121
trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh
View File

@ -7,12 +7,15 @@ namespace SA {
double X_out; // Independent Variable Out
double default_h; // Default step-size
void* user_data; // User data
void advanceIndyVar();
public:
Integrator(double h, void* udata): X_in(0.0), X_out(0.0), default_h(h), user_data(udata) {};
virtual ~Integrator() {}
virtual void step();
virtual void load();
virtual void unload();
virtual void step();
void integrate();
virtual double undo_integrate();
double getIndyVar();
void setIndyVar(double v);
};
@ -27,30 +30,47 @@ namespace SA {
double** inVars; // Array of pointers to the variables whose values comprise the input state vector.
double** outVars; // Array of pointers to the variables to which the output state vector values are to be disseminated.
DerivsFunc derivs_func; // Pointer to the function that calculates the derivatives to be integrated.
double last_h; // What was the last step-size? For undo_step().
// void advanceIndyVar(); // Inherited from Integrator::advanceIndyVar()
public:
FirstOrderODEIntegrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
virtual ~FirstOrderODEIntegrator();
virtual double undo_step();
void load();
void unload();
void load_from_outState();
double** set_in_vars( double* in_vars[]);
double** set_out_vars( double* out_vars[]);
void load(); // Overrides Integrator::load(), Final
void unload(); // Overrides Integrator::unload(), Final
virtual void step(); // Overrides Integrator::step(). Virtual
// integrate(); Inherited from Integrator::integrate(), Final
virtual double undo_integrate(); // Overrides Integrator::undo_integrate(), Virtual
void load_from_outState(); // New, Final
double** set_in_vars( double* in_vars[]); // New, Final
double** set_out_vars( double* out_vars[]);// New, Final
// double getIndyVar();
// void setIndyVar(double v);
};
typedef double (*RootErrorFunc)( double x, double state[], RootFinder* root_finder, void* udata);
class FirstOrderODEVariableStepIntegrator : public FirstOrderODEIntegrator {
private:
RootFinder* root_finder;
RootErrorFunc root_error_func;
protected:
double last_h;
// void advanceIndyVar(); // Inherited from Integrator::advanceIndyVar(), Final
void advanceIndyVar( double h);
public:
FirstOrderODEVariableStepIntegrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
~FirstOrderODEVariableStepIntegrator();
virtual void variable_step( double h);
void add_Rootfinder( RootFinder* root_finder, RootErrorFunc rfunc);
void step();
// load(); Inherited from FirstOrderODEIntegrator, Final
// unload(); Inherited from FirstOrderODEIntegrator, Final
void step(); // Overrides FirstOrderODEIntegrator::step(), Final
// integrate(); // Inherited from Integrator
double undo_integrate(); // Overrides FirstOrderODEIntegrator::undo_integrate(), Final
// load_from_outState(); // Inherited from FirstOrderODEIntegrator, Final
// set_in_vars(); // Inherited from FirstOrderODEIntegrator, Final
// set_out_vars(); // Inherited from FirstOrderODEIntegrator, Final
virtual void variable_step( double h); // New, Should be overridden in derived classes, Virtual
void add_Rootfinder( RootFinder* root_finder, RootErrorFunc rfunc); // New, Final
// double getIndyVar(); // Inherited from Integrator
// void setIndyVar(double v); // Inherited from Integrator
private:
void find_roots(double h, unsigned int depth);
};
@ -60,7 +80,19 @@ namespace SA {
double *derivs;
EulerIntegrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
~EulerIntegrator();
void variable_step( double h);
// load(); Inherited from FirstOrderODEIntegrator
// unload(); Inherited from FirstOrderODEIntegrator
// void step(); Inherited from FirstOrderODEVariableStepIntegrator
// void integrate() Inherited from FirstOrderODEVariableStepIntegrator
// double undo_integrate(); Inherited from FirstOrderODEVariableStepIntegrator
// load_from_outState(); Inherited from FirstOrderODEIntegrator
// set_in_vars(); Inherited from FirstOrderODEIntegrator
// set_out_vars(); Inherited from FirstOrderODEIntegrator
void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
// void integrate(double h) Inherited from FirstOrderODEVariableStepIntegrator
// void add_Rootfinder( RootFinder* root_finder, RootErrorFunc rfunc); Inherited from FirstOrderODEVariableStepIntegrator
// double getIndyVar(); Inherited from Integrator
// void setIndyVar(double v); Inherited from Integrator
};
class HeunsMethod : public FirstOrderODEVariableStepIntegrator {
@ -69,7 +101,7 @@ namespace SA {
double *derivs[2];
HeunsMethod( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
~HeunsMethod();
void variable_step( double h);
void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
};
class RK2Integrator : public FirstOrderODEVariableStepIntegrator {
@ -78,7 +110,7 @@ namespace SA {
double *derivs[2];
RK2Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc derivs_func, void* udata);
~RK2Integrator();
void variable_step( double h);
void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
};
class RK4Integrator : public FirstOrderODEVariableStepIntegrator {
@ -87,7 +119,7 @@ namespace SA {
double *derivs[4];
RK4Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
~RK4Integrator();
void variable_step( double h);
void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
};
class RK3_8Integrator : public FirstOrderODEVariableStepIntegrator {
@ -96,18 +128,18 @@ namespace SA {
double *derivs[4];
RK3_8Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
~RK3_8Integrator();
void variable_step( double h);
void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
};
class EulerCromerIntegrator : public Integrator {
protected:
unsigned int state_size;
double** x_p;
double** v_p;
double* x_in;
double* v_in;
double* x_out;
double* v_out;
unsigned int nDimensions;
double** pos_p;
double** vel_p;
double* pos_in;
double* vel_in;
double* pos_out;
double* vel_out;
double* g_out;
double* f_out;
DerivsFunc gderivs;
@ -120,25 +152,26 @@ namespace SA {
void step();
void load();
void unload();
double undo_integrate();
};
// class ABMIntegrator : public FirstOrderODEIntegrator {
// protected:
// double ** deriv_history;
// double * composite_deriv;
// unsigned int bufferSize;
// unsigned int algorithmHistorySize;
// unsigned int hix;
// unsigned int currentHistorySize;
//
// public:
// ABMIntegrator(int history_size, double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
// ~ABMIntegrator();
// void step();
// void undo_step();
// };
// AdamsBashforthMoulton Self Priming
class ABMIntegrator : public FirstOrderODEIntegrator {
protected:
double ** deriv_history;
double * composite_deriv;
unsigned int bufferSize;
unsigned int algorithmHistorySize;
unsigned int hix;
unsigned int currentHistorySize;
public:
ABMIntegrator(int history_size, double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
~ABMIntegrator();
void step();
double undo_integrate();
};
// AdamsBashforthMoulton 2
// AdamsBashforthMoulton 2 with RK2 Priming
class ABM2Integrator : public FirstOrderODEIntegrator {
protected:
double ** deriv_history;
@ -147,15 +180,14 @@ namespace SA {
unsigned int hix;
unsigned int currentHistorySize;
FirstOrderODEVariableStepIntegrator* priming_integrator;
public:
ABM2Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
~ABM2Integrator();
void step();
double undo_step();
double undo_integrate();
};
// AdamsBashforthMoulton 4
// AdamsBashforthMoulton 4 with RK4 Priming
class ABM4Integrator : public FirstOrderODEIntegrator {
protected:
double ** deriv_history;
@ -164,11 +196,10 @@ namespace SA {
unsigned int hix;
unsigned int currentHistorySize;
FirstOrderODEVariableStepIntegrator* priming_integrator;
public:
ABM4Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
~ABM4Integrator();
void step();
double undo_step();
double undo_integrate();
};
}

302
trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp
View File

@ -3,11 +3,33 @@
#include <iostream>
#include "SAIntegrator.hh"
void SA::Integrator::step() { X_out = X_in + default_h; }
void SA::Integrator::load() { X_in = X_out; }
void SA::Integrator::unload() {}
double SA::Integrator::getIndyVar() { return X_out; }
void SA::Integrator::setIndyVar(double v) { X_out = v; }
void SA::Integrator::advanceIndyVar() {
X_out = X_in + default_h;
}
void SA::Integrator::step() {
advanceIndyVar();
}
void SA::Integrator::load() {
X_in = X_out;
}
void SA::Integrator::unload() {
std::cout << "SA::Integrator : Nothing to unload." << std::endl;
}
void SA::Integrator::integrate() {
load();
step();
unload();
}
double SA::Integrator::undo_integrate() {
X_out = X_in;
return (default_h);
}
double SA::Integrator::getIndyVar() {
return X_out;
}
void SA::Integrator::setIndyVar(double v) {
X_out = v;
}
// ------------------------------------------------------------
SA::FirstOrderODEIntegrator::FirstOrderODEIntegrator(
@ -59,13 +81,24 @@ double** SA::FirstOrderODEIntegrator::set_out_vars( double* out_vars[]){
outVars = out_vars;
return (ret);
}
double SA::FirstOrderODEIntegrator::undo_step() {
for (unsigned int i=0 ; i<state_size; i++ ) {
*(inVars[i]) = inState[i];
}
X_out = X_in;
return (last_h);
double SA::FirstOrderODEIntegrator::undo_integrate() {
for (unsigned int i=0 ; i<state_size; i++ ) {
*(inVars[i]) = inState[i];
}
for (unsigned int i=0 ; i<state_size; i++ ) {
outState[i] = inState[i];
}
return (SA::Integrator::undo_integrate());
}
void SA::FirstOrderODEIntegrator::step() {
double derivs[state_size];
(*derivs_func)( X_in, inState, derivs, user_data);
for (unsigned int i=0; i<state_size; i++) {
outState[i] = inState[i] + derivs[i] * default_h;
}
advanceIndyVar();
}
// ------------------------------------------------------------
SA::FirstOrderODEVariableStepIntegrator::FirstOrderODEVariableStepIntegrator(
double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
@ -73,6 +106,10 @@ SA::FirstOrderODEVariableStepIntegrator::FirstOrderODEVariableStepIntegrator(
root_finder = (RootFinder*) NULL;
root_error_func = (RootErrorFunc)NULL;
}
double SA::FirstOrderODEVariableStepIntegrator::undo_integrate() {
SA::FirstOrderODEIntegrator::undo_integrate();
return (last_h);
}
SA::FirstOrderODEVariableStepIntegrator::~FirstOrderODEVariableStepIntegrator() {
// User is responsible for deleting the root_finder.
}
@ -81,15 +118,23 @@ void SA::FirstOrderODEVariableStepIntegrator::add_Rootfinder(
root_finder = rootFinder;
root_error_func = rfunc;
}
void SA::FirstOrderODEVariableStepIntegrator::variable_step( double h) {
void SA::FirstOrderODEVariableStepIntegrator::advanceIndyVar(double h) {
last_h = h; X_out = X_in + h;
}
void SA::FirstOrderODEVariableStepIntegrator::variable_step( double h) {
double derivs[state_size];
(*derivs_func)( X_in, inState, derivs, user_data);
for (unsigned int i=0; i<state_size; i++) {
outState[i] = inState[i] + derivs[i] * h;
}
advanceIndyVar(h);
}
void SA::FirstOrderODEVariableStepIntegrator::find_roots(double h, unsigned int depth) {
double new_h;
double h_correction = (*root_error_func)( getIndyVar(), outState, root_finder, user_data );
if (h_correction < 0.0) {
while (h_correction != 0.0) {
new_h = undo_step() + h_correction;
new_h = undo_integrate() + h_correction;
variable_step(new_h);
h_correction = (*root_error_func)( getIndyVar(), outState, root_finder, user_data );
}
@ -106,6 +151,7 @@ void SA::FirstOrderODEVariableStepIntegrator::step() {
find_roots( default_h, 5 );
}
}
// ------------------------------------------------------------
SA::EulerIntegrator::EulerIntegrator(
double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
@ -120,7 +166,7 @@ void SA::EulerIntegrator::variable_step( double h) {
for (unsigned int i=0; i<state_size; i++) {
outState[i] = inState[i] + derivs[i] * h;
}
SA::FirstOrderODEVariableStepIntegrator::variable_step(h);
advanceIndyVar(h);
}
// ------------------------------------------------------------
SA::HeunsMethod::HeunsMethod(
@ -146,7 +192,7 @@ void SA::HeunsMethod::variable_step( double h) {
for (unsigned int i=0; i<state_size; i++) {
outState[i] = inState[i] + (h/2) * ( derivs[0][i] + derivs[1][i] );
}
SA::FirstOrderODEVariableStepIntegrator::variable_step(h);
advanceIndyVar(h);
}
// ------------------------------------------------------------
SA::RK2Integrator::RK2Integrator(
@ -172,7 +218,7 @@ void SA::RK2Integrator::variable_step( double h) {
for (unsigned int i=0; i<state_size; i++) {
outState[i] = inState[i] + h * derivs[1][i];
}
SA::FirstOrderODEVariableStepIntegrator::variable_step(h);
advanceIndyVar(h);
}
// ------------------------------------------------------------
SA::RK4Integrator::RK4Integrator(
@ -213,7 +259,7 @@ void SA::RK4Integrator::variable_step( double h) {
(1/3.0)* derivs[2][i] +
(1/6.0)* derivs[3][i]) * h;
}
SA::FirstOrderODEVariableStepIntegrator::variable_step(h);
advanceIndyVar(h);
}
// ------------------------------------------------------------
SA::RK3_8Integrator::RK3_8Integrator(
@ -254,138 +300,146 @@ void SA::RK3_8Integrator::variable_step( double h) {
(3/8.0)* derivs[2][i] +
(1/8.0)* derivs[3][i]) * h;
}
SA::FirstOrderODEVariableStepIntegrator::variable_step(h);
advanceIndyVar(h);
}
// ------------------------------------------------------------
SA::EulerCromerIntegrator::EulerCromerIntegrator(
double dt, int N, double* xp[], double* vp[], DerivsFunc gfunc, DerivsFunc ffunc, void* udata)
:Integrator(dt, udata) {
state_size = N;
x_p = new double*[N];
v_p = new double*[N];
x_in = new double[N];
v_in = new double[N];
g_out = new double[N];
f_out = new double[N];
x_out = new double[N];
v_out = new double[N];
nDimensions = N;
pos_p = new double*[N];
pos_in = new double[N];
pos_out = new double[N];
vel_p = new double*[N];
vel_in = new double[N];
vel_out = new double[N];
g_out = new double[N];
f_out = new double[N];
for (int i=0 ; i<N; i++ ) {
x_p[i] = xp[i];
v_p[i] = vp[i];
pos_p[i] = xp[i];
vel_p[i] = vp[i];
}
gderivs = gfunc;
fderivs = ffunc;
}
SA::EulerCromerIntegrator::~EulerCromerIntegrator() {
delete x_p;
delete v_p;
delete x_in;
delete v_in;
delete pos_in;
delete pos_out;
delete vel_in;
delete vel_out;
delete g_out;
delete f_out;
delete x_out;
delete v_out;
}
void SA::EulerCromerIntegrator::step( double dt) {
(*gderivs)( X_in, x_in, g_out, user_data);
for (int i=0 ; i<state_size; i++ ) {
v_out[i] = v_in[i] + g_out[i] * dt;
(*gderivs)( X_in, pos_in, g_out, user_data);
for (int i=0 ; i<nDimensions; i++ ) {
vel_out[i] = vel_in[i] + g_out[i] * dt;
}
(*fderivs)( X_in, v_in, f_out, user_data);
for (int i=0 ; i<state_size; i++ ) {
x_out[i] = x_in[i] + f_out[i] * dt;
(*fderivs)( X_in, vel_in, f_out, user_data);
for (int i=0 ; i<nDimensions; i++ ) {
pos_out[i] = pos_in[i] + f_out[i] * dt;
}
X_out = X_in + dt;
}
void SA::EulerCromerIntegrator::step() {
step(default_h);
// if ((root_finder != NULL) && (root_error_func != NULL)) {
// find_roots( default_h, 5 );
// }
}
void SA::EulerCromerIntegrator::load() {
for (int i=0 ; i<state_size; i++ ) {
x_in[i] = *(x_p[i]);
v_in[i] = *(v_p[i]);
for (int i=0 ; i<nDimensions; i++ ) {
pos_in[i] = *(pos_p[i]);
vel_in[i] = *(vel_p[i]);
}
}
void SA::EulerCromerIntegrator::unload(){
for (int i=0 ; i<state_size; i++ ) {
*(x_p[i]) = x_out[i];
*(v_p[i]) = v_out[i];
for (int i=0 ; i<nDimensions; i++ ) {
*(pos_p[i]) = pos_out[i];
*(vel_p[i]) = vel_out[i];
}
}
double SA::EulerCromerIntegrator::undo_integrate() {
for (int i=0 ; i<nDimensions; i++ ) {
*(pos_p[i]) = pos_in[i];
*(vel_p[i]) = vel_in[i];
}
return (SA::Integrator::undo_integrate());
}
// // Adams-Bashforth Coefficients
// static const double ABCoeffs[5][5] = {
// {1.0, 0.0, 0.0, 0.0, 0.0},
// {(3/2.0), (-1/2.0), 0.0, 0.0, 0.0},
// {(23/12.0), (-16/12.0), (5/12.0), 0.0, 0.0},
// {(55/24.0), (-59/24.0), (37/24.0), (-9/24.0), 0.0},
// {(1901/720.0), (-2774/720.0), (2616/720.0), (-1274/720.0), (251/720.0)}
// };
static const double ABCoeffs[5][5] = {
{1.0, 0.0, 0.0, 0.0, 0.0},
{(3/2.0), (-1/2.0), 0.0, 0.0, 0.0},
{(23/12.0), (-16/12.0), (5/12.0), 0.0, 0.0},
{(55/24.0), (-59/24.0), (37/24.0), (-9/24.0), 0.0},
{(1901/720.0), (-2774/720.0), (2616/720.0), (-1274/720.0), (251/720.0)}
};
//
// // Adams-Moulton Coefficients
// static const double AMCoeffs[5][5] = {
// {1.0, 0.0, 0.0, 0.0, 0.0},
// {(1/2.0), (1/2.0), 0.0, 0.0, 0.0},
// {(5/12.0), (8/12.0), (-1/12.0), 0.0, 0.0},
// {(9/24.0), (19/24.0), (-5/24.0), (1/24.0), 0.0},
// {(251/720.0), (646/720.0), (-264/720.0), (106/720.0), (-19/720.0)}
// };
// SA::ABMIntegrator::ABMIntegrator ( int history_size, double h, int N, double* in_vars[], double* out_vars[], DerivsFunc func, void* udata)
// : FirstOrderODEIntegrator(h, N, in_vars, out_vars ,func, udata) {
//
// algorithmHistorySize=history_size; // The amount of history that we intend to use in this ABM integrator.
// bufferSize=algorithmHistorySize+1; // The size of the buffer needs to be one more than the history so that an integration step can be reset.
// hix=0;
// currentHistorySize=0; // How much derivative history is currently stored int the history buffer. Initially there will be none until we've primed the integrator.
// deriv_history = new double*[bufferSize];
// for (unsigned int i=0; i<bufferSize ; i++) {
// deriv_history[i] = new double[state_size];
// }
// composite_deriv = new double[state_size];
// }
// SA::ABMIntegrator::~ABMIntegrator() {
// for (int i=0; i<bufferSize ; i++) {
// delete (deriv_history[i]);
// }
// delete(deriv_history);
// delete(composite_deriv);
// }
// void SA::ABMIntegrator::step() {
//
// hix = (hix+1)%bufferSize; // Move history index forward
// // (*derivs_func)( indyVar, inState, deriv_history[hix], user_data);
// (*derivs_func)( X_in, inState, deriv_history[hix], user_data); // Calculated and store the deriv for current, corrected state.
// // Increase the size of the stored history, up to the limit specified by the user.
// currentHistorySize = (currentHistorySize < algorithmHistorySize) ? currentHistorySize+1 : algorithmHistorySize;
// // (Predictor) Predict the next state using the Adams-Bashforth explicit method.
// for (int i=0; i<state_size; i++) {
// composite_deriv[i] = 0.0;
// for (int n=0,j=hix; n<currentHistorySize ; n++,j=(j+bufferSize-1)%bufferSize) {
// composite_deriv[i] += ABCoeffs[currentHistorySize-1][n] * deriv_history[j][i];
// }
// outState[i] = inState[i] + default_h * composite_deriv[i];
// }
// // Move history index forward, so we can temporarily store the derivative of the predicted next state.
// // We do not increase the currentHistorySize.
// hix = (hix+1)%bufferSize;
// (*derivs_func)( X_out, outState, deriv_history[hix], user_data); // Calc deriv for predicted next state.
// // (Corrector) Refine the next state using the Adams-Moulton implicit method. This is the corrected next state.
// for (int i=0; i<state_size; i++) {
// composite_deriv[i] = 0.0;
// for (int n=0,j=hix; n<currentHistorySize ; n++,j=(j+bufferSize-1)%bufferSize) {
// composite_deriv[i] += AMCoeffs[currentHistorySize-1][n] * deriv_history[j][i];
// }
// outState[i] = inState[i] + default_h * composite_deriv[i];
// }
// // Move history index backward, so we over-write the predicted state with the corrected state on our next step().
// hix = (hix+bufferSize-1)%bufferSize;
// SA::Integrator::step();
// }
//
// double SA::ABMIntegrator::undo_step() {
// hix = (hix+bufferSize-1)%bufferSize;
// return (FirstOrderODEIntegrator::undo_step());
// }
static const double AMCoeffs[5][5] = {
{1.0, 0.0, 0.0, 0.0, 0.0},
{(1/2.0), (1/2.0), 0.0, 0.0, 0.0},
{(5/12.0), (8/12.0), (-1/12.0), 0.0, 0.0},
{(9/24.0), (19/24.0), (-5/24.0), (1/24.0), 0.0},
{(251/720.0), (646/720.0), (-264/720.0), (106/720.0), (-19/720.0)}
};
SA::ABMIntegrator::ABMIntegrator ( int history_size, double h, int N, double* in_vars[], double* out_vars[], DerivsFunc func, void* udata)
: FirstOrderODEIntegrator(h, N, in_vars, out_vars ,func, udata) {
algorithmHistorySize=history_size; // The amount of history that we intend to use in this ABM integrator.
bufferSize=algorithmHistorySize+1; // The size of the buffer needs to be one more than the history so that an integration step can be reset.
hix=0;
currentHistorySize=0; // How much derivative history is currently stored int the history buffer. Initially there will be none until we've primed the integrator.
deriv_history = new double*[bufferSize];
for (unsigned int i=0; i<bufferSize ; i++) {
deriv_history[i] = new double[state_size];
}
composite_deriv = new double[state_size];
}
SA::ABMIntegrator::~ABMIntegrator() {
for (int i=0; i<bufferSize ; i++) {
delete (deriv_history[i]);
}
delete(deriv_history);
delete(composite_deriv);
}
void SA::ABMIntegrator::step() {
hix = (hix+1)%bufferSize; // Move history index forward
(*derivs_func)( X_in, inState, deriv_history[hix], user_data); // Calculated and store the deriv for current, corrected state.
// Increase the size of the stored history, up to the limit specified by the user.
currentHistorySize = (currentHistorySize < algorithmHistorySize) ? currentHistorySize+1 : algorithmHistorySize;
// (Predictor) Predict the next state using the Adams-Bashforth explicit method.
for (int i=0; i<state_size; i++) {
composite_deriv[i] = 0.0;
for (int n=0,j=hix; n<currentHistorySize ; n++,j=(j+bufferSize-1)%bufferSize) {
composite_deriv[i] += ABCoeffs[currentHistorySize-1][n] * deriv_history[j][i];
}
outState[i] = inState[i] + default_h * composite_deriv[i];
}
// Move history index forward, so we can temporarily store the derivative of the predicted next state.
// We do not increase the currentHistorySize.
hix = (hix+1)%bufferSize;
(*derivs_func)( X_out, outState, deriv_history[hix], user_data); // Calc deriv for predicted next state.
// (Corrector) Refine the next state using the Adams-Moulton implicit method. This is the corrected next state.
for (int i=0; i<state_size; i++) {
composite_deriv[i] = 0.0;
for (int n=0,j=hix; n<currentHistorySize ; n++,j=(j+bufferSize-1)%bufferSize) {
composite_deriv[i] += AMCoeffs[currentHistorySize-1][n] * deriv_history[j][i];
}
outState[i] = inState[i] + default_h * composite_deriv[i];
}
// Move history index backward, so we over-write the predicted state with the corrected state on our next step().
hix = (hix+bufferSize-1)%bufferSize;
SA::Integrator::step();
}
double SA::ABMIntegrator::undo_integrate() {
hix = (hix+bufferSize-1)%bufferSize;
return (FirstOrderODEIntegrator::undo_integrate());
}
// =======================================================================
static const double AB2Coeffs[2] = {(3/2.0), (-1/2.0)};
@ -404,7 +458,7 @@ SA::ABM2Integrator::ABM2Integrator (
}
composite_deriv = new double[state_size];
// NEW: Create a priming integrator.
// Create a priming integrator.
double** priming_integrator_in_vars = new double*[state_size];
for (unsigned int i=0; i<N ; i++) {
priming_integrator_in_vars[i] = &(inState[i]);
@ -458,11 +512,11 @@ void SA::ABM2Integrator::step() {
// Move history index backward, so we over-write the predicted state with the corrected state on our next step().
hix = (hix+bufferSize-1)%bufferSize;
}
SA::Integrator::step();
advanceIndyVar();
}
double SA::ABM2Integrator::undo_step() {
double SA::ABM2Integrator::undo_integrate() {
hix = (hix+bufferSize-1)%bufferSize;
return ( FirstOrderODEIntegrator::undo_step());
return ( FirstOrderODEIntegrator::undo_integrate());
}
// =======================================================================
static const double AB4Coeffs[4] = {(55/24.0), (-59/24.0), (37/24.0), (-9/24.0)};
@ -534,9 +588,9 @@ void SA::ABM4Integrator::step() {
// Move history index backward, so we over-write the predicted state with the corrected state on our next step().
hix = (hix+bufferSize-1)%bufferSize;
}
SA::Integrator::step();
advanceIndyVar();
}
double SA::ABM4Integrator::undo_step() {
double SA::ABM4Integrator::undo_integrate() {
hix = (hix+bufferSize-1)%bufferSize;
return (FirstOrderODEIntegrator::undo_step());
return (FirstOrderODEIntegrator::undo_integrate());
}

159
trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEIntegrator_unittest.cc Normal file
View File

@ -0,0 +1,159 @@
#include <gtest/gtest.h>
#include <iostream>
#include "SAIntegrator.hh"
#include <math.h>
#define EXCEPTABLE_ERROR 0.00000000001
void calc_derivs( double x __attribute__((unused)),
double state[] __attribute__((unused)),
double derivs[],
void* udata __attribute__((unused))) {
derivs[0] = 1.0;
derivs[1] = 1.0;
derivs[2] = 1.0;
derivs[3] = 1.0;
}
/*
The default integration algorithm for FirstOrderODEIntegrator is Euler integration.
Tests:
* getIndyVar
* integrate
* load
* step
* unload
*/
TEST(FirstOrderODEIntegrator_unittest, test_1) {
double h = 0.1;
unsigned int N = 4;
double vars[4] = {1.0, 2.0, 3.0, 4.0} ;
double* varptrs[4] = {&vars[0], &vars[1], &vars[2], &vars[3]};
double x;
SA::FirstOrderODEIntegrator integ( h, N, varptrs, varptrs, calc_derivs, NULL);
x = integ.getIndyVar();
EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
integ.integrate();
EXPECT_NEAR(vars[0], 1.1, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[1], 2.1, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[2], 3.1, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[3], 4.1, EXCEPTABLE_ERROR);
x = integ.getIndyVar();
EXPECT_NEAR(x, 0.1, EXCEPTABLE_ERROR);
}
/*
Tests:
* getIndyVar
* integrate
* load
* step
* unload
* undo_integrate
*/
TEST(FirstOrderODEIntegrator_unittest, test_2) {
double h = 0.1;
unsigned int N = 4;
double vars[4] = {1.0, 2.0, 3.0, 4.0} ;
double* varptrs[4] = {&vars[0], &vars[1], &vars[2], &vars[3]};
double x;
SA::FirstOrderODEIntegrator integ( h, N, varptrs, varptrs, calc_derivs, NULL);
integ.integrate();
EXPECT_NEAR(vars[0], 1.1, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[1], 2.1, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[2], 3.1, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[3], 4.1, EXCEPTABLE_ERROR);
x = integ.getIndyVar();
EXPECT_NEAR(x, 0.1, EXCEPTABLE_ERROR);
integ.undo_integrate();
EXPECT_NEAR(vars[0], 1.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[1], 2.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[2], 3.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[3], 4.0, EXCEPTABLE_ERROR);
x = integ.getIndyVar();
EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
}
/*
Tests:
* getIndyVar
* integrate
* load
* step
* unload
* load_from_outState
*/
TEST(FirstOrderODEIntegrator_unittest, test_3) {
double h = 0.1;
unsigned int N = 4;
double vars[4] = {1.0, 2.0, 3.0, 4.0};
double* varptrs[4] = {&vars[0], &vars[1], &vars[2], &vars[3]};
SA::FirstOrderODEIntegrator integ( h, N, varptrs, varptrs, calc_derivs, NULL);
integ.integrate();
for (int i=0; i<5 ; i++) {
integ.load_from_outState();
integ.step();
}
integ.unload();
EXPECT_NEAR(vars[0], 1.6, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[1], 2.6, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[2], 3.6, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[3], 4.6, EXCEPTABLE_ERROR);
double x = integ.getIndyVar();
EXPECT_NEAR(x, 0.6, EXCEPTABLE_ERROR);
}
/*
Tests:
* getIndyVar
* integrate
* load
* step
* unload
* load_from_outState
* set_out_vars
*/
TEST(FirstOrderODEIntegrator_unittest, test_4) {
double h = 0.1;
unsigned int N = 4;
double vars[4] = {1.0, 2.0, 3.0, 4.0};
double* varptrs[4] = {&vars[0], &vars[1], &vars[2], &vars[3]};
double altvars[4];
double* altvarptrs[4] = {&altvars[0], &altvars[1], &altvars[2], &altvars[3]};
SA::FirstOrderODEIntegrator integ( h, N, varptrs, varptrs, calc_derivs, NULL);
integ.integrate();
for (int i=0; i<5 ; i++) {
integ.load_from_outState();
integ.step();
}
integ.set_out_vars( altvarptrs);
integ.unload();
EXPECT_NEAR(altvars[0], 1.6, EXCEPTABLE_ERROR);
EXPECT_NEAR(altvars[1], 2.6, EXCEPTABLE_ERROR);
EXPECT_NEAR(altvars[2], 3.6, EXCEPTABLE_ERROR);
EXPECT_NEAR(altvars[3], 4.6, EXCEPTABLE_ERROR);
integ.set_out_vars( varptrs);
integ.unload();
EXPECT_NEAR(vars[0], 1.6, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[1], 2.6, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[2], 3.6, EXCEPTABLE_ERROR);
EXPECT_NEAR(vars[3], 4.6, EXCEPTABLE_ERROR);
double x = integ.getIndyVar();
EXPECT_NEAR(x, 0.6, EXCEPTABLE_ERROR);
}

172
trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEVariableStepIntegrator_unittest.cc Normal file
View File

@ -0,0 +1,172 @@
#include <gtest/gtest.h>
#include <iostream>
#include "SAIntegrator.hh"
#include <math.h>
void deriv4( double t,
double state[] __attribute__((unused)),
double derivs[],
void* udata __attribute__((unused))) {
double t2 = t*t;
double t3 = t2*t;
derivs[0] = 3.0*t3 - 3.0*t2 - 3.0*t + 4.0;
derivs[1] = 2.0*t3 + 1.0*t2 - 5.0*t + 7.0;
derivs[2] = t3 + t2 + 5.0*t + 4.0;
derivs[3] = t3 - 6.0*t2 + 4.0*t + 12.0;
}
void deriv2( double t,
double state[] __attribute__((unused)),
double derivs[],
void* udata __attribute__((unused))) {
derivs[0] = -2.0*t + 3.0;
derivs[1] = 5.0*t + 4.0;
derivs[2] = 3.0*t - 3.0;
derivs[3] = 5.0*t + 4.0;
}
void deriv1( double t __attribute__((unused)),
double state[] __attribute__((unused)),
double derivs[],
void* udata __attribute__((unused))) {
derivs[0] = 4.0;
derivs[1] = -4.0;
derivs[2] = M_PI;
derivs[3] = -M_PI;
}
#define EXCEPTABLE_ERROR 0.00000000001
TEST(FirstOrderODEVariableStepIntegrator_unittest, EulerIntegrator_1) {
double state[4] = {0.0, 0.0, 0.0, 0.0};
double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3])};
double dt = 0.01;
unsigned int count = 0;
double t = count * dt;
SA::EulerIntegrator integ(dt, 4, state_var_p, state_var_p, deriv1, NULL);
while (t < 2.0) {
integ.load();
integ.step();
integ.unload();
count ++;
t = count * dt;
}
EXPECT_NEAR(state[0], 8.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], -8.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], (2*M_PI), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], (-2*M_PI), EXCEPTABLE_ERROR);
}
TEST(FirstOrderODEVariableStepIntegrator_unittest, RungeKutta2_1) {
double state[4] = {0.0, 0.0, 0.0, 0.0};
double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3]) };
double dt = 0.01;
unsigned int count = 0;
double t = count * dt;
SA::RK2Integrator integ(dt, 4, state_var_p, state_var_p, deriv2, NULL);
while (t < 2.0) {
integ.load();
integ.step();
integ.unload();
count ++;
t = count * dt;
}
EXPECT_NEAR(state[0], 2.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], 18.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], 18.0, EXCEPTABLE_ERROR);
}
TEST(FirstOrderODEVariableStepIntegrator_unittest, RungeKutta4_1) {
double state[4] = {0.0, 0.0, 0.0, 0.0};
double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3]) };
double dt = 0.01;
unsigned int count = 0;
double t = count * dt;
SA::RK4Integrator integ(dt, 4, state_var_p, state_var_p, deriv4, NULL);
while (t < 2.0) {
integ.load();
integ.step();
integ.unload();
count ++;
t = count * dt;
}
EXPECT_NEAR(state[0], 6.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], (14.0+2/3.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], (24+2/3.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], 20.0, EXCEPTABLE_ERROR);
}
TEST(FirstOrderODEVariableStepIntegrator_unittest, RungeKutta38_1) {
double state[4] = {0.0, 0.0, 0.0, 0.0};
double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3]) };
double dt = 0.01;
unsigned int count = 0;
double t = count * dt;
SA::RK3_8Integrator integ(dt, 4, state_var_p, state_var_p, deriv4, NULL);
while (t < 2.0) {
integ.load();
integ.step();
integ.unload();
count ++;
t = count * dt;
}
EXPECT_NEAR(state[0], 6.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], (14.0+2/3.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], (24+2/3.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], 20.0, EXCEPTABLE_ERROR);
}
TEST(FirstOrderODEVariableStepIntegrator_unittest, EulerIntegrator_undo_integrate) {
double state[4] = {0.0, 0.0, 0.0, 0.0};
double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3])};
double dx = 0.1;
double x;
SA::EulerIntegrator integ(dx, 4, state_var_p, state_var_p, deriv1, NULL);
integ.load();
integ.step();
integ.unload();
x = integ.getIndyVar();
EXPECT_NEAR(state[0], (4.0/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], (-4.0/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], (M_PI/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], (-M_PI/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(x, 0.1, EXCEPTABLE_ERROR);
integ.undo_integrate();
x = integ.getIndyVar();
EXPECT_NEAR(state[0], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
integ.undo_integrate();
x = integ.getIndyVar();
EXPECT_NEAR(state[0], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
integ.load();
integ.step();
integ.unload();
x = integ.getIndyVar();
EXPECT_NEAR(state[0], (4.0/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], (-4.0/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], (M_PI/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], (-M_PI/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(x, 0.1, EXCEPTABLE_ERROR);
}

18
trick_source/trick_utils/SAIntegrator/unittest/Makefile
View File

@ -15,8 +15,10 @@ LIBDIRS = -L${SAI_LIBDIR} -L${GTEST_HOME}/lib64 -L${GTEST_HOME}/lib
all: test
test: SAIntegrator_unittest RootFinder_unittest
test: SAIntegrator_unittest FirstOrderODEIntegrator_unittest FirstOrderODEVariableStepIntegrator_unittest RootFinder_unittest
./SAIntegrator_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/SAIntegrator_unittest.xml
./FirstOrderODEIntegrator_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/FirstOrderODEIntegrator_unittest.xml
./FirstOrderODEVariableStepIntegrator_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/FirstOrderODEVariableStepIntegrator_unittest.xml
./RootFinder_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/RootFinder_unittest.xml
SAIntegrator_unittest.o : SAIntegrator_unittest.cc
@ -25,6 +27,18 @@ SAIntegrator_unittest.o : SAIntegrator_unittest.cc
SAIntegrator_unittest : ${SAI_LIBDIR}/${SAI_LIBNAME} SAIntegrator_unittest.o
$(TRICK_CXX) $(TRICK_CPPFLAGS) -o $@ $^ ${LIBDIRS} -lSAInteg -lgtest -lgtest_main -lpthread
FirstOrderODEIntegrator_unittest.o : FirstOrderODEIntegrator_unittest.cc
$(TRICK_CXX) $(TRICK_CPPFLAGS) $(INCLUDE_DIRS) -c $<
FirstOrderODEIntegrator_unittest : ${SAI_LIBDIR}/${SAI_LIBNAME} FirstOrderODEIntegrator_unittest.o
$(TRICK_CXX) $(TRICK_CPPFLAGS) -o $@ $^ ${LIBDIRS} -lSAInteg -lgtest -lgtest_main -lpthread
FirstOrderODEVariableStepIntegrator_unittest.o : FirstOrderODEVariableStepIntegrator_unittest.cc
$(TRICK_CXX) $(TRICK_CPPFLAGS) $(INCLUDE_DIRS) -c $<
FirstOrderODEVariableStepIntegrator_unittest : ${SAI_LIBDIR}/${SAI_LIBNAME} FirstOrderODEVariableStepIntegrator_unittest.o
$(TRICK_CXX) $(TRICK_CPPFLAGS) -o $@ $^ ${LIBDIRS} -lSAInteg -lgtest -lgtest_main -lpthread
RootFinder_unittest.o : RootFinder_unittest.cc
$(TRICK_CXX) $(TRICK_CPPFLAGS) $(INCLUDE_DIRS) -c $<
@ -38,4 +52,4 @@ clean:
${RM} *.o
spotless: clean
${RM} SAIntegrator_unittest
${RM} FirstOrderODEVariableStepIntegrator_unittest

216
trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc
View File

@ -3,170 +3,72 @@
#include "SAIntegrator.hh"
#include <math.h>
void deriv4( double t,
double state[] __attribute__((unused)),
double derivs[],
void* udata __attribute__((unused))) {
double t2 = t*t;
double t3 = t2*t;
derivs[0] = 3.0*t3 - 3.0*t2 - 3.0*t + 4.0;
derivs[1] = 2.0*t3 + 1.0*t2 - 5.0*t + 7.0;
derivs[2] = t3 + t2 + 5.0*t + 4.0;
derivs[3] = t3 - 6.0*t2 + 4.0*t + 12.0;
}
void deriv2( double t,
double state[] __attribute__((unused)),
double derivs[],
void* udata __attribute__((unused))) {
derivs[0] = -2.0*t + 3.0;
derivs[1] = 5.0*t + 4.0;
derivs[2] = 3.0*t - 3.0;
derivs[3] = 5.0*t + 4.0;
}
void deriv1( double t __attribute__((unused)),
double state[] __attribute__((unused)),
double derivs[],
void* udata __attribute__((unused))) {
derivs[0] = 4.0;
derivs[1] = -4.0;
derivs[2] = M_PI;
derivs[3] = -M_PI;
}
#define EXCEPTABLE_ERROR 0.00000000001
TEST(SAIntegrator_unittest, EulerIntegrator_1) {
TEST(Integrator, test_getIndyVar_and_setIndyVar) {
double h = 0.1;
void* udata = NULL;
double x;
double state[4] = {0.0, 0.0, 0.0, 0.0};
double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3])};
double dt = 0.01;
unsigned int count = 0;
double t = count * dt;
SA::EulerIntegrator integ(dt, 4, state_var_p, state_var_p, deriv1, NULL);
while (t < 2.0) {
integ.load();
integ.step();
integ.unload();
count ++;
t = count * dt;
}
EXPECT_NEAR(state[0], 8.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], -8.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], (2*M_PI), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], (-2*M_PI), EXCEPTABLE_ERROR);
SA::Integrator integ( h, udata);
//Confirm that initial value of the independent variable is 0.0;
x = integ.getIndyVar();
EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
//Set the independent variable, and then confirm that it's the value we set.;
integ.setIndyVar(3.0);
x = integ.getIndyVar();
EXPECT_NEAR(x, 3.0, EXCEPTABLE_ERROR);
}
TEST(Integrator, test_load_and_step) {
double h = 0.1;
void* udata = NULL;
double x;
SA::Integrator integ( h, udata);
integ.setIndyVar(4.0);
integ.load();
integ.step();
x = integ.getIndyVar();
EXPECT_NEAR(x, 4.1, EXCEPTABLE_ERROR);
integ.load();
integ.step();
x = integ.getIndyVar();
EXPECT_NEAR(x, 4.2, EXCEPTABLE_ERROR);
}
TEST(SAIntegrator_unittest, RungeKutta2_1) {
TEST(Integrator, test_integrate) {
double h = 0.1;
void* udata = NULL;
double x;
double state[4] = {0.0, 0.0, 0.0, 0.0};
double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3]) };
double dt = 0.01;
unsigned int count = 0;
double t = count * dt;
SA::RK2Integrator integ(dt, 4, state_var_p, state_var_p, deriv2, NULL);
while (t < 2.0) {
integ.load();
integ.step();
integ.unload();
count ++;
t = count * dt;
}
EXPECT_NEAR(state[0], 2.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], 18.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], 18.0, EXCEPTABLE_ERROR);
SA::Integrator integ( h, udata);
integ.setIndyVar(4.0);
integ.integrate();
x = integ.getIndyVar();
EXPECT_NEAR(x, 4.1, EXCEPTABLE_ERROR);
integ.integrate();
x = integ.getIndyVar();
EXPECT_NEAR(x, 4.2, EXCEPTABLE_ERROR);
}
TEST(SAIntegrator_unittest, RungeKutta4_1) {
TEST(Integrator, test_undo_integrate) {
double h = 0.1;
void* udata = NULL;
double x;
double state[4] = {0.0, 0.0, 0.0, 0.0};
double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3]) };
double dt = 0.01;
unsigned int count = 0;
double t = count * dt;
SA::RK4Integrator integ(dt, 4, state_var_p, state_var_p, deriv4, NULL);
while (t < 2.0) {
integ.load();
integ.step();
integ.unload();
count ++;
t = count * dt;
}
EXPECT_NEAR(state[0], 6.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], (14.0+2/3.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], (24+2/3.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], 20.0, EXCEPTABLE_ERROR);
}
TEST(SAIntegrator_unittest, RungeKutta38_1) {
double state[4] = {0.0, 0.0, 0.0, 0.0};
double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3]) };
double dt = 0.01;
unsigned int count = 0;
double t = count * dt;
SA::RK3_8Integrator integ(dt, 4, state_var_p, state_var_p, deriv4, NULL);
while (t < 2.0) {
integ.load();
integ.step();
integ.unload();
count ++;
t = count * dt;
}
EXPECT_NEAR(state[0], 6.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], (14.0+2/3.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], (24+2/3.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], 20.0, EXCEPTABLE_ERROR);
}
TEST(SAIntegrator_unittest, EulerIntegrator_undo_step) {
double state[4] = {0.0, 0.0, 0.0, 0.0};
double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3])};
double dx = 0.1;
double x;
SA::EulerIntegrator integ(dx, 4, state_var_p, state_var_p, deriv1, NULL);
integ.load();
integ.step();
integ.unload();
x = integ.getIndyVar();
EXPECT_NEAR(state[0], (4.0/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], (-4.0/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], (M_PI/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], (-M_PI/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(x, 0.1, EXCEPTABLE_ERROR);
integ.undo_step();
x = integ.getIndyVar();
EXPECT_NEAR(state[0], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
integ.undo_step();
x = integ.getIndyVar();
EXPECT_NEAR(state[0], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], 0.0, EXCEPTABLE_ERROR);
EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
integ.load();
integ.step();
integ.unload();
x = integ.getIndyVar();
EXPECT_NEAR(state[0], (4.0/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[1], (-4.0/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[2], (M_PI/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(state[3], (-M_PI/10.0), EXCEPTABLE_ERROR);
EXPECT_NEAR(x, 0.1, EXCEPTABLE_ERROR);
SA::Integrator integ( h, udata);
integ.setIndyVar(4.0);
integ.integrate();
x = integ.getIndyVar();
EXPECT_NEAR(x, 4.1, EXCEPTABLE_ERROR);
integ.undo_integrate();
x = integ.getIndyVar();
EXPECT_NEAR(x, 4.0, EXCEPTABLE_ERROR);
}