0

I am struggling with -- apparently -- a well known problem. I have an ODE system which is defined by a class member function, and I want to solve/integrate it by one of the GSL solvers. I.e., say my model is defined by:

class my_model{
...
public:
    int NEQ = 4;
    double y[4], dydt[4];
    double params[25]; 
    int ode_gsl(double t, const double y[], double dydt[], void * params);
    ...
};

int my_model::int ode_gsl(double t, const double y[], double dydt[], void * params){
    dydt[0] = params[1]*params[0]*y[0] - y[1];
    dydt[1] = ...;
    ...
    return GSL_SUCCESS;    
}

Then in my integration routine:

int main(){
    my_model chemosc;

   // Parameters for the Integrator
   double HSTART = 1e-3;
   double ATOL = 1e-8;
   double RTOL = 1e-6;

   // Instantiate GSL solver
   gsl_odeiv2_system sys = {&chemosc.ode_gsl, nullptr, chemosc.NEQ, chemosc.params};
   gsl_odeiv2_driver * d;
   d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, HSTART, ATOL, RTOL);

   double twin[2] = {0.,60.};
   double dt = 1e-3;
   double t = twin[0], t1 = twin[1];
   long int NSTEP = (long int)((t1-t)/dt)+1; // +1 if you start counting from zero...
   int NEQ = 4;
   long int NUMEL = (NEQ+1)*NSTEP; // number of elements for solution

   int i = 0,j;
   do{
      double ti = t + dt; // t is automatically updated by the driver
      printf("\n%.3f\t%.3f\t%.3f t%.3f",astro.y[0],astro.y[1],astro.y[2],astro.y[3]);
      int status = gsl_odeiv2_driver_apply (d, &t, ti, astro.y);
      ...
      }
...
}

Compiling the above code creates the somehow well known error that GSL requires a pointer-to-function whereas I am passing a pointer-to-member function, i.e.:

error: cannot convert ‘int (chemosc::*)(double, const double*, double*, void*)’ to ‘int (*)(double, const double*, double*, void*)’

I found several questions related to this topic: Q1, Q2, Q3, Q4, Q5, Q6, but seriously, the answer there are hard to follow. Declaring as static my member function, has the downside that the compiler requires me to also declare as static all the member parameters. Using static_cast as suggested here, result in all bunch of segmentation errors (but I assume that I did something wrong in the implementation, inasmuch as the directions in that post are quite cryptic). It would be nice to have a once-for-all clear and as simple as possible working solution for this problem (possibly without using boost libraries).

maurizio
  • 585
  • 6
  • 19
  • Function signature seem to contain `void *` parameter that you can use to pass a pointer to a class instance (and whatever parameters you need). So you can declare a static member function (or non-member function) then cast `void *` parameter to the class instance and then invoke class member function that will do all the computations. – user7860670 Nov 01 '17 at 08:31
  • @VTT Yes, I kinda figure out that this is the way to proceed. But could you please give me come exemplifying code? What you say was mentioned in https://stackoverflow.com/questions/10593726/a-function-pointer-issue-how-to-efficiently-interface-with-c-api-ie-gsl-from, but again, it is related to integrators, rather than ODE solvers. And the explanation remains somewhat confusing. – maurizio Nov 01 '17 at 08:36

1 Answers1

2

Something like this:

class my_model
{
    private: int
    ode_gsl_impl(double const t, double const * const y, double * const dydt);

    public: static int
    ode_gsl(double const t, double const * const y, double * const dydt, void * const opaque)
    {
        assert(opaque);
        return(static_cast<my_model *>(opaque)->ode_gsl_impl(t, y, dydt));
    }
};

gsl_odeiv2_system sys
{
    &my_model::ode_gsl
,   nullptr
,   chemosc.NEQ
,   reinterpret_cast<void *>(::std::addressof(chemosc))
};

I would like to mention that your original code suffers from name collisions between callback parameter names and class field names.

user7860670
  • 32,142
  • 4
  • 44
  • 75