Add Recursive Dynamic Event Processing. Ref #995

This commit is contained in:
Penn, John M 047828115 2020-04-29 17:26:02 -05:00
parent 119b60d7c5
commit 479f3458cd
22 changed files with 238 additions and 118 deletions

View File

@ -404,14 +404,17 @@ namespace Trick {
*/
virtual int integrate_dt (double beg_time, double del_time);
/**
* Process dynamic events.
*
* @return Zero/non-zero success indicator.
* Out-of-sync integrators cause a non-zero return.
* @param end_time Time at the end of the integration interval.
* @param start_t Time at the beginning of the integration interval.
* @param end_t Time at the end of the integration interval.
*/
int process_dynamic_events (double end_time);
int process_dynamic_events (double start_t, double end_t, unsigned int depth=0);
};
}

View File

@ -65,6 +65,7 @@ namespace Trick {
#ifndef SWIGPYTHON
void state_reset ();
void state_element_in (unsigned int index , double* state_p);
#endif
#ifndef SWIGPYTHON

View File

@ -13,6 +13,9 @@ int integrate_2nd_order_ode(const double* acc, double* vel, double * pos);
double get_integ_time(void);
void set_integ_time(double time_value);
void reset_state();
#ifndef USE_ER7_UTILS_INTEGRATORS
void load_state_element(unsigned int index , double* state);
#endif
void load_state(double* arg1, ... );
void load_deriv( double* arg1, ...);
void load_deriv2( double* arg1, ...);

View File

@ -6,8 +6,9 @@ trick.var_set_copy_mode(2)
dyn.contact.nballs = 2
dyn.contact.balls = trick.TMM_declare_var_1d("Ball*", dyn.contact.nballs)
# x, y, vx, vy, r, m
dyn.contact.balls[0] = trick.make_Ball(0.0, 0.0, 0.5, 0.0, 0.5, 1.0)
dyn.contact.balls[1] = trick.make_Ball(3.0, 0.0, -0.5, 0.0, 0.5, 1.0)
dyn.contact.balls[1] = trick.make_Ball(3.1415926535, 0.0, -0.5, 0.0, 0.5, 1.0)
dyn_integloop.getIntegrator(trick.Euler, 2*dyn.contact.nballs)

View File

@ -6,13 +6,13 @@ trick.var_set_copy_mode(2)
dyn.contact.nballs = 7
dyn.contact.balls = trick.TMM_declare_var_1d("Ball*", dyn.contact.nballs)
# x, y, vx, vy, r
dyn.contact.balls[0] = trick.make_Ball(-5.0, 0.1, 0.5, 0.0, 0.5, 1.0)
dyn.contact.balls[0] = trick.make_Ball(-5.0, 0.1, 5.0, 0.0, 0.5, 1.0)
dyn.contact.balls[1] = trick.make_Ball( 0.0, 0.0, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[2] = trick.make_Ball( 1.2, 0.6, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[3] = trick.make_Ball( 1.2,-0.6, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[4] = trick.make_Ball( 2.4, 1.2, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[5] = trick.make_Ball( 2.4, 0.0, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[6] = trick.make_Ball( 2.4,-1.2, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[2] = trick.make_Ball( 1.05, 0.505, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[3] = trick.make_Ball( 1.05,-0.505, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[4] = trick.make_Ball( 2.10, 1.010, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[5] = trick.make_Ball( 2.10, 0.0, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[6] = trick.make_Ball( 2.10,-1.005, 0.0, 0.0, 0.5, 1.0)
dyn_integloop.getIntegrator(trick.Euler, 2*dyn.contact.nballs)

View File

@ -0,0 +1,34 @@
execfile("Modified_data/realtime.py")
# Variable Server Data should be copied at top of frame.
trick.var_set_copy_mode(2)
dyn.contact.nballs = 7
dyn.contact.balls = trick.TMM_declare_var_1d("Ball*", dyn.contact.nballs)
# x, y, vx, vy, r, m
dyn.contact.balls[0] = trick.make_Ball(-4.00, 0.0, 2.0, 0.0, 0.5, 1.0)
dyn.contact.balls[1] = trick.make_Ball(-1.00, 0.0, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[2] = trick.make_Ball( 0.01, 0.0, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[3] = trick.make_Ball( 1.02, 0.0, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[4] = trick.make_Ball( 2.03, 0.0, 0.0, 0.0, 0.5, 1.0)
dyn.contact.balls[5] = trick.make_Ball( 7.00, 0.0, 0.0, 0.0, 1.0, 1000000.0)
dyn.contact.balls[6] = trick.make_Ball(-7.00, 0.0, 0.0, 0.0, 1.0, 1000000.0)
dyn_integloop.getIntegrator(trick.Euler, 2*dyn.contact.nballs)
#==========================================
# Start the Satellite Graphics Client
#==========================================
varServerPort = trick.var_server_get_port();
BallDisplay_path = "models/graphics/dist/BallDisplay.jar"
if (os.path.isfile(BallDisplay_path)) :
BallDisplay_cmd = "java -jar " \
+ BallDisplay_path \
+ " " + str(varServerPort) + " &" ;
print(BallDisplay_cmd)
os.system( BallDisplay_cmd);
else :
print('==================================================================================')
print('BallDisplay needs to be built. Please \"cd\" into ../models/graphics and type \"make\".')
print('==================================================================================')

View File

@ -14,9 +14,10 @@ class Contact {
Contact(){}
Ball ** balls;
unsigned int nballs;
REGULA_FALSI rf ;
unsigned int numAssociations ;
REGULA_FALSI* ballAssociations ;
void ballCollision(Ball &b1, Ball &b2);
int default_data();
int state_init();
int state_deriv();

View File

@ -7,6 +7,7 @@ LIBRARY DEPENDENCY:
#include <math.h>
#include <iostream>
#include "trick/integrator_c_intf.h"
#include "trick/memorymanager_c_intf.h"
#include "../include/Contact.hh"
int Contact::default_data() {
@ -16,8 +17,20 @@ int Contact::default_data() {
}
int Contact::state_init() {
rf.mode = Decreasing;
rf.error_tol = 0.001;
double now ;
numAssociations = (nballs*(nballs-1))/2;
ballAssociations = (REGULA_FALSI*)TMM_declare_var_1d("REGULA_FALSI", numAssociations);
unsigned int ii,jj;
for (ii=1; ii<nballs; ii++) {
for (jj=0; jj<ii; jj++) {
unsigned int association_index = ii*(ii-1)/2+jj;
ballAssociations[association_index].mode = Decreasing;
ballAssociations[association_index].error_tol = 0.0000001;
now = get_integ_time() ;
reset_regula_falsi( now, &ballAssociations[association_index] );
}
}
return (0);
}
@ -28,9 +41,15 @@ int Contact::state_deriv() {
int Contact::state_integ() {
int integration_step;
for (unsigned int ii = 0; ii < nballs; ii++) {
// Be sure to fix load_indexed_state() so it only loads @ intermediate_step == 0
#ifdef USE_ER7_UTILS_INTEGRATORS
load_indexed_state( 2*ii, balls[ii]->pos[0]);
load_indexed_state( 2*ii+1, balls[ii]->pos[1]);
#else
load_state_element( 2*ii, &balls[ii]->pos[0]);
load_state_element( 2*ii+1, &balls[ii]->pos[1]);
#endif
}
for (unsigned int ii = 0; ii < nballs; ii++) {
load_indexed_deriv( 2*ii, balls[ii]->vel[0]);
@ -70,32 +89,55 @@ void Contact::ballCollision(Ball &b1, Ball &b2) {
}
double Contact::collision() {
double tgo ; /* time-to-go */
//double tgo ; /* time-to-go */
double now ; /* current integration time. */
int first, second;
unsigned int first, second;
double min_dist = 1000000.0;
for (unsigned int ii = 0; ii < nballs; ii++) {
for (unsigned int jj = ii+1; jj < nballs; jj++) {
unsigned int association_index;
double event_tgo;
unsigned int ii,jj;
now = get_integ_time() ;
event_tgo = BIG_TGO;
for (ii=1; ii<nballs; ii++) {
for (jj=0; jj<ii; jj++) {
association_index = ii*(ii-1)/2+jj;
double xdiff = balls[ii]->pos[0] - balls[jj]->pos[0];
double ydiff = balls[ii]->pos[1] - balls[jj]->pos[1];
double dist = sqrt( xdiff * xdiff + ydiff * ydiff) - (balls[ii]->radius + balls[jj]->radius);
if (dist < min_dist) {
min_dist = dist;
double distance = sqrt( xdiff * xdiff + ydiff * ydiff) - (balls[ii]->radius + balls[jj]->radius);
ballAssociations[association_index].error = distance;
double tgo = regula_falsi( now, &ballAssociations[association_index] ) ;
if (tgo < BIG_TGO) {
first = ii;
second = jj;
event_tgo = tgo;
}
}
}
rf.error = min_dist;
now = get_integ_time() ; /* Get the current integration time */
tgo = regula_falsi( now, &rf ) ; /* Estimate remaining integration time. */
if (tgo == 0.0) { /* If we are at the event, it's action time! */
if (event_tgo == 0.0) {
now = get_integ_time() ;
ballCollision(*balls[first], *balls[second]);
std::cout << "Ball["<<first<<"] and Ball["<<second<<"] collide at t = "<<now<<"." << std::endl;
reset_regula_falsi( now, &rf );
std::cout << "Balls "<<second<<" and "<<first<<" collide at t = "<<now<<"." << std::endl;
//
for (ii=1; ii<nballs; ii++) {
for (jj=0; jj<ii; jj++) {
association_index = ii*(ii-1)/2+jj;
if ((ii==first) && (jj==second)) {
reset_regula_falsi( now, &ballAssociations[association_index] );
} else {
//Set X_upper and t_upper to event values.
double xdiff = balls[ii]->pos[0] - balls[jj]->pos[0];
double ydiff = balls[ii]->pos[1] - balls[jj]->pos[1];
double distance = sqrt( xdiff * xdiff + ydiff * ydiff) - (balls[ii]->radius + balls[jj]->radius);
ballAssociations[association_index].error = distance;
ballAssociations[association_index].x_upper = distance;
ballAssociations[association_index].t_upper = now;
ballAssociations[association_index].upper_set = 1;
}
return (tgo) ;
}
}
}
return (event_tgo) ;
}

View File

@ -10,11 +10,13 @@
#include "trick/exec_proto.h"
#include "trick/message_proto.h"
#include "trick/message_type.h"
#include "trick/JobData.hh"
// System includes
#include <iostream>
#include <iomanip>
#include <cstdarg>
#include <math.h>
// Anonymous namespace for local functions
@ -329,22 +331,24 @@ void Trick::IntegLoopScheduler::call_deriv_jobs ()
*/
int Trick::IntegLoopScheduler::integrate()
{
double end_time = exec_get_sim_time();
double beg_time = end_time - next_cycle;
double t_end = exec_get_sim_time();
double t_start = t_end - next_cycle; // This is the time of the current state vector.
int status;
// Call all of the jobs in the pre-integration job queue.
call_jobs (pre_integ_jobs);
// Integrate sim objects to the current time.
status = integrate_dt (beg_time, next_cycle);
status = integrate_dt (t_start, next_cycle);
if (status != 0) {
return status;
}
// Process dynamic events, if any.
if (! dynamic_event_jobs.empty()) {
status = process_dynamic_events (end_time);
status = process_dynamic_events (t_start, t_end);
if (status != 0) {
return status;
}
@ -434,10 +438,8 @@ bool Trick::IntegLoopScheduler::get_last_step_deriv()
/**
Integrate over the specified time interval.
*/
int Trick::IntegLoopScheduler::integrate_dt (
double beg_time,
double dt)
{
int Trick::IntegLoopScheduler::integrate_dt ( double t_start, double dt) {
int ipass = 0;
int ex_pass = 0;
bool need_derivs = get_first_step_deriv_from_integrator();
@ -455,8 +457,8 @@ int Trick::IntegLoopScheduler::integrate_dt (
integ_jobs.reset_curr_index();
while ((curr_job = integ_jobs.get_next_job()) != NULL) {
void* sup_class_data = curr_job->sup_class_data;
void* sup_class_data = curr_job->sup_class_data;
// Jobs without supplemental data use the default integrator.
if (sup_class_data == NULL) {
trick_curr_integ = integ_ptr;
@ -468,6 +470,7 @@ int Trick::IntegLoopScheduler::integrate_dt (
*(static_cast<Trick::Integrator**>(sup_class_data));
}
if (trick_curr_integ == NULL) {
message_publish (
MSG_ERROR,
@ -477,13 +480,13 @@ int Trick::IntegLoopScheduler::integrate_dt (
}
if (ex_pass == 1) {
trick_curr_integ->time = beg_time;
trick_curr_integ->time = t_start;
trick_curr_integ->dt = dt;
}
if (verbosity || trick_curr_integ->verbosity) {
message_publish (MSG_DEBUG, "Job: %s, time: %f, dt: %f\n",
curr_job->name.c_str(), beg_time, dt);
curr_job->name.c_str(), t_start, dt);
}
ipass = curr_job->call();
@ -505,51 +508,55 @@ int Trick::IntegLoopScheduler::integrate_dt (
return 0;
}
/**
Process the dynamic event queue.
@param end_time The time the curr_time variable is set and be used for integration time.
*/
int Trick::IntegLoopScheduler::process_dynamic_events (double end_time)
{
int Trick::IntegLoopScheduler::process_dynamic_events ( double t_start, double t_end, unsigned int depth) {
bool fired = false;
double end_offset = 1e-15 * nominal_cycle;
double curr_time = end_time;
double t_to = t_end;
double t_from = t_start;
Trick::JobData * curr_job;
dynamic_event_jobs.reset_curr_index();
while ((curr_job = dynamic_event_jobs.get_next_job()) != NULL) {
// Call the dynamic event job to calculate an estimated time-to-go.
double tgo = curr_job->call_double();
// If there is a root in this interval ...
if (tgo < end_offset) {
if (tgo < 0.0) { // If there is a root in this interval ...
fired = true;
// Search for the event.
// Dynamic events return 0.0, exactly, to indicate that the
// event has been found and has been triggered.
while (tgo != 0.0) {
// Integrate to the estimated event time.
int status = integrate_dt (curr_time, tgo);
if (status != 0) {
return status;
if (tgo > 0.0) {
t_from = t_to;
t_to = t_to + tgo;
int status = integrate_dt( t_from, (t_to-t_from));
if (status != 0) { return status; }
} else {
#define FORWARD_INTEGRATION 0
#if FORWARD_INTEGRATION
integ_ptr->state_reset(); // We always want to integrate forward.
t_to = t_to + tgo;
#else
t_from = t_to;
t_to = t_from + tgo;
#endif
int status = integrate_dt( t_from, (t_to-t_from));
if (status != 0) { return status; }
}
// Refine the estimate of the time to the event time.
end_offset -= tgo;
curr_time += tgo;
tgo = curr_job->call_double();
}
}
}
// The last call to the event job should have changed derivatives
// in an event-dependent manner. Integrate back to the end of
// the integration cycle using the updated derivatives.
/*
Integrate to the end of the integration cycle using the updated derivatives.
*/
if (fired) {
return integrate_dt (curr_time, end_time-curr_time);
}
else {
t_from = t_to;
t_to = t_end;
/*
Because the derivatives have likely changed, we need to recursively
check for new events in the remaining interval.
*/
int status = integrate_dt( t_from, (t_to-t_from));
if (status != 0) { return status; }
status = process_dynamic_events(t_from, t_to, depth+1);
if (status != 0) { return status; }
} else {
return 0;
}
}

View File

@ -92,16 +92,11 @@ void Trick::Integrator::state_reset () {
#ifdef USE_ER7_UTILS_INTEGRATORS
#else
if (intermediate_step == 0) {
int i = 0;
double* next_arg = state_origin[i];
if (verbosity) std::cout << "RESET STATE: \n";
while (next_arg != (double*) NULL) {
*next_arg = state[i];
if (verbosity) std::cout << " " << *next_arg << "\n";
i++;
next_arg = state_origin[i];
if (verbosity) message_publish(MSG_DEBUG, "STATE RESET");
for (int i=0; i<num_state ; i++) {
*state_origin[i] = state[i];
}
if (verbosity) std::cout << std::endl;
}
#endif
}
@ -137,6 +132,14 @@ void Trick::Integrator::state_in (double* arg1, va_list argp) {
if (verbosity) std::cout << std::endl;
}
}
void Trick::Integrator::state_element_in (unsigned int index , double* state_p) {
if (intermediate_step == 0) {
state_origin[index] = state_p;
state[index] = *state_p;
}
}
#endif
void Trick::Integrator::state_in (double* arg1, ...) {

View File

@ -57,6 +57,18 @@ extern "C" void load_indexed_state(unsigned int index , double state) {
}
}
#ifndef USE_ER7_UTILS_INTEGRATORS
// Warning: state_p should never point to an automatic local variable.
extern "C" void load_state_element(unsigned int index , double* state_p) {
if (trick_curr_integ != NULL) {
trick_curr_integ->state_element_in (index, state_p);
} else {
message_publish(MSG_ERROR, "Integ load_indexed_state ERROR: trick_curr_integ is not set.\n") ;
}
}
#endif
extern "C" void load_deriv( double* arg1, ...) {
va_list argp;

View File

@ -10,7 +10,9 @@ double regula_falsi(
double time, /* In: Current time - input (seconds) */
REGULA_FALSI * R) /* Inout: Regula Falsi parameters */
{
if (R->iterations > 0 &&
if (
R->iterations > 0 &&
((M_ABS(R->error) < R->error_tol) ||
(M_ABS(R->last_error - R->error) < R->error_tol)))
{
@ -25,6 +27,7 @@ double regula_falsi(
return (0.0);
}
}
if (R->error < 0.0) {
/* Set lower bounds */
R->x_lower = R->error;
@ -40,16 +43,17 @@ double regula_falsi(
/* Have now got upper and lower bounds of zero point */
if (R->upper_set == 1 && R->lower_set == 1) {
/* Calculate time to error function zero point */
if (M_ABS(R->error) < R->error_tol)
if (M_ABS(R->error) < R->error_tol) {
R->delta_time = 0.0;
else {
R->delta_time = -R->error /
((R->x_upper - R->x_lower) / (R->t_upper -
R->t_lower));
if (R->iterations > 20)
} else {
double slope = (R->x_upper - R->x_lower) / (R->t_upper - R->t_lower);
R->delta_time = -R->error / slope;
if (R->iterations > 20) {
R->delta_time = 0.0;
}
}
/* Now check for increasing or decreasing constraint */
switch (R->mode) {
@ -58,35 +62,25 @@ double regula_falsi(
R->last_tgo = R->delta_time;
return (R->delta_time);
case Increasing: /* Increasing slope */
if (R->function_slope == Increasing) {
R->last_error = R->error;
R->last_tgo = R->delta_time;
// If the slope of the error function is positive.
if (R->t_upper > R->t_lower) {
return (R->delta_time);
} else {
R->lower_set = 0;
R->function_slope = Increasing;
}
break;
case Decreasing: /* Decreasing slope */
if (R->function_slope == Decreasing) {
R->last_error = R->error;
R->last_tgo = R->delta_time;
// If the slope of the error function is negative.
if (R->t_upper < R->t_lower) {
return (R->delta_time);
} else {
R->upper_set = 0;
R->function_slope = Decreasing;
}
break;
}
R->function_slope = Any;
} else if (R->lower_set == 1) {
R->function_slope = Increasing;
} else if (R->upper_set == 1) {
R->function_slope = Decreasing;
}
R->iterations = 0;
R->last_tgo = BIG_TGO;
return (BIG_TGO);
}

View File

@ -12,7 +12,9 @@ void Trick::ABM_Integrator::initialize(int State_size, double Dt) {
num_state = State_size;
state_origin = INTEG_ALLOC( double*, num_state );
state_origin[0] = (double*)NULL;
for(i=0; i<num_state ; i++) {
state_origin[i] = (double*)NULL;
}
/** Allocate the state vector.*/
state = INTEG_ALLOC( double, num_state );

View File

@ -14,7 +14,9 @@ void Trick::Euler_Cromer_Integrator::initialize(int State_size, double Dt) {
num_state = State_size;
state_origin = INTEG_ALLOC( double*, num_state );
state_origin[0] = (double*)NULL;
for(i=0; i<num_state ; i++) {
state_origin[i] = (double*)NULL;
}
/** Allocate the state vector.*/
state = INTEG_ALLOC( double, num_state );

View File

@ -16,7 +16,9 @@ void Trick::Euler_Integrator::initialize(int State_size, double Dt) {
num_state = State_size;
state_origin = INTEG_ALLOC( double*, num_state );
state_origin[0] = (double*)NULL;
for(int i=0; i<num_state ; i++) {
state_origin[i] = (double*)NULL;
}
/** Allocate the state vector.*/
state = INTEG_ALLOC( double, num_state );

View File

@ -13,7 +13,9 @@ void Trick::MM4_Integrator::initialize(int State_size, double Dt) {
num_state = State_size;
state_origin = INTEG_ALLOC( double*, num_state );
state_origin[0] = (double*)NULL;
for(i=0; i<num_state ; i++) {
state_origin[i] = (double*)NULL;
}
/** Allocate the state vector.*/
state = INTEG_ALLOC( double, num_state );

View File

@ -14,7 +14,9 @@ void Trick::NL2_Integrator::initialize(int State_size, double Dt) {
num_state = State_size;
state_origin = INTEG_ALLOC( double*, num_state );
state_origin[0] = (double*)NULL;
for(i=0; i<num_state ; i++) {
state_origin[i] = (double*)NULL;
}
/** Allocate the state vector.*/
state = INTEG_ALLOC( double, num_state );

View File

@ -15,7 +15,9 @@ void Trick::RK2_Integrator::initialize(int State_size, double Dt) {
dt = Dt;
state_origin = INTEG_ALLOC( double*, num_state );
state_origin[0] = (double*)NULL;
for(i=0; i<num_state ; i++) {
state_origin[i] = (double*)NULL;
}
/** Allocate the state vector.*/
state = INTEG_ALLOC( double, num_state );

View File

@ -13,8 +13,9 @@ void Trick::RK4_Integrator::initialize(int State_size, double Dt) {
num_state = State_size;
state_origin = INTEG_ALLOC( double*, num_state );
state_origin[0] = (double*)NULL;
for(i=0; i<num_state ; i++) {
state_origin[i] = (double*)NULL;
}
/** Allocate the state vector.*/
state = INTEG_ALLOC( double, num_state );

View File

@ -13,7 +13,9 @@ void Trick::RKF45_Integrator::initialize(int State_size, double Dt) {
num_state = State_size;
state_origin = INTEG_ALLOC( double*, num_state );
state_origin[0] = (double*)NULL;
for(i=0; i<num_state ; i++) {
state_origin[i] = (double*)NULL;
}
/** Allocate the state vector.*/
state = INTEG_ALLOC( double, num_state );

View File

@ -13,7 +13,9 @@ void Trick::RKF78_Integrator::initialize(int State_size, double Dt) {
num_state = State_size;
state_origin = INTEG_ALLOC( double*, num_state );
state_origin[0] = (double*)NULL;
for(i=0; i<num_state ; i++) {
state_origin[i] = (double*)NULL;
}
/** Allocate the state vector.*/
state = INTEG_ALLOC( double, num_state );

View File

@ -13,7 +13,9 @@ void Trick::RKG4_Integrator::initialize(int State_size, double Dt) {
num_state = State_size;
state_origin = INTEG_ALLOC( double*, num_state );
state_origin[0] = (double*)NULL;
for(i=0; i<num_state ; i++) {
state_origin[i] = (double*)NULL;
}
/** Allocate the state vector.*/
state = INTEG_ALLOC( double, num_state );