diff --git a/trick_source/trick_utils/SAIntegrator/examples/DoubleIntegral/DoubleIntegral.cpp b/trick_source/trick_utils/SAIntegrator/examples/DoubleIntegral/DoubleIntegral.cpp index 33e5480b..91822b07 100644 --- a/trick_source/trick_utils/SAIntegrator/examples/DoubleIntegral/DoubleIntegral.cpp +++ b/trick_source/trick_utils/SAIntegrator/examples/DoubleIntegral/DoubleIntegral.cpp @@ -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(); diff --git a/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh b/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh index f58b170c..41f41614 100644 --- a/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh +++ b/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh @@ -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(); diff --git a/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp b/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp index 7e1abf85..e9ae2659 100644 --- a/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp +++ b/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp @@ -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; diff --git a/trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc b/trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc index 59ef1175..3234385b 100644 --- a/trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc +++ b/trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc @@ -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;