diff --git a/trick_source/trick_utils/SAIntegrator/README.md b/trick_source/trick_utils/SAIntegrator/README.md new file mode 100644 index 00000000..b7036247 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/README.md @@ -0,0 +1,255 @@ +# Stand-Alone Integration Library + +## Contents + +* [Introduction](#Introduction) +* [class Integrator](#class-Integrator) +* [typedef derivsFunc](#typedef-derivsFunc) +* [class FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator) +* [class FirstOrderODEVariableStepIntegrator](#class-FirstOrderODEVariableStepIntegrator) +* [class EulerIntegrator](#class-EulerIntegrator) +* [class HeunsMethod](#class-HeunsMethod) +* [class RK2Integrator](#class-RK2Integrator) +* [class RK4Integrator](#class-RK4Integrator) +* [class RK3_8Integrator](#class-RK3_8Integrator) +* [class EulerCromerIntegrator](#class-EulerCromerIntegrator) +* [class ABM2Integrator](#class-ABM2Integrator) +* [class ABM4Integrator](#class-ABM4Integrator) + + +## Introduction +The Stand-Alone Integration Library can be used within a Trick simulation, or independent of it. + +Some examples of using these integrators can be found in the **examples/** directory. + +* [CannonBall](examples/CannonBall/README.md) uses the RK2Integrator. +* [MassSpringDamper](examples/MassSpringDamper/README.md) uses the EulerCromerIntegrator. +* [Orbit](examples/Orbit/README.md) uses the EulerCromerIntegrator. + + +## class Integrator + +### Description +This class represents an **integrator**. + +### Constructor + +#### ```Integrator(double dt, void* user_data);``` + + +|Parameter|Type |Description| +|---------|------------|--------------------| +|dt |```double```|Default time step value| +|user_data|```void*```| A pointer to user defined data that will be passed to a derivsFunc when called by the Integrator. | +### Member Functions + +#### ```virtual void step();``` + +Increment time by the default time-step. + +#### ```virtual void load() = 0;``` + +Load the input state. + +#### ```virtual void unload() = 0;``` + +Load the output state. + +#### ```double getTime();``` +Return the integrator time. + + +## typedef derivsFunc + +``` +typedef void (*derivsFunc)( double t, double state[], double derivs[], void* user_data); +``` +where: + +|Parameter|Type|Description| +|---|---|---| +|t|```double```| time| +|state|```double*```| (input) an array of state variable values | +|derivs|```double*```| (output) an array into which derivatives are to be returned| +|user_data|```void*```| a pointer to user_data | + +#### Example +``` +void my_derivs( double t, double state[], double deriv[], void* udata) { ... } +``` + + +## class FirstOrderODEIntegrator +Derived from [Integrator](#class-Integrator). + +### Description +This class represents an integrator for a first order [Ordinary Differential Equation (ODE)]([https://en.wikipedia.org/wiki/Ordinary_differential_equation). + +### Constructor + +``` +FirstOrderODEIntegrator( double dt, + int N, + double* in_vars[], + double* out_vars[], + derivsFunc func, + void* user_data); +``` +where: + + + +|Parameter|Type|Description| +|---|---|---| +|dt|```double```|Default time step value| +|N|```int```|Number of state variables to be integrated| +|in_vars|```double*```|Array of pointers to the state variables from which we ```load()``` the integrator state (```in_vars``` and ```out_vars``` will generally point to the same array of pointers.)| +|out_vars|```double*```|Array of pointers to the state variables to which we ```unload()``` the integrator state (```in_vars``` and ```out_vars``` will generally point to the same array of pointers.)| +| derivs_func |[derivsFunc](#typedef-derivsFunc)| A function that returns | +|user_data|```void*```| A pointer to user defined data that will be passed to a derivsFunc when called by the Integrator. | + +### Member Functions + +#### void step( double dt ) +Integrate over time step specified by **dt**. + +|Parameter|Type|Description| +|---|---|---| +|dt | double| A variable time-step | + +#### virtual void undo_step() + +If the integrator has not already been reset then : + +1. Decrement time by the last time-step dt. +2. Copy the input state back to the origin variables. +3. Set the integrator's 'reset' mode to true. + +Calling load() sets the 'reset' mode to false. + +#### void step() +Integrate over the default time-step. + +#### void load() +Load the integrator's initial state from the variables specified by **in_vars**. + +#### void unload(); + +Unload the integrator's result state to the variables specified by **out_vars**. + + +## class FirstOrderODEVariableStepIntegrator + +### Description +This class represents an integrator for a first order ODE, that be + +### Constructor +### Member Functions +#### virtual void step( double dt ); + + +## class EulerIntegrator +Derived from [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). + +### Description +The Euler method is a first order numerical integration method. It is the simplest, explicit Runge–Kutta method. +### Constructor +``` +EulerIntegrator(double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data) +``` +[Constructor Parameters](#FOODEConstructorParameters) are those [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). + + +## class HeunsMethod +Derived from [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). +### Description +This integrator implements +[Heun's Method](https://en.wikipedia.org/wiki/Heun%27s_method). + +### Constructor +``` +HeunsMethod( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data) +``` +[Constructor Parameters](#FOODEConstructorParameters) are those [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). + + +## class RK2Integrator +Derived from [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). +### Description +The Runga-Kutta-2 method is a second order, explicit, numerical integration method. +### Constructor +``` +RK2Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data) +``` +[Constructor Parameters](#FOODEConstructorParameters) are those [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). + + +## class RK4Integrator +Derived from [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). +### Description +The Runga-Kutta-4 method is a fourth order, explicit, numerical integration method. +### Constructor +``` +RK4Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data) +``` +[Constructor Parameters](#FOODEConstructorParameters) are those [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). + + +## class RK3_8Integrator +Derived from [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). +### Description +The Runga-Kutta-3/8 method is a fourth order, explicit, numerical integration method. +### Constructor +``` +RK3_8Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data) +``` +[Constructor Parameters](#FOODEConstructorParameters) are those [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). + + +## class ABM2Integrator +Derived from [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). +### Description + +The ABM2Integrator implements the second-order Adams-Bashforth-Moulton predictor/corrector method. Adams methods maintain a history of derivatives rather than calculating intermediate values like Runga-Kutta methods. + +### Constructor +``` +ABM2Integrator ( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data) +``` +[Constructor Parameters](#FOODEConstructorParameters) are those [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). + + +## class ABM4Integrator +Derived from [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). + +### Description +The ABM2Integrator implements the second-order Adams-Bashforth-Moulton predictor/corrector method. Adams methods maintain a history of derivatives rather than calculating intermediate values like Runga-Kutta methods. + +### Constructor +``` +ABM4Integrator ( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data) +``` + +[Constructor Parameters](#FOODEConstructorParameters) are those [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator). + + +## EulerCromerIntegrator +Derived from [class-Integrator](#class-Integrator). + +### Description +EulerCromer is integration method that conserves energy in oscillatory systems better than Runge-Kutta. So, it's good for mass-spring-damper systems, and orbital systems. + +### Constructor +``` +SemiImplicitEuler(double dt, int N, double* xp[], double* vp[], derivsFunc gfunc, derivsFunc ffunc, void* user_data) +``` + +|Parameter|Type |Description| +|-----------|-------------|-----------------------| +| dt |```double``` |Default time step value| +| N |```int``` |Number of state variables to be integrated| +| xp |```double*```|Array of pointers to the variables from which we ```load()``` and to which we ```unload()``` the integrator's position values .| +| vp |```double*```|Array of pointers to the variables from which we ```load()``` and to which we ```unload()``` the integrator's velocity values .| +| gfunc |[derivsFunc](#typedef-derivsFunc)| A function that returns acceleration | +| ffunc |[derivsFunc](#typedef-derivsFunc)| A function that returns velocity | +|user_data |```void*```| A pointer to user defined data that will be passed to a derivsFunc when called by the Integrator. | \ No newline at end of file diff --git a/trick_source/trick_utils/SAIntegrator/examples/CannonBall/CannonBall.cpp b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/CannonBall.cpp new file mode 100644 index 00000000..36512766 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/CannonBall.cpp @@ -0,0 +1,45 @@ +#include +#include +#include "SAIntegrator.hh" + +struct Cannon { + double pos[2]; + double vel[2]; +}; +void init_state ( Cannon& cannon ) { + cannon.pos[0] = 0.0; + cannon.pos[1] = 0.0; + cannon.vel[0] = 50.0 * cos(M_PI/6.0); + cannon.vel[1] = 50.0 * sin(M_PI/6.0); +} +void deriv( double t, double state[], double derivs[], void* udata) { + derivs[0] = state[2]; // vel_x + derivs[1] = state[3]; // vel_y + derivs[2] = 0.0; // acc_x + derivs[3] = -9.81; // acc_y +} +void print_header() { + printf ("t, cannon.pos[0], cannon.pos[1], cannon.vel[0], cannon.vel[1]\n"); +} +void print_state( double t, Cannon& cannon) { + printf ("%5.3f, %10.10f, %10.10f, %10.10f, %10.10f\n", t, cannon.pos[0], cannon.pos[1], cannon.vel[0], cannon.vel[1]); +} +int main ( int argc, char* argv[]) { + Cannon cannon; + double* state_var_p[4] = { &(cannon.pos[0]), &(cannon.pos[1]), &(cannon.vel[0]), &(cannon.vel[1]) }; + unsigned count = 0; + double dt = 0.01; + double t = count * dt; + SA::RK2Integrator integ(dt, 4, state_var_p, state_var_p, deriv, NULL); + init_state(cannon); + print_header(); + print_state(t, cannon ); + while (t < 5.1) { + integ.load(); + integ.step(); + integ.unload(); + count ++; + t = count * dt; + print_state(t, cannon ); + } +} diff --git a/trick_source/trick_utils/SAIntegrator/examples/CannonBall/README.md b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/README.md new file mode 100644 index 00000000..dfb30fcf --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/README.md @@ -0,0 +1,11 @@ +# CannonBall + +This is an example of using the RK2Integrator to create a simple cannon ball simulation. + +''' +$ make +$ ./CannonBall > cannon.csv +$ python plot_trajectory.py +''' + +![Cannon](images/Cannon.png) diff --git a/trick_source/trick_utils/SAIntegrator/examples/CannonBall/images/Cannon.png b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/images/Cannon.png new file mode 100644 index 00000000..866afeb2 Binary files /dev/null and b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/images/Cannon.png differ diff --git a/trick_source/trick_utils/SAIntegrator/examples/CannonBall/makefile b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/makefile new file mode 100644 index 00000000..3bf6eb8d --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/makefile @@ -0,0 +1,20 @@ + +RM = rm -rf +CC = cc +CPP = c++ + +CXXFLAGS = -g -Wall +INCLUDE_DIRS = -I../../include +LIBDIR = ../../lib + +all: CannonBall + +CannonBall: + $(CPP) $(CXXFLAGS) CannonBall.cpp ${INCLUDE_DIRS} -L${LIBDIR} -lSAInteg -o CannonBall + +clean: + ${RM} CannonBall.dSYM + +spotless: clean + ${RM} CannonBall + ${RM} cannon.csv diff --git a/trick_source/trick_utils/SAIntegrator/examples/CannonBall/plot_trajectory.py b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/plot_trajectory.py new file mode 100755 index 00000000..942d1fe8 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/plot_trajectory.py @@ -0,0 +1,18 @@ +#!/usr/bin/env python + +import matplotlib.pyplot as plt +import numpy as np +data = np.genfromtxt('cannon.csv', + delimiter=',', + skip_header=1, + skip_footer=1, + names=['t', 'px', 'py', 'vx', 'vx'], + dtype=(float, float, float, float, float) + ) + +curve1 = plt.plot(data['px'], data['py'], 'C1-') +plt.title('Cannonball Trajectory') +plt.xlabel('pos-x') +plt.ylabel('pos-y') +plt.grid(True) +plt.show() diff --git a/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/MassSpringDamper.cpp b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/MassSpringDamper.cpp new file mode 100644 index 00000000..a1aa2fc8 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/MassSpringDamper.cpp @@ -0,0 +1,52 @@ +#include +#include +#include "SAIntegrator.hh" +#define N_DIMENSIONS 1 +struct MassSpringDamper { + double pos; + double vel; + double k; // Spring constant + double c; // Damping constant + double mass; +}; +void init_state ( MassSpringDamper& msd ) { + msd.pos = 2.0; + msd.vel = 0.0; + msd.k = 1.0; + msd.c = 0.0; + msd.mass = 5.0; +} +void print_header() { + printf ("time, msd.pos, msd.vel\n"); +} +void print_state( double t, MassSpringDamper& msd ) { + printf ("%10.10f, %10.10f, %10.10f\n", t, msd.pos, msd.vel); +} +void G( double t, double x[], double g_out[], void* udata) { + MassSpringDamper* msd = (MassSpringDamper*)udata; + g_out[0] = -(msd->k/msd->mass) * x[0]; +} +void F( double t, double v[], double f_out[], void* udata) { + f_out[0] = v[0]; + f_out[1] = v[1]; +} +int main ( int argc, char* argv[]) { + MassSpringDamper msd; + double* x_p[N_DIMENSIONS] = { &msd.pos }; + double* v_p[N_DIMENSIONS] = { &msd.vel }; + unsigned count = 0; + double dt = 0.001; + double time = count * dt; + init_state(msd); + print_header(); + print_state(time, msd); + SA::EulerCromerIntegrator integ(dt, N_DIMENSIONS, x_p, v_p, G, F, &msd); + while (time < 50) { + integ.load(); + integ.step(); + integ.unload(); + count ++; + time = count * dt; + print_state(time, msd); + } +} diff --git a/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/README.md b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/README.md new file mode 100644 index 00000000..da535eca --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/README.md @@ -0,0 +1,22 @@ +# MassSpringDamper + +This program uses the SemiImplicitEuler integrator to simulate a simple mass-spring-damper. + +Generate the results as follows: +''' +$ make +$ ./MassSpringDamper > msd.csv +''' + +Plot the results as follows: + +''' +$ python plot_position.py +''' + +and + +''' +$ python plot_velocity.py +''' +![MSD](images/MSD.png) diff --git a/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/images/MSD.png b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/images/MSD.png new file mode 100644 index 00000000..d0c894f8 Binary files /dev/null and b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/images/MSD.png differ diff --git a/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/makefile b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/makefile new file mode 100644 index 00000000..bda3abb5 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/makefile @@ -0,0 +1,20 @@ + +RM = rm -rf +CC = cc +CPP = c++ + +CXXFLAGS = -g -Wall +INCLUDE_DIRS = -I../../include +LIBDIR = ../../lib + +all: MassSpringDamper + +MassSpringDamper: MassSpringDamper.cpp + $(CPP) $(CXXFLAGS) MassSpringDamper.cpp ${INCLUDE_DIRS} -L${LIBDIR} -lSAInteg -o MassSpringDamper + +clean: + ${RM} MassSpringDamper.dSYM + +spotless: clean + ${RM} MassSpringDamper + ${RM} msd.csv diff --git a/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/plot_position.py b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/plot_position.py new file mode 100644 index 00000000..359e2cdc --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/plot_position.py @@ -0,0 +1,18 @@ +#!/usr/bin/env python + +import matplotlib.pyplot as plt +import numpy as np +data = np.genfromtxt('msd.csv', + delimiter=',', + skip_header=1, + skip_footer=1, + names=['t', 'pos', 'vel'], + dtype=(float, float, float) + ) + +curve1 = plt.plot(data['t'], data['pos'], 'C1-') +plt.title('MSD Position') +plt.xlabel('position') +plt.ylabel('time') +plt.grid(True) +plt.show() diff --git a/trick_source/trick_utils/SAIntegrator/examples/Orbit/Orbit.cpp b/trick_source/trick_utils/SAIntegrator/examples/Orbit/Orbit.cpp new file mode 100644 index 00000000..85c56f0b --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/Orbit/Orbit.cpp @@ -0,0 +1,57 @@ +#include +#include +#include "SAIntegrator.hh" + +#define GRAVITATIONAL_CONSTANT 6.674e-11 +#define EARTH_MASS 5.9723e24 +#define EARTH_RADIUS 6367500.0 + +struct Orbit { + double pos[2]; + double vel[2]; + double planet_mass; +}; +void init_state ( Orbit& orbit ) { + orbit.pos[0] = 0.0; + orbit.pos[1] = 6578000.0; + orbit.vel[0] = 7905.0; + orbit.vel[1] = 0.0; + orbit.planet_mass = EARTH_MASS; +} +void print_header() { + printf ("time, orbit.pos[0], orbit.pos[1], orbit.vel[0], orbit.vel[1]\n"); +} +void print_state( double t, Orbit& orbit ) { + printf ("%10.10f, %10.10f, %10.10f, %10.10f, %10.10f\n", + t, orbit.pos[0], orbit.pos[1], orbit.vel[0], orbit.vel[1]); +} +void G( double t, double x[], double g_out[], void* udata) { + Orbit* orbit = (Orbit*)udata; + double d = sqrt( x[0]*x[0] + x[1]*x[1]); + g_out[0] = -x[0] * GRAVITATIONAL_CONSTANT * orbit->planet_mass / (d*d*d); + g_out[1] = -x[1] * GRAVITATIONAL_CONSTANT * orbit->planet_mass / (d*d*d); +} +void F( double t, double v[], double f_out[], void* udata) { + f_out[0] = v[0]; + f_out[1] = v[1]; +} +int main ( int argc, char* argv[]) { + Orbit orbit; + double* x_p[2] = { &orbit.pos[0], &orbit.pos[1] }; + double* v_p[2] = { &orbit.vel[0], &orbit.vel[1] }; + unsigned count = 0; + double dt = 0.1; + double time = count * dt; + init_state(orbit); + print_header(); + print_state(time, orbit); + SA::EulerCromerIntegrator integ(dt, 2, x_p, v_p, G, F, &orbit); + while (time < 6000) { + integ.load(); + integ.step(); + integ.unload(); + count ++; + time = count * dt; + print_state(time, orbit); + } +} diff --git a/trick_source/trick_utils/SAIntegrator/examples/Orbit/README.md b/trick_source/trick_utils/SAIntegrator/examples/Orbit/README.md new file mode 100644 index 00000000..6f87c0a3 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/Orbit/README.md @@ -0,0 +1,17 @@ +# Orbit + +This program uses the *SemiImplicitEuler* integrator to simulate the orbit of an Earth satellite. + +Generate the results as follows: +''' +$ make +$ ./Orbit > orbit.csv +''' + +Plot the results as follows: + +''' +$ python plot_position.py +''' + +![Orbit](images/Orbit.png) diff --git a/trick_source/trick_utils/SAIntegrator/examples/Orbit/images/Orbit.png b/trick_source/trick_utils/SAIntegrator/examples/Orbit/images/Orbit.png new file mode 100644 index 00000000..45e73e98 Binary files /dev/null and b/trick_source/trick_utils/SAIntegrator/examples/Orbit/images/Orbit.png differ diff --git a/trick_source/trick_utils/SAIntegrator/examples/Orbit/makefile b/trick_source/trick_utils/SAIntegrator/examples/Orbit/makefile new file mode 100644 index 00000000..5d155c26 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/Orbit/makefile @@ -0,0 +1,20 @@ + +RM = rm -rf +CC = cc +CPP = c++ + +CXXFLAGS = -g -Wall +INCLUDE_DIRS = -I../../include +LIBDIR = ../../lib + +all: Orbit + +Orbit: Orbit.cpp + $(CPP) $(CXXFLAGS) Orbit.cpp ${INCLUDE_DIRS} -L${LIBDIR} -lSAInteg -o Orbit + +clean: + ${RM} Orbit.dSYM + +spotless: clean + ${RM} Orbit + ${RM} orbit.csv diff --git a/trick_source/trick_utils/SAIntegrator/examples/Orbit/plot_position.py b/trick_source/trick_utils/SAIntegrator/examples/Orbit/plot_position.py new file mode 100755 index 00000000..84beaee3 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/Orbit/plot_position.py @@ -0,0 +1,18 @@ +#!/usr/bin/env python + +import matplotlib.pyplot as plt +import numpy as np +data = np.genfromtxt('orbit.csv', + delimiter=',', + skip_header=1, + skip_footer=1, + names=['t', 'posx','posy','velx','vely'], + dtype=(float, float, float, float, float) + ) + +curve1 = plt.plot(data['posx'], data['posy'], 'C1-') +plt.title('Orbit') +plt.xlabel('position') +plt.ylabel('position') +plt.grid(True) +plt.show() diff --git a/trick_source/trick_utils/SAIntegrator/examples/makefile b/trick_source/trick_utils/SAIntegrator/examples/makefile new file mode 100644 index 00000000..01db8baa --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/makefile @@ -0,0 +1,15 @@ + +examples: + @make -C CannonBall + @make -C MassSpringDamper + @make -C Orbit + +clean: + @make -C CannonBall clean + @make -C MassSpringDamper clean + @make -C Orbit clean + +spotless: + @make -C CannonBall spotless + @make -C MassSpringDamper spotless + @make -C Orbit spotless diff --git a/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh b/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh new file mode 100644 index 00000000..f1081cdf --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh @@ -0,0 +1,177 @@ + +namespace SA { + + class Integrator { + protected: + double time; + double default_dt; + void * udata; + public: + Integrator(double dt, void* user_data): time(0.0), default_dt(dt), udata(user_data) {}; + virtual ~Integrator() {} + virtual void step() { time += default_dt; }; // Integrate over implicit (default) time step (default: default_dt) + virtual void load() = 0; + virtual void unload() = 0; + double getTime() {return time;} + }; + + typedef void (*derivsFunc)( double t, double state[], double derivs[], void* user_data); + + class FirstOrderODEIntegrator : public Integrator { + protected: + unsigned int state_size; // size of the state vector. + double* istate; + double* ostate; + double** orig_vars; + double** dest_vars; + derivsFunc derivs_func; + + bool reset; + double last_dt; // 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(); + // void init(double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func); + virtual void undo_step(); + void load(); + void unload(); + }; + + class FirstOrderODEVariableStepIntegrator : public FirstOrderODEIntegrator { + + public: + FirstOrderODEVariableStepIntegrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data); + using Integrator::step; + virtual void step( double dt ); + }; + + 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(); + void step( double dt); + void step(); + }; + + class HeunsMethod : public FirstOrderODEVariableStepIntegrator { + public: + double *wstate; + double *derivs[2]; + + HeunsMethod( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data); + ~HeunsMethod(); + void step( double dt); + void step(); + }; + + class RK2Integrator : public FirstOrderODEVariableStepIntegrator { + public: + double *wstate; + double *derivs[2]; + + RK2Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data); + ~RK2Integrator(); + void step( double dt); + void step(); + }; + + class RK4Integrator : public FirstOrderODEVariableStepIntegrator { + public: + double *wstate[3]; + double *derivs[4]; + + RK4Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data); + ~RK4Integrator(); + void step( double dt); + void step(); + + }; + + class RK3_8Integrator : public FirstOrderODEVariableStepIntegrator { + public: + 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(); + void step( double dt); + void 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; + double* g_out; + double* f_out; + derivsFunc gderivs; + derivsFunc fderivs; + + public: + EulerCromerIntegrator(double dt, int N, double* xp[], double* vp[], derivsFunc g, derivsFunc f , void* user_data); + ~EulerCromerIntegrator(); + void step( double dt); + void step(); + void load(); + void unload(); + }; + + // AdamsBashforthMoulton 4 + // 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 dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data); + // ~ABMIntegrator(); + // void step(); + // void undo_step(); + // }; + + // AdamsBashforthMoulton 2 + class ABM2Integrator : public FirstOrderODEIntegrator { + protected: + double ** deriv_history; + double * composite_deriv; + unsigned int bufferSize; + unsigned int hix; + unsigned int currentHistorySize; + FirstOrderODEVariableStepIntegrator* priming_integrator; + + public: + ABM2Integrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc derivs_func, void* user_data); + ~ABM2Integrator(); + void step(); + void undo_step(); + }; + + // AdamsBashforthMoulton 4 + class ABM4Integrator : public FirstOrderODEIntegrator { + protected: + double ** deriv_history; + double * composite_deriv; + unsigned int bufferSize; + unsigned int hix; + unsigned int currentHistorySize; + FirstOrderODEVariableStepIntegrator* priming_integrator; + + public: + ABM4Integrator( double dt, 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/makefile b/trick_source/trick_utils/SAIntegrator/makefile new file mode 100644 index 00000000..f299e378 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/makefile @@ -0,0 +1,44 @@ + +RM = rm -rf +CC = cc +CPP = c++ + +CFLAGS = -g -Wall +INCLUDE_DIRS = -Iinclude +OBJDIR = obj +LIBDIR = lib +LIBNAME = libSAInteg.a +LIBOBJS = ${OBJDIR}/SAIntegrator.o + +all: test examples + +test: ${LIBDIR}/${LIBNAME} + ${MAKE} -C test + +examples: ${LIBDIR}/${LIBNAME} + ${MAKE} -C examples + +$(LIBOBJS): $(OBJDIR)/%.o : src/%.cpp | $(OBJDIR) + $(CPP) $(CFLAGS) ${INCLUDE_DIRS} -c $< -o $@ + +${LIBDIR}/${LIBNAME}: ${LIBOBJS} | ${LIBDIR} + ar crs $@ ${LIBOBJS} + +${OBJDIR}: + mkdir -p ${OBJDIR} + +${LIBDIR}: + mkdir -p ${LIBDIR} + +clean: + ${RM} *~ + ${RM} ${OBJDIR} + ${MAKE} -C test clean + ${MAKE} -C examples clean + +spotless: + ${RM} *~ + ${RM} ${OBJDIR} + ${RM} ${LIBDIR} + ${MAKE} -C test spotless + ${MAKE} -C examples spotless diff --git a/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp b/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp new file mode 100644 index 00000000..a5756819 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp @@ -0,0 +1,511 @@ + +#include +#include +#include "SAIntegrator.hh" + +// ------------------------------------------------------------ + +SA::FirstOrderODEIntegrator::FirstOrderODEIntegrator(double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data) +: Integrator(dt, user_data) { + state_size = N; + orig_vars = in_vars; + dest_vars = out_vars; + istate = new double[state_size]; + ostate = new double[state_size]; + derivs_func = func; + reset=false; +} +SA::FirstOrderODEIntegrator::~FirstOrderODEIntegrator() { + delete istate; + delete ostate; +} +void SA::FirstOrderODEIntegrator::load() { + reset = false; + for (int i=0 ; iload(); + priming_integrator->step(); + priming_integrator->unload(); + } else { + // (Predictor) Predict the next state using the Adams-Bashforth explicit method. + for (int i=0; iload(); + priming_integrator->step(); + priming_integrator->unload(); + } else { + // (Predictor) Predict the next state using the Adams-Bashforth explicit method. + for (int i=0; i +#include +#include "SAIntegrator.hh" +#include + +void deriv4( double t, double state[], double derivs[], void* udata) { + 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) { + + 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, double state[], double derivs[], void* udata) { + + derivs[0] = 4.0; + derivs[1] = -4.0; + derivs[2] = M_PI; + derivs[3] = -M_PI; +} + +#define EXCEPTABLE_ERROR 0.00000000001 + +TEST(Integ_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(Integ_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(Integ_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(Integ_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); +} diff --git a/trick_source/trick_utils/SAIntegrator/test/makefile b/trick_source/trick_utils/SAIntegrator/test/makefile new file mode 100644 index 00000000..a2438e30 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/test/makefile @@ -0,0 +1,28 @@ + +RM = rm -rf +CC = cc +CPP = c++ + +CFLAGS = -g -Wall +INCLUDE_DIRS = -I../include +OBJDIR = obj +LIBDIR = ../lib +LIBNAME = libSAInteg.a +LIBOBJS = ${OBJDIR}/Integrator.o + +all: Integ_unittest + +Integ_unittest.o : Integ_unittest.cc + $(CPP) $(CFLAGS) $(INCLUDE_DIRS) -c $< + +Integ_unittest : ${LIBDIR}/${LIBNAME} Integ_unittest.o + $(CPP) $(CFLAGS) -o $@ $^ -lgtest -lgtest_main -lpthread + +${LIBDIR}/${LIBNAME} : + $(MAKE) -C .. + +clean: + ${RM} *.o + +spotless: clean + ${RM} Integ_unittest