mardi 10 septembre 2019

How to pass a member of a class to gsl with many parameters

I'm trying to solve a system of ordinary differential equations with gsl. Everything should happen within a class, so the function that is passed to the gsl system is a member of that class. I have found some successful suggestions about how to wrap the functions into a static cast. This works perfectly for the function itself, but all the parameters that are initialised when I instantiate the class are lost and the gsl system does not see them, and results in segmentation fault.

I've found useful tips in

how to avoid static member function when using gsl with c++

Calling GSL function inside a class in a shared library

how to avoid static member function when using gsl with c++

That's where I learnt how to wrap the rhs and jac. Like I said this wrapping works well, but the parameters that are initialised by the constructor gmf() do not get passed along and the ode solver fails.

Here's a minimal header file:

#ifndef TEST_H_
#define TEST_H_

#include <vector>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

//integration parameters
#define STEPSIZE 1e-3 //1e-3
#define ERR_ABS 1e-6 //1e-6
#define ERR_REL 1e-10 //1e-10
#define ENDPOINT 100.0 //100.0

//gmf
class gmf {
    double a, b, c;

public:
    gmf();

    std::vector<double> value(const double r[]);
    int rhs(const double t, const double y[], double fy[]);
    static int rhs_wrap(const double t, const double y[], double fy[], void * wrap);
    int jac(const double t, const double y[], double *dfdy, double dfdt[]);
    static int jac_wrap(const double t, const double y[], double *dfdy, double dfdt[], void * wrap);
    void original_direction(double l, double b);
};

#endif

that goes along with the implementation

#include <iostream>
#include <cmath>

#include "test.h"

//constructor
gmf::gmf() {

    a = 1.;
    b = 1.;
    c = 1.;
}

//field value
std::vector<double> gmf::value(const double r[]) {

    //I shouldn't have to explicitly overwrite these parameters here, but I don't know how to pass them properly
//  double a = 1.;
//  double b = 1.;
//  double c = 1.;
    std::vector<double> result(3);
    double rx = r[0];
    double ry = r[1];
    double rz = r[2];
    result[0] = a*rx;
    result[1] = b*ry;
    result[1] = c*rz;
    return result;
}

//integration
int gmf::rhs(const double t, const double y[], double fy[]){

    std::vector<double> result(3);
    result = this->value(y);
    fy[0] = y[3];
    fy[1] = y[4];
    fy[2] = y[5];
    fy[3] = (y[4]*result[2]-y[5]*result[1]);
    fy[4] = (y[5]*result[0]-y[3]*result[2]);
    fy[5] = (y[3]*result[1]-y[4]*result[0]);
    return GSL_SUCCESS;
}

int gmf::rhs_wrap(const double t, const double y[], double fy[], void * wrap) {
        return static_cast<gmf*>(wrap)->rhs(t,y,fy);
}

int gmf::jac(const double t, const double y[], double *dfdy, double dfdt[]) {

    gsl_matrix_view dfdy_mat = gsl_matrix_view_array(dfdy,6,6);
    gsl_matrix * m = &dfdy_mat.matrix; 
    gsl_matrix_set_zero(m);
    dfdt[0] = 0.;
    dfdt[1] = 0.;
    dfdt[2] = 0.;
    dfdt[3] = 0.;
    dfdt[4] = 0.;
    dfdt[5] = 0.;
    return GSL_SUCCESS;
}

int gmf::jac_wrap(const double t, const double y[], double *dfdy, double dfdt[], void * wrap) {
        return static_cast<gmf*>(wrap)->jac(t,y,dfdy,dfdt);
}

void gmf::original_direction(double l, double b) {

    double y[6] = {0.,0.,0.,cos(b)*cos(l),cos(b)*sin(l),sin(b)};
    //double h = STEPSIZE;
    double eps_abs = ERR_ABS;
    double eps_rel = ERR_REL;
    gsl_odeiv2_system sys = {rhs_wrap,jac_wrap,6}; //we need to wrap these up in static to run from within class
    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new(&sys,gsl_odeiv2_step_rk8pd,eps_abs,eps_rel,0.);
    int iter;
    double t = 0., tend = 100.;
    int nsteps = 100;
    for (iter=1; iter<=nsteps; iter++) {
        double t1 = iter*tend/((double) nsteps);
        int status = gsl_odeiv2_driver_apply(d,&t,t1,y);
        if (status != GSL_SUCCESS) {
            std::cout << "  gmf::original_direction: gsl integration failed\n";
            exit(EXIT_FAILURE);
        }
    }
    gsl_odeiv2_driver_free(d);
    double xv = y[3];
    double yv = y[4];
    double zv = y[5];
    std::cout<<xv<<"\t"<<yv<<"\t"<<zv<<"\n";
}

/////////////////////////
//MAIN
/////////////////////////

//main engine
int main(void){
    std::cout<<"\n";

    gmf myGMF; //here I should have already instantiated the values of a,b,c, but somehow they don't get passed to the gsl system
    myGMF.original_direction(0.1,0.1);

    return EXIT_SUCCESS;
}

Like I said, when I instantiate gmf myGMF in main(void) this object automatically has the correct values for the parameters a,b,c (I checked). But when I use the member gmf::original_direction for the gsl system these parameters are lost and I get segmentation fault.

If however I explicitly declare these parameters inside the member function gmf::value (in the snippet above they are commented out) then everything goes as it should.

So, the question is: how do I make sure to pass everything from the instantiation of the class to the gsl?

In the links above I have seen people suggesting the use of lambdas to collect the pointer to the member function and all the parameters at once, but I don't know lambdas enough to understand their suggestions.

Thanks!

Aucun commentaire:

Enregistrer un commentaire