Skip to content

Commit

Permalink
grid: Avoid floating point exception when cell diagonal is zero
Browse files Browse the repository at this point in the history
  • Loading branch information
oschuett committed Jun 4, 2021
1 parent 50ea084 commit a6cabf1
Showing 1 changed file with 20 additions and 19 deletions.
39 changes: 20 additions & 19 deletions src/grid/gpu/grid_gpu_task_list.cu
Original file line number Diff line number Diff line change
Expand Up @@ -156,25 +156,26 @@ create_tasks(const bool orthorhombic, const int ntasks,
task->ab_block_offset = block_offset + subblock_offset;

// Stuff for the ortho kernel ----------------------------------------------

// Discretize the radius.
const double drmin =
fmin(dh[level][0][0], fmin(dh[level][1][1], dh[level][2][2]));
const int imr = imax(1, (int)ceil(task->radius / drmin));
task->disr_radius = drmin * imr;

// Position of the gaussian product:
// this is the actual definition of the position on the grid
// i.e. a point rp(:) gets here grid coordinates
// MODULO(rp(:)/dr(:),npts_global(:))+1
// hence (0.0,0.0,0.0) in real space is rsgrid%lb on the rsgrid in Fortran
// and (1,1,1) on grid here in C.
//
// cubecenter(:) = FLOOR(MATMUL(dh_inv, rp))
for (int i = 0; i < 3; i++) {
const int cubecenter = floor(task->gp[i]);
task->cube_center_shifted[i] = cubecenter - shift_local[level][i];
task->cube_offset[i] = cubecenter * dh[level][i][i] - task->rp[i];
if (orthorhombic) {
// Discretize the radius.
const double drmin =
fmin(dh[level][0][0], fmin(dh[level][1][1], dh[level][2][2]));
const int imr = imax(1, (int)ceil(task->radius / drmin));
task->disr_radius = drmin * imr;

// Position of the gaussian product:
// this is the actual definition of the position on the grid
// i.e. a point rp(:) gets here grid coordinates
// MODULO(rp(:)/dr(:),npts_global(:))+1
// hence (0.0,0.0,0.0) in real space is rsgrid%lb on the rsgrid in Fortran
// and (1,1,1) on grid here in C.
//
// cubecenter(:) = FLOOR(MATMUL(dh_inv, rp))
for (int i = 0; i < 3; i++) {
const int cubecenter = floor(task->gp[i]);
task->cube_center_shifted[i] = cubecenter - shift_local[level][i];
task->cube_offset[i] = cubecenter * dh[level][i][i] - task->rp[i];
}
}

// Stuff for the general kernel --------------------------------------------
Expand Down

0 comments on commit a6cabf1

Please sign in to comment.