Commit cf9c803e authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

Merge branch 'feature/rFFTW-distributed-particles' into develop

parents d4227289 38717d96
...@@ -39,7 +39,7 @@ ...@@ -39,7 +39,7 @@
extern int myrank, nprocs; extern int myrank, nprocs;
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
distributed_particles<particle_type, rnumber, interp_neighbours>::distributed_particles( distributed_particles<particle_type, rnumber, interp_neighbours>::distributed_particles(
const char *NAME, const char *NAME,
const hid_t data_file_id, const hid_t data_file_id,
...@@ -56,21 +56,25 @@ distributed_particles<particle_type, rnumber, interp_neighbours>::distributed_pa ...@@ -56,21 +56,25 @@ distributed_particles<particle_type, rnumber, interp_neighbours>::distributed_pa
this->vel = FIELD; this->vel = FIELD;
this->rhs.resize(INTEGRATION_STEPS); this->rhs.resize(INTEGRATION_STEPS);
this->integration_steps = INTEGRATION_STEPS; this->integration_steps = INTEGRATION_STEPS;
this->state.reserve(2*this->nparticles / this->nprocs);
for (unsigned int i=0; i<this->rhs.size(); i++)
this->rhs[i].reserve(2*this->nparticles / this->nprocs);
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
distributed_particles<particle_type, rnumber, interp_neighbours>::~distributed_particles() distributed_particles<particle_type, rnumber, interp_neighbours>::~distributed_particles()
{ {
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
void distributed_particles<particle_type, rnumber, interp_neighbours>::sample( void distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
interpolator<rnumber, interp_neighbours> *field, interpolator<rnumber, interp_neighbours> *field,
const std::unordered_map<int, single_particle_state<particle_type>> &x,
std::unordered_map<int, single_particle_state<POINT3D>> &y) std::unordered_map<int, single_particle_state<POINT3D>> &y)
{ {
double *yy = new double[3]; double *yy = new double[3];
y.clear(); y.clear();
for (auto &pp: this->state) for (auto &pp: x)
{ {
(*field)(pp.second.data, yy); (*field)(pp.second.data, yy);
y[pp.first] = yy; y[pp.first] = yy;
...@@ -78,39 +82,41 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::sample( ...@@ -78,39 +82,41 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
delete[] yy; delete[] yy;
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
void distributed_particles<particle_type, rnumber, interp_neighbours>::get_rhs( void distributed_particles<particle_type, rnumber, interp_neighbours>::get_rhs(
const std::unordered_map<int, single_particle_state<particle_type>> &x, const std::unordered_map<int, single_particle_state<particle_type>> &x,
std::unordered_map<int, single_particle_state<particle_type>> &y) std::unordered_map<int, single_particle_state<particle_type>> &y)
{ {
double *yy = new double[this->ncomponents]; std::unordered_map<int, single_particle_state<POINT3D>> yy;
switch(particle_type)
{
case VELOCITY_TRACER:
this->sample(this->vel, this->state, yy);
y.clear(); y.clear();
for (auto &pp: x) for (auto &pp: x)
{ y[pp.first] = yy[pp.first].data;
(*this->vel)(pp.second.data, yy); break;
y[pp.first] = yy;
} }
delete[] yy;
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
void distributed_particles<particle_type, rnumber, interp_neighbours>::sample( void distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
interpolator<rnumber, interp_neighbours> *field, interpolator<rnumber, interp_neighbours> *field,
const char *dset_name) const char *dset_name)
{ {
std::unordered_map<int, single_particle_state<POINT3D>> y; std::unordered_map<int, single_particle_state<POINT3D>> y;
this->sample(field, y); this->sample(field, this->state, y);
this->write(dset_name, y); this->write(dset_name, y);
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
void distributed_particles<particle_type, rnumber, interp_neighbours>::roll_rhs() void distributed_particles<particle_type, rnumber, interp_neighbours>::roll_rhs()
{ {
for (int i=this->integration_steps-2; i>=0; i--) for (int i=this->integration_steps-2; i>=0; i--)
rhs[i+1] = rhs[i]; rhs[i+1] = rhs[i];
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
void distributed_particles<particle_type, rnumber, interp_neighbours>::redistribute( void distributed_particles<particle_type, rnumber, interp_neighbours>::redistribute(
std::unordered_map<int, single_particle_state<particle_type>> &x, std::unordered_map<int, single_particle_state<particle_type>> &x,
std::vector<std::unordered_map<int, single_particle_state<particle_type>>> &vals) std::vector<std::unordered_map<int, single_particle_state<particle_type>>> &vals)
...@@ -168,7 +174,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib ...@@ -168,7 +174,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib
buffer_size = (buffer_size > npr[0])? buffer_size : npr[0]; buffer_size = (buffer_size > npr[0])? buffer_size : npr[0];
buffer_size = (buffer_size > npr[1])? buffer_size : npr[1]; buffer_size = (buffer_size > npr[1])? buffer_size : npr[1];
//DEBUG_MSG("buffer size is %d\n", buffer_size); //DEBUG_MSG("buffer size is %d\n", buffer_size);
double *buffer = new double[buffer_size*this->ncomponents*(1+vals.size())]; double *buffer = new double[buffer_size*state_dimension(particle_type)*(1+vals.size())];
for (rsrc = 0; rsrc<this->nprocs; rsrc++) for (rsrc = 0; rsrc<this->nprocs; rsrc++)
for (int i=0; i<2; i++) for (int i=0; i<2; i++)
{ {
...@@ -186,21 +192,21 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib ...@@ -186,21 +192,21 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib
for (int p: ps[i]) for (int p: ps[i])
{ {
std::copy(x[p].data, std::copy(x[p].data,
x[p].data + this->ncomponents, x[p].data + state_dimension(particle_type),
buffer + pcounter*(1+vals.size())*this->ncomponents); buffer + pcounter*(1+vals.size())*state_dimension(particle_type));
x.erase(p); x.erase(p);
for (int tindex=0; tindex<vals.size(); tindex++) for (int tindex=0; tindex<vals.size(); tindex++)
{ {
std::copy(vals[tindex][p].data, std::copy(vals[tindex][p].data,
vals[tindex][p].data + this->ncomponents, vals[tindex][p].data + state_dimension(particle_type),
buffer + (pcounter*(1+vals.size()) + tindex+1)*this->ncomponents); buffer + (pcounter*(1+vals.size()) + tindex+1)*state_dimension(particle_type));
vals[tindex].erase(p); vals[tindex].erase(p);
} }
pcounter++; pcounter++;
} }
MPI_Send( MPI_Send(
buffer, buffer,
nps[i]*(1+vals.size())*this->ncomponents, nps[i]*(1+vals.size())*state_dimension(particle_type),
MPI_DOUBLE, MPI_DOUBLE,
rdst, rdst,
2*(rsrc*this->nprocs + rdst)+1, 2*(rsrc*this->nprocs + rdst)+1,
...@@ -218,7 +224,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib ...@@ -218,7 +224,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib
MPI_STATUS_IGNORE); MPI_STATUS_IGNORE);
MPI_Recv( MPI_Recv(
buffer, buffer,
npr[1-i]*(1+vals.size())*this->ncomponents, npr[1-i]*(1+vals.size())*state_dimension(particle_type),
MPI_DOUBLE, MPI_DOUBLE,
rsrc, rsrc,
2*(rsrc*this->nprocs + rdst)+1, 2*(rsrc*this->nprocs + rdst)+1,
...@@ -227,10 +233,10 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib ...@@ -227,10 +233,10 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib
int pcounter = 0; int pcounter = 0;
for (int p: pr[1-i]) for (int p: pr[1-i])
{ {
x[p] = buffer + (pcounter*(1+vals.size()))*this->ncomponents; x[p] = buffer + (pcounter*(1+vals.size()))*state_dimension(particle_type);
for (int tindex=0; tindex<vals.size(); tindex++) for (int tindex=0; tindex<vals.size(); tindex++)
{ {
vals[tindex][p] = buffer + (pcounter*(1+vals.size()) + tindex+1)*this->ncomponents; vals[tindex][p] = buffer + (pcounter*(1+vals.size()) + tindex+1)*state_dimension(particle_type);
} }
pcounter++; pcounter++;
} }
...@@ -255,13 +261,13 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib ...@@ -255,13 +261,13 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
void distributed_particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth( void distributed_particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
const int nsteps) const int nsteps)
{ {
this->get_rhs(this->state, this->rhs[0]); this->get_rhs(this->state, this->rhs[0]);
for (auto &pp: this->state) for (auto &pp: this->state)
for (int i=0; i<this->ncomponents; i++) for (int i=0; i<state_dimension(particle_type); i++)
switch(nsteps) switch(nsteps)
{ {
case 1: case 1:
...@@ -303,7 +309,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::AdamsBash ...@@ -303,7 +309,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::AdamsBash
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
void distributed_particles<particle_type, rnumber, interp_neighbours>::step() void distributed_particles<particle_type, rnumber, interp_neighbours>::step()
{ {
this->AdamsBashforth((this->iteration < this->integration_steps) ? this->AdamsBashforth((this->iteration < this->integration_steps) ?
...@@ -313,10 +319,10 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::step() ...@@ -313,10 +319,10 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::step()
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
void distributed_particles<particle_type, rnumber, interp_neighbours>::read() void distributed_particles<particle_type, rnumber, interp_neighbours>::read()
{ {
double *temp = new double[this->chunk_size*this->ncomponents]; double *temp = new double[this->chunk_size*state_dimension(particle_type)];
for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++) for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
{ {
//read state //read state
...@@ -324,14 +330,14 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::read() ...@@ -324,14 +330,14 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::read()
this->read_state_chunk(cindex, temp); this->read_state_chunk(cindex, temp);
MPI_Bcast( MPI_Bcast(
temp, temp,
this->chunk_size*this->ncomponents, this->chunk_size*state_dimension(particle_type),
MPI_DOUBLE, MPI_DOUBLE,
0, 0,
this->comm); this->comm);
for (int p=0; p<this->chunk_size; p++) for (int p=0; p<this->chunk_size; p++)
{ {
if (this->vel->get_rank(temp[this->ncomponents*p+2]) == this->myrank) if (this->vel->get_rank(temp[state_dimension(particle_type)*p+2]) == this->myrank)
this->state[p+cindex*this->chunk_size] = temp + this->ncomponents*p; this->state[p+cindex*this->chunk_size] = temp + state_dimension(particle_type)*p;
} }
//read rhs //read rhs
if (this->iteration > 0) if (this->iteration > 0)
...@@ -341,7 +347,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::read() ...@@ -341,7 +347,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::read()
this->read_rhs_chunk(cindex, i, temp); this->read_rhs_chunk(cindex, i, temp);
MPI_Bcast( MPI_Bcast(
temp, temp,
this->chunk_size*this->ncomponents, this->chunk_size*state_dimension(particle_type),
MPI_DOUBLE, MPI_DOUBLE,
0, 0,
this->comm); this->comm);
...@@ -349,7 +355,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::read() ...@@ -349,7 +355,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::read()
{ {
auto pp = this->state.find(p+cindex*this->chunk_size); auto pp = this->state.find(p+cindex*this->chunk_size);
if (pp != this->state.end()) if (pp != this->state.end())
this->rhs[i][p+cindex*this->chunk_size] = temp + this->ncomponents*p; this->rhs[i][p+cindex*this->chunk_size] = temp + state_dimension(particle_type)*p;
} }
} }
} }
...@@ -357,7 +363,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::read() ...@@ -357,7 +363,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::read()
delete[] temp; delete[] temp;
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
void distributed_particles<particle_type, rnumber, interp_neighbours>::write( void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
const char *dset_name, const char *dset_name,
std::unordered_map<int, single_particle_state<POINT3D>> &y) std::unordered_map<int, single_particle_state<POINT3D>> &y)
...@@ -389,28 +395,28 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::write( ...@@ -389,28 +395,28 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
delete[] data; delete[] data;
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
void distributed_particles<particle_type, rnumber, interp_neighbours>::write( void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
const bool write_rhs) const bool write_rhs)
{ {
double *temp0 = new double[this->chunk_size*this->ncomponents]; double *temp0 = new double[this->chunk_size*state_dimension(particle_type)];
double *temp1 = new double[this->chunk_size*this->ncomponents]; double *temp1 = new double[this->chunk_size*state_dimension(particle_type)];
for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++) for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
{ {
//write state //write state
std::fill_n(temp0, this->ncomponents*this->chunk_size, 0); std::fill_n(temp0, state_dimension(particle_type)*this->chunk_size, 0);
for (int p=0; p<this->chunk_size; p++) for (int p=0; p<this->chunk_size; p++)
{ {
auto pp = this->state.find(p + cindex*this->chunk_size); auto pp = this->state.find(p + cindex*this->chunk_size);
if (pp != this->state.end()) if (pp != this->state.end())
std::copy(pp->second.data, std::copy(pp->second.data,
pp->second.data + this->ncomponents, pp->second.data + state_dimension(particle_type),
temp0 + pp->first*this->ncomponents); temp0 + pp->first*state_dimension(particle_type));
} }
MPI_Allreduce( MPI_Allreduce(
temp0, temp0,
temp1, temp1,
this->ncomponents*this->chunk_size, state_dimension(particle_type)*this->chunk_size,
MPI_DOUBLE, MPI_DOUBLE,
MPI_SUM, MPI_SUM,
this->comm); this->comm);
...@@ -420,19 +426,19 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::write( ...@@ -420,19 +426,19 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
if (write_rhs) if (write_rhs)
for (int i=0; i<this->integration_steps; i++) for (int i=0; i<this->integration_steps; i++)
{ {
std::fill_n(temp0, this->ncomponents*this->chunk_size, 0); std::fill_n(temp0, state_dimension(particle_type)*this->chunk_size, 0);
for (int p=0; p<this->chunk_size; p++) for (int p=0; p<this->chunk_size; p++)
{ {
auto pp = this->rhs[i].find(p + cindex*this->chunk_size); auto pp = this->rhs[i].find(p + cindex*this->chunk_size);
if (pp != this->rhs[i].end()) if (pp != this->rhs[i].end())
std::copy(pp->second.data, std::copy(pp->second.data,
pp->second.data + this->ncomponents, pp->second.data + state_dimension(particle_type),
temp0 + pp->first*this->ncomponents); temp0 + pp->first*state_dimension(particle_type));
} }
MPI_Allreduce( MPI_Allreduce(
temp0, temp0,
temp1, temp1,
this->ncomponents*this->chunk_size, state_dimension(particle_type)*this->chunk_size,
MPI_DOUBLE, MPI_DOUBLE,
MPI_SUM, MPI_SUM,
this->comm); this->comm);
......
...@@ -39,7 +39,7 @@ ...@@ -39,7 +39,7 @@
#define DISTRIBUTED_PARTICLES #define DISTRIBUTED_PARTICLES
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
class distributed_particles: public particles_io_base<particle_type> class distributed_particles: public particles_io_base<particle_type>
{ {
private: private:
...@@ -74,6 +74,7 @@ class distributed_particles: public particles_io_base<particle_type> ...@@ -74,6 +74,7 @@ class distributed_particles: public particles_io_base<particle_type>
const char *dset_name); const char *dset_name);
void sample( void sample(
interpolator<rnumber, interp_neighbours> *field, interpolator<rnumber, interp_neighbours> *field,
const std::unordered_map<int, single_particle_state<particle_type>> &x,
std::unordered_map<int, single_particle_state<POINT3D>> &y); std::unordered_map<int, single_particle_state<POINT3D>> &y);
void get_rhs( void get_rhs(
const std::unordered_map<int, single_particle_state<particle_type>> &x, const std::unordered_map<int, single_particle_state<particle_type>> &x,
......
...@@ -42,6 +42,7 @@ class interpolator:public interpolator_base<rnumber, interp_neighbours> ...@@ -42,6 +42,7 @@ class interpolator:public interpolator_base<rnumber, interp_neighbours>
rnumber *field; rnumber *field;
public: public:
using interpolator_base<rnumber, interp_neighbours>::operator();
ptrdiff_t buffer_size; ptrdiff_t buffer_size;
/* descriptor for buffered field */ /* descriptor for buffered field */
...@@ -67,17 +68,6 @@ class interpolator:public interpolator_base<rnumber, interp_neighbours> ...@@ -67,17 +68,6 @@ class interpolator:public interpolator_base<rnumber, interp_neighbours>
const double *__restrict__ x, const double *__restrict__ x,
double *__restrict__ y, double *__restrict__ y,
const int *deriv = NULL); const int *deriv = NULL);
/* interpolate 1 point */
inline void operator()(
const double *__restrict__ x,
double *__restrict__ dest,
const int *deriv = NULL)
{
int xg[3];
double xx[3];
this->get_grid_coordinates(x, xg, xx);
(*this)(xg, xx, dest, deriv);
}
void operator()( void operator()(
const int *__restrict__ xg, const int *__restrict__ xg,
const double *__restrict__ xx, const double *__restrict__ xx,
......
...@@ -87,6 +87,18 @@ class interpolator_base ...@@ -87,6 +87,18 @@ class interpolator_base
const double *__restrict__ xx, const double *__restrict__ xx,
double *__restrict__ dest, double *__restrict__ dest,
const int *deriv = NULL) = 0; const int *deriv = NULL) = 0;
/* interpolate 1 point */
inline void operator()(
const double *__restrict__ x,
double *__restrict__ dest,
const int *deriv = NULL)
{
int xg[3];
double xx[3];
this->get_grid_coordinates(x, xg, xx);
(*this)(xg, xx, dest, deriv);
}
}; };
#endif//INTERPOLATOR_BASE #endif//INTERPOLATOR_BASE
......
...@@ -39,7 +39,7 @@ ...@@ -39,7 +39,7 @@
extern int myrank, nprocs; extern int myrank, nprocs;
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
particles<particle_type, rnumber, interp_neighbours>::particles( particles<particle_type, rnumber, interp_neighbours>::particles(
const char *NAME, const char *NAME,
const hid_t data_file_id, const hid_t data_file_id,
...@@ -55,7 +55,7 @@ particles<particle_type, rnumber, interp_neighbours>::particles( ...@@ -55,7 +55,7 @@ particles<particle_type, rnumber, interp_neighbours>::particles(
(INTEGRATION_STEPS >= 1)); (INTEGRATION_STEPS >= 1));
this->vel = FIELD; this->vel = FIELD;
this->integration_steps = INTEGRATION_STEPS; this->integration_steps = INTEGRATION_STEPS;
this->array_size = this->nparticles * this->ncomponents; this->array_size = this->nparticles * state_dimension(particle_type);
this->state = new double[this->array_size]; this->state = new double[this->array_size];
std::fill_n(this->state, this->array_size, 0.0); std::fill_n(this->state, this->array_size, 0.0);
for (int i=0; i < this->integration_steps; i++) for (int i=0; i < this->integration_steps; i++)
...@@ -65,7 +65,7 @@ particles<particle_type, rnumber, interp_neighbours>::particles( ...@@ -65,7 +65,7 @@ particles<particle_type, rnumber, interp_neighbours>::particles(
} }
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
particles<particle_type, rnumber, interp_neighbours>::~particles() particles<particle_type, rnumber, interp_neighbours>::~particles()
{ {
delete[] this->state; delete[] this->state;
...@@ -75,18 +75,18 @@ particles<particle_type, rnumber, interp_neighbours>::~particles() ...@@ -75,18 +75,18 @@ particles<particle_type, rnumber, interp_neighbours>::~particles()
} }
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
void particles<particle_type, rnumber, interp_neighbours>::get_rhs(double *x, double *y) void particles<particle_type, rnumber, interp_neighbours>::get_rhs(double *x, double *y)
{ {
switch(particle_type) switch(particle_type)
{ {
case VELOCITY_TRACER: case VELOCITY_TRACER:
this->vel->sample(this->nparticles, this->ncomponents, x, y); this->vel->sample(this->nparticles, state_dimension(particle_type), x, y);
break; break;
} }
} }
template <int particle_type, class rnumber, int interp_neighbours> template <particle_types particle_type, class rnumber, int interp_neighbours>
void particles<particle_type, rnumber, interp_neighbours>::roll_rhs() void particles<particle_type, rnumber, interp_neighbours>::roll_rhs()
{ {
for (int i=this->integration_steps-2; i>=0; i--) for (int i=this->integration_steps-2; i>=0; i--)