Cuda: 11.2 GPU: NVIDIA T4 OS: Ubuntu 16.04 C++14
When I tried to run my project, this error occurred.
In real code, I did very simply funtion but it failed:
terminate called after throwing an instance of 'thrust::system::system_error'
what(): for_each: failed to synchronize: cudaErrorIllegalAddress: an illegal memory access was encountered
Aborted (core dumped)
My project only contains four major files (main.cu, Grid.cu, Grid.h, Particle.h):
Compile command:
nvcc -o out -g -std=c++14 -O2 --expt-extended-lambda --expt-relaxed-constexpr -DTHRUST_DEBUG -arch=sm_70 -Wno-deprecated-declarations -Xcudafe --diag_suppress=esa_on_defaulted_function_ignored -rdc=true main.cu Grid.cu Matrix2D.cu Particle.cu Vector2D.cu
main.cu:
#include <iostream>
#include <vector>
#include "Constants.h"
#include "Particle.h"
#include "Grid.h"
using namespace std;
Grid *grid;
__host__ void init() {
vector<Particle> p;
for (int i = 0; i < 10000; ++i)
p.push_back(Particle());
grid = new Grid(Vector2D(0, 0), Vector2D(WIN_METERS_X, WIN_METERS_Y), Vector2D(256, 128), &p);
grid->initGridMassVel();
}
__host__ void update() {
}
int main()
{
init();
return 0;
}
Grid.h:
#ifndef GRID_H
#define GRID_H
#include "Constants.h"
#include "Vector2D.h"
#include "Particle.h"
#include <cmath>
#include <cstring>
#include <vector>
#include <thrust/functional.h>
#include <thrust/execution_policy.h>
#include <thrust/for_each.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <thrust/copy.h>
#include <thrust/device_vector.h>
using namespace std;
struct Node {
double mass;
int active;
Vector2D pos, vel, vel_new, force;
__host__ __device__ Node() {
mass = 0;
active = 0;
pos = vel = vel_new = force = Vector2D();
}
};
class Grid
{
public:
// Grid origin and size
Vector2D origin, size;
// Grid nodes
Vector2D node_size;
double node_area;
int nodes_length;
__host__ Grid(Vector2D pos, Vector2D dims, Vector2D window_size, vector<Particle>* _particles) :
origin(pos), size(window_size) {
node_size = dims / window_size;
size.x += 1, size.y += 1;
nodes_length = int(size.x * size.y);
node_area = node_size.x * node_size.y;
particles.resize(_particles->size());
thrust::copy(_particles->begin(), _particles->end(), particles.begin());
nodes.resize(nodes_length);
vector<Node> temp_nodes;
temp_nodes.resize(nodes_length);
for (int y = 0; y < size.y; ++y) {
for (int x = 0; x < size.x; ++x) {
int node_id = int(y * size.x + x);
temp_nodes[node_id].pos = Vector2D(x, y);
}
}
thrust::copy(temp_nodes.begin(), temp_nodes.end(), nodes.begin());
}
// Map particles to Grid
__host__ void initGridMassVel();
private:
thrust::device_vector<Particle> particles;
thrust::device_vector<Node> nodes;
};
#endif
Grid.cu:
#include "Grid.h"
#include <cassert>
#include <cstdio>
#define MATRIX_EPSILON 1e-6
__host__ __device__ double bspline(double x) {
x = fabs(x);
double w;
if (x < 1)
w = x * x * (x / 2 - 1) + 2 / 3.0;
else if (x < 2)
w = x * (x * (-x / 6 + 1) - 2) + 4 / 3.0;
else return 0;
return w;
}
//Slope of interpolation function
__host__ __device__ double bsplineSlope(double x) {
double abs_x = fabs(x);
if (abs_x < 1)
return 1.5 * x * abs_x - 2 * x;
else if (x < 2)
return -x * abs_x / 2 + 2 * x - 2 * x / abs_x;
else return 0;
}
__device__ double my_atomicAdd(double* address, double val)
{
unsigned long long int* address_as_ull =
(unsigned long long int*)address;
unsigned long long int old = *address_as_ull, assumed;
do {
assumed = old;
old = atomicCAS(address_as_ull, assumed,
__double_as_longlong(val +
__longlong_as_double(assumed)));
// Note: uses integer comparison to avoid hang in case of NaN (since NaN != NaN)
} while (assumed != old);
return __longlong_as_double(old);
}
__host__ void Grid::initGridMassVel() {
// Map particle to grid
Node* grid_ptr = thrust::raw_pointer_cast(&nodes[0]);
auto func = [=] __device__ (Particle & p) {
// get the index of the grid cross point corresponding to the particle (it is on the bottom left of the particle)
p.grid_p = (p.pos - origin) / node_size;
int p_x = (int)p.grid_p.x; // x coord index in grid
int p_y = (int)p.grid_p.y; // y coord index in grid
//printf("P: %d %d\n", p_x, p_y);
// Map from (p_x - 1, p_y - 1) to (p_x + 2, p_y + 2)
// The origin is bottom left, which means node_id = y * size.x + x
for (int it = 0, y = p_y - 1; y <= p_y + 2; ++y) {
if (y < 0 || y >= size.y) // here size.y has already been added by 1
continue;
// Y interpolation
double weight_y = bspline(p.grid_p.y - y);
double dy = bsplineSlope(p.grid_p.y - y);
for (int x = p_x - 1; x <= p_x + 2; ++x, ++it) {
if (x < 0 || x >= size.x)
continue;
// X interpolation
double weight_x = bspline(p.grid_p.x - x);
double dx = bsplineSlope(p.grid_p.x - x);
// set weight of particles related nodes
double w = weight_x * weight_y;
p.weights[it] = w;
// set weight gradient
p.weight_gradient[it] = Vector2D(dx * weight_y, dy * weight_x);
p.weight_gradient[it] /= node_size;
// set node weighted mass and velocity
int node_id = int(y * size.x + x);
//nodes[node_id].mass += w * p.mass;
my_atomicAdd(&(grid_ptr[node_id].mass), w * p.mass);
//nodes[node_id].vel += p.vel * w * p.mass;
Vector2D temp = p.vel * w * p.mass;
my_atomicAdd(&(grid_ptr[node_id].vel.x), temp.x);
my_atomicAdd(&(grid_ptr[node_id].vel.y), temp.y);
//nodes[node_id].active = true;
atomicAdd(&(grid_ptr[node_id].active), 1);
}
}
};
thrust::for_each(thrust::device, particles.begin(), particles.end(), func);
thrust::for_each(
thrust::device,
nodes.begin(),
nodes.end(),
[=] __device__(Node& n) {
if (n.active)
n.vel /= n.mass;
}
);
}
Particle.h:
#include <cstring>
#include "Vector2D.h"
#include "Matrix2D.h"
#include <cuda.h>
#include <cuda_runtime.h>
using namespace std;
class Particle
{
public:
double mass, density, volume;
Vector2D pos, vel;
Matrix2D velocity_gradient;
Matrix2D deformation_gradient, plastic_deformation, elastic_deformation;
Matrix2D stress;
// assigned grid index
Vector2D grid_p;
// weight value for nearest 16 grid nodes
double weights[16];
Vector2D weight_gradient[16];
// TODO: need more variables here like deformation;
__host__ __device__ Particle() {};
__host__ __device__ Particle(const Vector2D& pos, const Vector2D& vel, double mass) :
pos(pos), vel(vel), mass(mass) {
// TODO: more varibles means re-write init function
this->deformation_gradient = Matrix2D(); // init as identity matrix
this->plastic_deformation = Matrix2D(); // init as identity matrix
this->elastic_deformation = Matrix2D(); // init as identity matrix
this->volume = 1; // init volume is 1 (a unit volume)
memset(weights, 0, sizeof(double) * 16);
memset(weight_gradient, 0, sizeof(Vector2D) * 16);
}
};