ExternalVendorCode/trick
1
0
Fork 0
mirror of https://github.com/nasa/trick.git synced 2025-04-14 14:36:46 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Refactoring: eliminate assumption that independent variable is time. #1058

Browse Source
This commit is contained in:
Penn, John M 047828115 2020-10-05 14:16:50 -05:00
parent fc83dbe6f3
commit 5884e64920
4 changed files with 125 additions and 116 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -21,8 +21,8 @@ void deriv_x( double x, double state_x[], double derivs_x[], void* udata) {
SA::RK2Integrator* integ_y = context->integ;
integ_y->set_in_vars(state_y);
integ_y->set_out_vars(state_y);
integ_y->setTime(context->start);
while (integ_y->getTime() <= context->end) {
integ_y->setIndyVar(context->start);
while (integ_y->getIndyVar() <= context->end) {
integ_y->load();
integ_y->step();
integ_y->unload();
@ -38,8 +38,8 @@ double doubleIntegral( double x_start, double x_end, double y_start, double y_en
SA::RK2Integrator integ_y (dy, 2, NULL, NULL, deriv_y, NULL);
IContext y_integ_context = {&integ_y, y_start, y_end};
SA::RK2Integrator integ_x (dx, 1, state_x, state_x, deriv_x, &y_integ_context);
integ_x.setTime(x_start);
while (integ_x.getTime()<=x_end) { // x-end
integ_x.setIndyVar(x_start);
while (integ_x.getIndyVar()<=x_end) {
integ_x.load();
integ_x.step();
integ_x.unload();

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

@ -3,20 +3,22 @@ namespace SA {
class Integrator {
protected:
double time;
double default_dt;
void * udata;
double indyVar; // Independent Variable
double default_h; // Default step-size
void * udata; // User data
public:
Integrator(double dt, void* user_data): time(0.0), default_dt(dt), udata(user_data) {};
Integrator(double h, void* user_data): indyVar(0.0), default_h(h), udata(user_data) {};
virtual ~Integrator() {}
virtual void step() { time += default_dt; }; // Integrate over implicit (default) time step (default: default_dt)
virtual void step() { indyVar += default_h; }; // Integrate over default step-size (default_h)
virtual void load() = 0;
virtual void unload() = 0;
double getTime() {return time;}
void setTime(double t) {time = t;}
double getIndyVar() {return indyVar;}
void setIndyVar(double t) {indyVar = t;}
double getTime() {return indyVar;}
void setTime(double t) {indyVar = t;}
};
typedef void (*derivsFunc)( double t, double state[], double derivs[], void* user_data);
typedef void (*derivsFunc)( double x, double state[], double derivs[], void* user_data);
class FirstOrderODEIntegrator : public Integrator {
protected:
@ -28,10 +30,10 @@ namespace SA {
derivsFunc derivs_func;
bool reset;
double last_dt; // What was the last step-size? For undo_step().
double last_h; // What was the last step-size? For undo_step().
public:
FirstOrderODEIntegrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data);
FirstOrderODEIntegrator( double h, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data);
~FirstOrderODEIntegrator();
virtual void undo_step();
void load();
@ -43,17 +45,17 @@ namespace SA {
class FirstOrderODEVariableStepIntegrator : public FirstOrderODEIntegrator {
public:
FirstOrderODEVariableStepIntegrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data);
FirstOrderODEVariableStepIntegrator( double h, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data);
using Integrator::step;
virtual void step( double dt );
virtual void step( double h );
};
class EulerIntegrator : public FirstOrderODEVariableStepIntegrator {
public:
double *derivs;
EulerIntegrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
EulerIntegrator( double h, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
~EulerIntegrator();
void step( double dt);
void step( double h);
void step();
};
@ -62,9 +64,9 @@ namespace SA {
double *wstate;
double *derivs[2];
HeunsMethod( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
HeunsMethod( double h, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
~HeunsMethod();
void step( double dt);
void step( double h);
void step();
};
@ -73,9 +75,9 @@ namespace SA {
double *wstate;
double *derivs[2];
RK2Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
RK2Integrator( double h, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
~RK2Integrator();
void step( double dt);
void step( double h);
void step();
};
@ -84,9 +86,9 @@ namespace SA {
double *wstate[3];
double *derivs[4];
RK4Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
RK4Integrator( double h, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
~RK4Integrator();
void step( double dt);
void step( double h);
void step();
};
@ -96,9 +98,9 @@ namespace SA {
double *wstate[4];
double *derivs[4];
RK3_8Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
RK3_8Integrator( double h, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
~RK3_8Integrator();
void step( double dt);
void step( double h);
void step();
};
@ -137,7 +139,7 @@ namespace SA {
// unsigned int currentHistorySize;
//
// public:
// ABMIntegrator(int history_size, double dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
// ABMIntegrator(int history_size, double h, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
// ~ABMIntegrator();
// void step();
// void undo_step();
@ -154,7 +156,7 @@ namespace SA {
FirstOrderODEVariableStepIntegrator* priming_integrator;
public:
ABM2Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
ABM2Integrator( double h, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
~ABM2Integrator();
void step();
void undo_step();
@ -171,7 +173,7 @@ namespace SA {
FirstOrderODEVariableStepIntegrator* priming_integrator;
public:
ABM4Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
ABM4Integrator( double h, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data);
~ABM4Integrator();
void step();
void undo_step();

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

@ -5,8 +5,8 @@
// ------------------------------------------------------------
SA::FirstOrderODEIntegrator::FirstOrderODEIntegrator(double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: Integrator(dt, user_data) {
SA::FirstOrderODEIntegrator::FirstOrderODEIntegrator(double h, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: Integrator(h, user_data) {
state_size = N;
orig_vars = in_vars;
dest_vars = out_vars;
@ -49,45 +49,45 @@ void SA::FirstOrderODEIntegrator::undo_step() {
for (int i=0 ; i<state_size; i++ ) {
*(orig_vars[i]) = istate[i]; // Copy istate values back to the state variables from whence they came.
}
time -= last_dt; // Undo the last time-step.
indyVar -= last_h; // Undo the last step.
}
reset = true;
}
// ------------------------------------------------------------
SA::FirstOrderODEVariableStepIntegrator::FirstOrderODEVariableStepIntegrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEIntegrator(dt, N, in_vars, out_vars, func, user_data) {}
SA::FirstOrderODEVariableStepIntegrator::FirstOrderODEVariableStepIntegrator( double h, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEIntegrator(h, N, in_vars, out_vars, func, user_data) {}
void SA::FirstOrderODEVariableStepIntegrator::step( double dt ) {
last_dt = dt; time += dt;
void SA::FirstOrderODEVariableStepIntegrator::step( double h ) {
last_h = h; indyVar += h;
}
// ------------------------------------------------------------
SA::EulerIntegrator::EulerIntegrator(double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEVariableStepIntegrator(dt, N, in_vars, out_vars, func, user_data)
SA::EulerIntegrator::EulerIntegrator(double h, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEVariableStepIntegrator(h, N, in_vars, out_vars, func, user_data)
{
derivs = new double[N];
}
SA::EulerIntegrator::~EulerIntegrator() {
delete derivs;
}
void SA::EulerIntegrator::step( double dt) {
(*derivs_func)( time, istate, derivs, udata);
void SA::EulerIntegrator::step( double h) {
(*derivs_func)( indyVar, istate, derivs, udata);
for (unsigned int i=0; i<state_size; i++) {
ostate[i] = istate[i] + derivs[i] * dt;
ostate[i] = istate[i] + derivs[i] * h;
}
SA::FirstOrderODEVariableStepIntegrator::step(dt);
SA::FirstOrderODEVariableStepIntegrator::step(h);
}
void SA::EulerIntegrator::step() {
step(default_dt);
step(default_h);
}
// ------------------------------------------------------------
SA::HeunsMethod::HeunsMethod( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEVariableStepIntegrator(dt, N, in_vars, out_vars, func, user_data)
SA::HeunsMethod::HeunsMethod( double h, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEVariableStepIntegrator(h, N, in_vars, out_vars, func, user_data)
{
wstate = new double[N];
for (unsigned int i = 0; i < 2 ; i++) {
@ -100,25 +100,25 @@ SA::HeunsMethod::~HeunsMethod() {
delete derivs[i];
}
}
void SA::HeunsMethod::step( double dt) {
(*derivs_func)( time, istate, derivs[0], udata);
void SA::HeunsMethod::step( double h) {
(*derivs_func)( indyVar, istate, derivs[0], udata);
for (unsigned int i=0; i<state_size; i++) {
wstate[i] = istate[i] + dt * derivs[0][i];
wstate[i] = istate[i] + h * derivs[0][i];
}
(*derivs_func)( time, wstate, derivs[1], udata);
(*derivs_func)( indyVar, wstate, derivs[1], udata);
for (unsigned int i=0; i<state_size; i++) {
ostate[i] = istate[i] + (dt/2) * ( derivs[0][i] + derivs[1][i] );
ostate[i] = istate[i] + (h/2) * ( derivs[0][i] + derivs[1][i] );
}
SA::FirstOrderODEVariableStepIntegrator::step(dt);
SA::FirstOrderODEVariableStepIntegrator::step(h);
}
void SA::HeunsMethod::step() {
step(default_dt);
step(default_h);
}
// ------------------------------------------------------------
SA::RK2Integrator::RK2Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEVariableStepIntegrator(dt, N, in_vars, out_vars, func, user_data)
SA::RK2Integrator::RK2Integrator( double h, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEVariableStepIntegrator(h, N, in_vars, out_vars, func, user_data)
{
wstate = new double[N];
for (unsigned int i = 0; i < 2 ; i++) {
@ -131,25 +131,25 @@ SA::RK2Integrator::~RK2Integrator() {
delete derivs[i];
}
}
void SA::RK2Integrator::step( double dt) {
(*derivs_func)( time, istate, derivs[0], udata);
void SA::RK2Integrator::step( double h) {
(*derivs_func)( indyVar, istate, derivs[0], udata);
for (unsigned int i=0; i<state_size; i++) {
wstate[i] = istate[i] + 0.5 * dt * derivs[0][i];
wstate[i] = istate[i] + 0.5 * h * derivs[0][i];
}
(*derivs_func)( time + 0.5 * dt , wstate, derivs[1], udata);
(*derivs_func)( indyVar + 0.5 * h , wstate, derivs[1], udata);
for (unsigned int i=0; i<state_size; i++) {
ostate[i] = istate[i] + dt * derivs[1][i];
ostate[i] = istate[i] + h * derivs[1][i];
}
SA::FirstOrderODEVariableStepIntegrator::step(dt);
SA::FirstOrderODEVariableStepIntegrator::step(h);
}
void SA::RK2Integrator::step() {
step(default_dt);
step(default_h);
}
// ------------------------------------------------------------
SA::RK4Integrator::RK4Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEVariableStepIntegrator(dt, N, in_vars, out_vars, func, user_data)
SA::RK4Integrator::RK4Integrator( double h, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEVariableStepIntegrator(h, N, in_vars, out_vars, func, user_data)
{
for (unsigned int i = 0; i < 3 ; i++) {
wstate[i] = new double[N];
@ -166,36 +166,36 @@ SA::RK4Integrator::~RK4Integrator() {
delete (derivs[i]);
}
}
void SA::RK4Integrator::step( double dt) {
(*derivs_func)( time, istate, derivs[0], udata);
void SA::RK4Integrator::step( double h) {
(*derivs_func)( indyVar, istate, derivs[0], udata);
for (unsigned int i=0; i<state_size; i++) {
wstate[0][i] = istate[i] + 0.5 * derivs[0][i] * dt;
wstate[0][i] = istate[i] + 0.5 * derivs[0][i] * h;
}
(*derivs_func)( time + 0.5 * dt , wstate[0], derivs[1], udata);
(*derivs_func)( indyVar + 0.5 * h , wstate[0], derivs[1], udata);
for (unsigned int i=0; i<state_size; i++) {
wstate[1][i] = istate[i] + 0.5 * derivs[1][i] * dt;
wstate[1][i] = istate[i] + 0.5 * derivs[1][i] * h;
}
(*derivs_func)( time + 0.5 * dt , wstate[1], derivs[2], udata);
(*derivs_func)( indyVar + 0.5 * h , wstate[1], derivs[2], udata);
for (unsigned int i=0; i<state_size; i++) {
wstate[2][i] = istate[i] + derivs[2][i] * dt;
wstate[2][i] = istate[i] + derivs[2][i] * h;
}
(*derivs_func)( time + dt , wstate[2], derivs[3], udata);
(*derivs_func)( indyVar + h , wstate[2], derivs[3], udata);
for (unsigned int i=0; i<state_size; i++) {
ostate[i] = istate[i] + ((1/6.0)* derivs[0][i] +
(1/3.0)* derivs[1][i] +
(1/3.0)* derivs[2][i] +
(1/6.0)* derivs[3][i]) * dt;
(1/6.0)* derivs[3][i]) * h;
}
SA::FirstOrderODEVariableStepIntegrator::step(dt);
SA::FirstOrderODEVariableStepIntegrator::step(h);
}
void SA::RK4Integrator::step() {
step(default_dt);
step(default_h);
}
// ------------------------------------------------------------
SA::RK3_8Integrator::RK3_8Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEVariableStepIntegrator(dt, N, in_vars, out_vars ,func, user_data)
SA::RK3_8Integrator::RK3_8Integrator( double h, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEVariableStepIntegrator(h, N, in_vars, out_vars ,func, user_data)
{
for (unsigned int i = 0; i < 3 ; i++) {
wstate[i] = new double[N];
@ -212,30 +212,30 @@ SA::RK3_8Integrator::~RK3_8Integrator() {
delete (derivs[i]);
}
}
void SA::RK3_8Integrator::step( double dt) {
(*derivs_func)( time, istate, derivs[0], udata);
void SA::RK3_8Integrator::step( double h) {
(*derivs_func)( indyVar, istate, derivs[0], udata);
for (unsigned int i=0; i<state_size; i++) {
wstate[0][i] = istate[i] + dt * (1/3.0) * derivs[0][i];
wstate[0][i] = istate[i] + h * (1/3.0) * derivs[0][i];
}
(*derivs_func)( time + (1/3.0) * dt , wstate[0], derivs[1], udata);
(*derivs_func)( indyVar + (1/3.0) * h , wstate[0], derivs[1], udata);
for (unsigned int i=0; i<state_size; i++) {
wstate[1][i] = istate[i] + dt * ((-1/3.0) * derivs[0][i] + derivs[1][i]);
wstate[1][i] = istate[i] + h * ((-1/3.0) * derivs[0][i] + derivs[1][i]);
}
(*derivs_func)( time + (2/3.0) * dt , wstate[1], derivs[2], udata);
(*derivs_func)( indyVar + (2/3.0) * h , wstate[1], derivs[2], udata);
for (unsigned int i=0; i<state_size; i++) {
wstate[2][i] = istate[i] + dt * (derivs[0][i] - derivs[1][i] + derivs[2][i]);
wstate[2][i] = istate[i] + h * (derivs[0][i] - derivs[1][i] + derivs[2][i]);
}
(*derivs_func)( time + dt , wstate[2], derivs[3], udata);
(*derivs_func)( indyVar + h , wstate[2], derivs[3], udata);
for (unsigned int i=0; i<state_size; i++) {
ostate[i] = istate[i] + ((1/8.0)* derivs[0][i] +
(3/8.0)* derivs[1][i] +
(3/8.0)* derivs[2][i] +
(1/8.0)* derivs[3][i]) * dt;
(1/8.0)* derivs[3][i]) * h;
}
SA::FirstOrderODEVariableStepIntegrator::step(dt);
SA::FirstOrderODEVariableStepIntegrator::step(h);
}
void SA::RK3_8Integrator::step() {
step(default_dt);
step(default_h);
}
// ------------------------------------------------------------
@ -269,18 +269,18 @@ SA::EulerCromerIntegrator::~EulerCromerIntegrator() {
delete v_out;
}
void SA::EulerCromerIntegrator::step( double dt) {
(*gderivs)( time, x_in, g_out, udata);
(*gderivs)( indyVar, x_in, g_out, udata);
for (int i=0 ; i<state_size; i++ ) {
v_out[i] = v_in[i] + g_out[i] * dt;
}
(*fderivs)( time, v_in, f_out, udata);
(*fderivs)( indyVar, v_in, f_out, udata);
for (int i=0 ; i<state_size; i++ ) {
x_out[i] = x_in[i] + f_out[i] * dt;
}
time += dt;
indyVar += dt;
}
void SA::EulerCromerIntegrator::step() {
step(default_dt);
step(default_h);
}
void SA::EulerCromerIntegrator::load() {
for (int i=0 ; i<state_size; i++ ) {
@ -312,8 +312,8 @@ void SA::EulerCromerIntegrator::unload(){
// {(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 dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
// : FirstOrderODEIntegrator(dt, N, in_vars, out_vars ,func, user_data) {
// SA::ABMIntegrator::ABMIntegrator ( int history_size, double h, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
// : FirstOrderODEIntegrator(h, N, in_vars, out_vars ,func, user_data) {
//
// 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.
@ -335,7 +335,7 @@ void SA::EulerCromerIntegrator::unload(){
// void SA::ABMIntegrator::step() {
//
// hix = (hix+1)%bufferSize; // Move history index forward
// (*derivs_func)( time, istate, deriv_history[hix], udata); // Calculated and store the deriv for current, corrected state.
// (*derivs_func)( indyVar, istate, deriv_history[hix], udata); // 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.
@ -344,19 +344,19 @@ void SA::EulerCromerIntegrator::unload(){
// 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];
// }
// ostate[i] = istate[i] + default_dt * composite_deriv[i];
// ostate[i] = istate[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)( time, ostate, deriv_history[hix], udata); // Calc deriv for predicted next state.
// (*derivs_func)( indyVar, ostate, deriv_history[hix], udata); // 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];
// }
// ostate[i] = istate[i] + default_dt * composite_deriv[i];
// ostate[i] = istate[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;
@ -372,8 +372,8 @@ void SA::EulerCromerIntegrator::unload(){
static const double AB2Coeffs[2] = {(3/2.0), (-1/2.0)};
static const double AM2Coeffs[2] = {(1/2.0), (1/2.0)};
SA::ABM2Integrator::ABM2Integrator ( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEIntegrator(dt, N, in_vars, out_vars ,func, user_data) {
SA::ABM2Integrator::ABM2Integrator ( double h, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEIntegrator(h, N, in_vars, out_vars ,func, user_data) {
// The amount of history that we intend to use in this ABM integrator.
bufferSize=2+1; // The size of the buffer needs to be one more than the history so that an integration step can be reset.
hix=0;
@ -393,7 +393,7 @@ SA::ABM2Integrator::ABM2Integrator ( double dt, int N, double* in_vars[], double
for (unsigned int i=0; i<N ; i++) {
priming_integrator_out_vars[i] = &(ostate[i]);
}
priming_integrator = new SA::RK2Integrator(dt, N, priming_integrator_in_vars, priming_integrator_out_vars, func, user_data);
priming_integrator = new SA::RK2Integrator(h, N, priming_integrator_in_vars, priming_integrator_out_vars, func, user_data);
}
SA::ABM2Integrator::~ABM2Integrator() {
for (int i=0; i<bufferSize ; i++) {
@ -406,7 +406,7 @@ SA::ABM2Integrator::~ABM2Integrator() {
void SA::ABM2Integrator::step() {
hix = (hix+1)%bufferSize; // Move history index forward
(*derivs_func)( time, istate, deriv_history[hix], udata); // Calculated and store the deriv for current, corrected state.
(*derivs_func)( indyVar, istate, deriv_history[hix], udata); // 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 < 2) ? currentHistorySize+1 : 2;
@ -421,19 +421,19 @@ void SA::ABM2Integrator::step() {
for (int n=0,j=hix; n<2 ; n++,j=(j+bufferSize-1)%bufferSize) {
composite_deriv[i] += AB2Coeffs[n] * deriv_history[j][i];
}
ostate[i] = istate[i] + default_dt * composite_deriv[i];
ostate[i] = istate[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)( time, ostate, deriv_history[hix], udata); // Calc deriv for predicted next state.
(*derivs_func)( indyVar, ostate, deriv_history[hix], udata); // 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<2 ; n++,j=(j+bufferSize-1)%bufferSize) {
composite_deriv[i] += AM2Coeffs[n] * deriv_history[j][i];
}
ostate[i] = istate[i] + default_dt * composite_deriv[i];
ostate[i] = istate[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;
@ -449,8 +449,8 @@ void SA::ABM2Integrator::undo_step() {
static const double AB4Coeffs[4] = {(55/24.0), (-59/24.0), (37/24.0), (-9/24.0)};
static const double AM4Coeffs[4] = {(9/24.0), (19/24.0), (-5/24.0), (1/24.0)};
SA::ABM4Integrator::ABM4Integrator ( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEIntegrator(dt, N, in_vars, out_vars ,func, user_data) {
SA::ABM4Integrator::ABM4Integrator ( double h, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data)
: FirstOrderODEIntegrator(h, N, in_vars, out_vars ,func, user_data) {
// The amount of history that we intend to use in this ABM integrator.
bufferSize=4+1; // The size of the buffer needs to be one more than the history so that an integration step can be reset.
@ -471,7 +471,7 @@ SA::ABM4Integrator::ABM4Integrator ( double dt, int N, double* in_vars[], double
for (unsigned int i=0; i<N ; i++) {
priming_integrator_out_vars[i] = &(ostate[i]);
}
priming_integrator = new SA::RK4Integrator(dt, N, priming_integrator_in_vars, priming_integrator_out_vars, func, user_data);
priming_integrator = new SA::RK4Integrator(h, N, priming_integrator_in_vars, priming_integrator_out_vars, func, user_data);
}
SA::ABM4Integrator::~ABM4Integrator() {
@ -486,7 +486,7 @@ SA::ABM4Integrator::~ABM4Integrator() {
void SA::ABM4Integrator::step() {
hix = (hix+1)%bufferSize; // Move history index forward
(*derivs_func)( time, istate, deriv_history[hix], udata); // Calculated and store the deriv for current, corrected state.
(*derivs_func)( indyVar, istate, deriv_history[hix], udata); // 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 < 4) ? currentHistorySize+1 : 4;
if ( currentHistorySize < 4 ) {
@ -500,19 +500,19 @@ void SA::ABM4Integrator::step() {
for (int n=0,j=hix; n<4 ; n++,j=(j+bufferSize-1)%bufferSize) {
composite_deriv[i] += AB4Coeffs[n] * deriv_history[j][i];
}
ostate[i] = istate[i] + default_dt * composite_deriv[i];
ostate[i] = istate[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)( time, ostate, deriv_history[hix], udata); // Calc deriv for predicted next state.
(*derivs_func)( indyVar, ostate, deriv_history[hix], udata); // 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<4 ; n++,j=(j+bufferSize-1)%bufferSize) {
composite_deriv[i] += AM4Coeffs[n] * deriv_history[j][i];
}
ostate[i] = istate[i] + default_dt * composite_deriv[i];
ostate[i] = istate[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;

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

@ -3,17 +3,22 @@
#include "SAIntegrator.hh"
#include <math.h>
void deriv4( double t, double state[], double derivs[], void* udata) {
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[], double derivs[], void* udata) {
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;
@ -21,8 +26,10 @@ void deriv2( double t, double state[], double derivs[], void* udata) {
derivs[3] = 5.0*t + 4.0;
}
void deriv1( double t, double state[], double derivs[], void* udata) {
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;