Skip to content

Commit

Permalink
Merge pull request #191 from ExCALIBUR-NEPTUNE/bugfix/sycl-buffers-at…
Browse files Browse the repository at this point in the history 
…tempt-two

Bugfix/sycl buffers attempt two
  • Loading branch information
@JamesEdgeley
JamesEdgeley authored Apr 10, 2024
2 parents 40c3dfe + 1cdc5c3 commit b78aac7
Show file tree
Hide file tree
Showing 7 changed files with 226 additions and 230 deletions.
83 changes: 39 additions & 44 deletions include/mesh.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,17 +28,13 @@ class Mesh {
int nintervals;
// number of grid points (including periodic point)
int nmesh;
sycl::buffer<int, 1> nmesh_d;
// grid spacing
double dx;
sycl::buffer<double, 1> dx_d;

// box length in units of Debye length
double normalized_box_length;
// mesh point vector
std::vector<double> mesh;
// mesh points (device buffer)
sycl::buffer<double, 1> mesh_d;
// mesh point vector staggered at half points
std::vector<double> mesh_staggered;
// Fourier wavenumbers corresponding to mesh
Expand All @@ -47,15 +43,11 @@ class Mesh {
std::vector<double> poisson_factor;
// Factor to use in combined field solve and E = -Grad(phi)
std::vector<Complex> poisson_E_factor;
sycl::buffer<Complex, 1> poisson_E_factor_d;

// charge density
std::vector<double> charge_density;
sycl::buffer<double, 1> charge_density_d;
// electric field
std::vector<double> electric_field;
// electric field (device)
sycl::buffer<double, 1> electric_field_d;
// electric field on a staggered grid
std::vector<double> electric_field_staggered;
// electrostatic potential
Expand All @@ -64,33 +56,6 @@ class Mesh {
// Calculate a particle's contribution to the electric field
double evaluate_electric_field(const double x);

/*
* Evaluate the electric field at x grid points by
* interpolating onto the grid
* SYCL note: this is a copy of evaluate_electric_field, but able to be called
* in sycl. This should become evaluate_electric_field eventually
*/
inline double
sycl_evaluate_electric_field(sycl::accessor<double> mesh_d,
sycl::accessor<double> electric_field_d,
double x) {

// Find grid cell that x is in
int index = sycl_get_left_index(x, mesh_d);

// now x is in the cell ( mesh[index-1], mesh[index] )

double cell_width = mesh_d[index + 1] - mesh_d[index];
double distance_into_cell = x - mesh_d[index];

// r is the proportion if the distance into the cell that the particle is at
// e.g. midpoint => r = 0.5
double r = distance_into_cell / cell_width;

return (1.0 - r) * electric_field_d[index] +
r * electric_field_d[index + 1];
};

// Deposit particle onto mesh
void deposit(Plasma &plasma);
void sycl_deposit(sycl::queue &Q, Plasma &plasma);
Expand Down Expand Up @@ -118,16 +83,46 @@ class Mesh {
// Given a point x and a grid, find the indices of the grid points
// either side of x
int get_left_index(const double x, const std::vector<double> mesh);
};

inline int sycl_get_left_index(const double x,
const sycl::accessor<double> mesh_d) {
namespace Mesh1D {

int index = 0;
while (mesh_d[index + 1] < x and index < int(mesh.size())) {
index++;
};
return index;
}
};
template <typename T>
inline int sycl_get_left_index(const double x, const T &mesh_d,
const int mesh_size) {

int index = 0;
while ((mesh_d[index + 1] < x) and (index < mesh_size)) {
index++;
};
return index;
}

/*
* Evaluate the electric field at x grid points by
* interpolating onto the grid
* SYCL note: this is a copy of evaluate_electric_field, but able to be called
* in sycl. This should become evaluate_electric_field eventually
*/
inline double sycl_evaluate_electric_field(
const sycl::accessor<double> &mesh_d, const int mesh_size,
const sycl::accessor<double> &electric_field_d, double x) {

// Find grid cell that x is in
int index = sycl_get_left_index(x, mesh_d, mesh_size);

// now x is in the cell ( mesh[index-1], mesh[index] )

double cell_width = mesh_d[index + 1] - mesh_d[index];
double distance_into_cell = x - mesh_d[index];

// r is the proportion if the distance into the cell that the particle is at
// e.g. midpoint => r = 0.5
double r = distance_into_cell / cell_width;

return (1.0 - r) * electric_field_d[index] + r * electric_field_d[index + 1];
}

} // namespace Mesh1D

#endif // __MESH_H__
13 changes: 0 additions & 13 deletions include/species.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,19 +30,10 @@ class Species {
double vth;
// particle position array
std::vector<double> x;
// particle positions (device)
sycl::buffer<double, 1> x_d;
// particle velocity structure of arrays
Velocity v;
// particle x velocity (device)
sycl::buffer<double, 1> vx_d;
// particle y velocity (device)
sycl::buffer<double, 1> vy_d;
// particle z velocity (device)
sycl::buffer<double, 1> vz_d;
// charge density of species (if adiabatic)
double charge_density;
sycl::buffer<double, 1> charge_density_d;
// particle position array at
// next timestep
std::vector<double> xnew;
Expand All @@ -51,8 +42,6 @@ class Species {
std::vector<double> vnew;
// particle weight
std::vector<double> w;
// particle weights (device)
sycl::buffer<double, 1> w_d;
// particle pusher
void push(sycl::queue &q, Mesh *mesh);
// set array sizes for particle properties
Expand All @@ -62,8 +51,6 @@ class Species {
// Coefficients for particle pusher
double dx_coef;
double dv_coef;
sycl::buffer<double, 1> dx_coef_d;
sycl::buffer<double, 1> dv_coef_d;
};

#endif // __SPECIES_H__
28 changes: 19 additions & 9 deletions src/diagnostics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,9 +41,11 @@ void Diagnostics::compute_field_energy(sycl::queue &Q, Mesh &mesh) {
// Create buffer using host allocated "data" array
sycl::buffer<double, 1> buf(data.data(), sycl::range<1>{data.size()});

sycl::buffer<double, 1> electric_field_d(
mesh.electric_field.data(), sycl::range<1>{mesh.electric_field.size()});
Q.submit([&](sycl::handler &h) {
auto writeresult = sycl::accessor(buf, h);
auto electric_field_a = sycl::accessor(mesh.electric_field_d, h);
auto electric_field_a = sycl::accessor(electric_field_d, h);
h.parallel_for(sycl::range<1>{size_t(num_steps)}, [=](auto idx) {
writeresult[idx[0]] = std::pow(electric_field_a[idx], 2);
});
Expand Down Expand Up @@ -80,17 +82,25 @@ void Diagnostics::compute_particle_energy(sycl::queue &Q, Plasma &plasma) {

// Create buffer using host allocated "data" array
sycl::buffer<double, 1> buf(data, sycl::range<1>{size_t(n)});
sycl::buffer<double, 1> vx_d(
plasma.kinetic_species.at(j).v.x.data(),
sycl::range<1>{plasma.kinetic_species.at(j).v.x.size()});
sycl::buffer<double, 1> vy_d(
plasma.kinetic_species.at(j).v.y.data(),
sycl::range<1>{plasma.kinetic_species.at(j).v.y.size()});
sycl::buffer<double, 1> vz_d(
plasma.kinetic_species.at(j).v.z.data(),
sycl::range<1>{plasma.kinetic_species.at(j).v.z.size()});
sycl::buffer<double, 1> w_d(
plasma.kinetic_species.at(j).w.data(),
sycl::range<1>{plasma.kinetic_species.at(j).w.size()});

Q.submit([&](sycl::handler &h) {
sycl::accessor species_energy(buf, h, sycl::write_only);
auto vx_a = plasma.kinetic_species.at(j)
.vx_d.get_access<sycl::access::mode::read_write>(h);
auto vy_a = plasma.kinetic_species.at(j)
.vy_d.get_access<sycl::access::mode::read_write>(h);
auto vz_a = plasma.kinetic_species.at(j)
.vz_d.get_access<sycl::access::mode::read_write>(h);
auto w_a = plasma.kinetic_species.at(j)
.w_d.get_access<sycl::access::mode::read_write>(h);
auto vx_a = vx_d.get_access<sycl::access::mode::read_write>(h);
auto vy_a = vy_d.get_access<sycl::access::mode::read_write>(h);
auto vz_a = vz_d.get_access<sycl::access::mode::read_write>(h);
auto w_a = w_d.get_access<sycl::access::mode::read_write>(h);

h.parallel_for(sycl::range<1>{size_t(n)}, [=](sycl::id<1> idx) {
species_energy[idx[0]] =
Expand Down
Loading

0 comments on commit b78aac7

Please sign in to comment.