I'm writing a C++ code to solve the Lorentz force ODE, and I'm using the GSL library. I set an integer as the number of particles to be simulated and solved by the ODE, and print the states (x,y,z,vx,vy,vz) at each step. The initial values of position and velocity values are configured as random values from a GSL library. The code works for 1 particle, but returns "ERROR: integration limits and/or step direction not consistent" for any number of primaries > 1. I have tried different options, but haven't so far figured out what's wrong with my script. Below is the full code. I'd really appreciate any feedback on how to overcome this. Thanks a lot in advance.
Br, Leonardo.
#include<iostream>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
using namespace std;
int lorentzODE (double t, const double s[], double f[], void *params);
int main(){
int noPrimaries = 2;
double q = 1.0; // Particle charge
double m = 1.0; // Particle mass
double B[3] = {0.0, 0.0, 1.0}; // Magnetic field
double E[3] = {0.1, 0.0, 0.5}; // Electric field
double x0[3]; // Initial position
double v0[3]; // Initial velocity
const gsl_rng_type * T = gsl_rng_ranlxs0;
gsl_rng * r = gsl_rng_alloc(T);
gsl_rng_set(r, (unsigned long) time(NULL));
double t0 = 0.0;
double tf = 100.0;
double dt = 0.05;
double system[6]; // Initial state vector
int status; // status of driver function
double paramsB[8];
paramsB[0] = q;
paramsB[1] = m;
paramsB[2] = B[0];
paramsB[3] = B[1];
paramsB[4] = B[2];
paramsB[5] = E[0];
paramsB[6] = E[1];
paramsB[7] = E[2];
const string &solverMethod = "RK8PD";
const double h = 1.0e-06;
const double epsAbs = 1.0e-08;
const double epsRel = 1.0e-10;
// Integrator configuration
double t, t_next;
gsl_odeiv2_system odeSystem;
odeSystem.function = lorentzODE;
odeSystem.dimension = 6;
odeSystem.params = paramsB;
// Algorithm option
gsl_odeiv2_driver *drv;
drv = gsl_odeiv2_driver_alloc_y_new(&odeSystem, gsl_odeiv2_step_rk8pd, h, epsAbs, epsRel);
for(int i=0;i<noPrimaries;i++){
x0[0] = gsl_rng_uniform (r);
x0[1] = gsl_rng_uniform (r);
x0[2] = gsl_rng_uniform (r);
v0[0] = gsl_rng_uniform (r);
v0[1] = gsl_rng_uniform (r);
v0[2] = gsl_rng_uniform (r);
system[0] = x0[0]; system[1] = x0[1]; system[2] = x0[2];
system[3] = v0[0]; system[4] = v0[1]; system[5] = v0[2];
for(t_next = t0+dt; t_next<tf; t_next += dt){
status = gsl_odeiv2_driver_apply(drv, &t, t_next, system);
if(status != GSL_SUCCESS){
printf("Error: status = %d\n", status);
break;
}
printf("%.5e %.5e %.5e %.5e %.5e %.5e %.5e\n", t, system[0], system[1], system[2], system[3], system[4], system[5]);
}
gsl_odeiv2_driver_free(drv);
}
}
int lorentzODE (double t, const double s[], double f[], void *params){
(void)(t); /* avoid unused parameter warning */
double *lparams = (double *)params;
double q = lparams[0];
double m = lparams[1];
double mu = q/m;
double Bx = lparams[2];
double By = lparams[3];
double Bz = lparams[4];
double Ex = lparams[5];
double Ey = lparams[6];
double Ez = lparams[7];
f[0] = s[3];
f[1] = s[4];
f[2] = s[5];
f[3] = mu*(Bz*s[4] - By*s[5] + Ex);
f[4] = mu*(Bx*s[5] - Bz*s[3] + Ey);
f[5] = mu*(By*s[3] - Bx*s[4] + Ez);
return GSL_SUCCESS;
}
The problem was solved :) There were two reasons for why it wasn't working:
t. Corrected by initializing it inside the first loop;drv has to be performed inside the loop as well.Thanks Lutz Lehmann for pointing it out. The below is the working code:
#include<iostream>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
using namespace std;
int lorentzODE (double t, const double s[], double f[], void *params);
int main(){
int noPrimaries = 10;
double q = 1.0; // Particle charge
double m = 1.0; // Particle mass
double B[3] = {0.0, 0.0, 1.0}; // Magnetic field
double E[3] = {0.1, 0.0, 0.5}; // Electric field
double x0[3]; // Initial position
double v0[3]; // Initial velocity
const gsl_rng_type * T = gsl_rng_ranlxs0;
gsl_rng * r = gsl_rng_alloc(T);
gsl_rng_set(r, (unsigned long) time(NULL));
double t0 = 0.0;
double tf = 10.0;
double dt = 0.5;
double system[6]; // Initial state vector
int status; // status of driver function
double paramsB[8];
paramsB[0] = q;
paramsB[1] = m;
paramsB[2] = B[0];
paramsB[3] = B[1];
paramsB[4] = B[2];
paramsB[5] = E[0];
paramsB[6] = E[1];
paramsB[7] = E[2];
const double h = 1.0e-06;
const double epsAbs = 1.0e-08;
const double epsRel = 1.0e-10;
double t, t_next;
gsl_odeiv2_system odeSystem;
odeSystem.function = lorentzODE;
odeSystem.dimension = 6;
odeSystem.params = paramsB;
gsl_odeiv2_driver *drv;
for(int i=0;i<noPrimaries;i++){
t = t0;
drv = gsl_odeiv2_driver_alloc_y_new(&odeSystem, gsl_odeiv2_step_rk8pd, h, epsAbs, epsRel);
x0[0] = gsl_rng_uniform (r);
x0[1] = gsl_rng_uniform (r);
x0[2] = gsl_rng_uniform (r);
v0[0] = gsl_rng_uniform (r);
v0[1] = gsl_rng_uniform (r);
v0[2] = gsl_rng_uniform (r);
system[0] = x0[0]; system[1] = x0[1]; system[2] = x0[2];
system[3] = v0[0]; system[4] = v0[1]; system[5] = v0[2];
for(t_next = t0+dt; t_next<tf; t_next += dt){
status = gsl_odeiv2_driver_apply(drv, &t, t_next, system);
if(status != GSL_SUCCESS){
printf("Error: status = %d\n", status);
break;
}
printf("%.5e %.5e %.5e %.5e %.5e %.5e %.5e\n", t, system[0], system[1], system[2], system[3], system[4], system[5]);
}
gsl_odeiv2_driver_free(drv);
}
}
int lorentzODE (double t, const double s[], double f[], void *params){
(void)(t); /* avoid unused parameter warning */
double *lparams = (double *)params;
double q = lparams[0];
double m = lparams[1];
double mu = q/m;
double Bx = lparams[2];
double By = lparams[3];
double Bz = lparams[4];
double Ex = lparams[5];
double Ey = lparams[6];
double Ez = lparams[7];
f[0] = s[3];
f[1] = s[4];
f[2] = s[5];
f[3] = mu*(Bz*s[4] - By*s[5] + Ex);
f[4] = mu*(Bx*s[5] - Bz*s[3] + Ey);
f[5] = mu*(By*s[3] - Bx*s[4] + Ez);
return GSL_SUCCESS;
}