Skip to content

Commit

Permalink
Rename StaggeredGridDirection -> StaggerGrid for conciseness.
Browse files Browse the repository at this point in the history
  • Loading branch information
@tranqui
tranqui committed Jul 8, 2024
1 parent 4481ef5 commit 162531c
Show file tree
Hide file tree
Showing 3 changed files with 87 additions and 29 deletions.
29 changes: 14 additions & 15 deletions src/finite_difference.cuh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,28 +122,27 @@ namespace kernel
* We cannot use half-integer indices internally for data representation, so we have
* to use an implicit offset from integral indices.
*
* @tparam StaggerDirection: direction of staggering of grid from integer indices.
* If StaggerDirection=Left : grad[i] -> grad(i - 1/2)
* If StaggerDirection=Right: grad[i] -> grad(i + 1/2)
* Generally: grad[i] -> grad(i + StaggerDirection/2)
* @tparam Stagger: direction of staggering of grid from integer indices.
* If Stagger=Left : grad[i] -> grad(i - 1/2)
* If Stagger=Right: grad[i] -> grad(i + 1/2)
* @tparam Order: Derivatives must be implemented for the specific
* Order of the numerical approximation, which is done below
* by complete class specialisation.
*/
template <StaggeredGridDirection StaggerDirection, int Order>
template <StaggerGrid Stagger, int Order>
struct StaggeredDerivative
: public BaseDerivative<StaggeredDerivative<StaggerDirection, Order>>
: public BaseDerivative<StaggeredDerivative<Stagger, Order>>
{ };

/// Second-order staggered finite-difference derivatives.
template <StaggeredGridDirection StaggerDirection>
struct StaggeredDerivative<StaggerDirection, 2>
: public BaseDerivative<StaggeredDerivative<StaggerDirection, 2>>
template <StaggerGrid Stagger>
struct StaggeredDerivative<Stagger, 2>
: public BaseDerivative<StaggeredDerivative<Stagger, 2>>
{
template <typename T>
static __device__ inline Scalar grad_x(T&& tile, int i, int j)
{
if constexpr (StaggerDirection == Right)
if constexpr (Stagger == Right)
return stencil.dxInv * (tile[i][j+1] - tile[i][j]);
else
return stencil.dxInv * (tile[i][j] - tile[i][j-1]);
Expand All @@ -152,7 +151,7 @@ namespace kernel
template <typename T>
static __device__ inline Scalar grad_y(T&& tile, int i, int j)
{
if constexpr (StaggerDirection == Right)
if constexpr (Stagger == Right)
return stencil.dyInv * (tile[i+1][j] - tile[i][j]);
else
return stencil.dyInv * (tile[i][j] - tile[i-1][j]);
Expand All @@ -161,7 +160,7 @@ namespace kernel
template <typename T>
static __device__ inline Scalar lap_x(T&& tile, int i, int j)
{
if constexpr (StaggerDirection == Right)
if constexpr (Stagger == Right)
return 0.25 * stencil.dxInv * stencil.dxInv * (
(tile[i][j+2] - tile[i][ j ] + tile[i+1][j+2] - tile[i+1][ j ])
- (tile[i][j+1] - tile[i][j-1] + tile[i+1][j+1] - tile[i+1][j-1]) );
Expand All @@ -174,7 +173,7 @@ namespace kernel
template <typename T>
static __device__ inline Scalar lap_y(T&& tile, int i, int j)
{
if constexpr (StaggerDirection == Right)
if constexpr (Stagger == Right)
return 0.25 * stencil.dyInv * stencil.dyInv * (
(tile[i+2][j] - tile[ i ][j] + tile[i+2][j+1] - tile[ i ][j+1])
- (tile[i+1][j] - tile[i-1][j] + tile[i+1][j+1] - tile[i-1][j+1]) );
Expand All @@ -187,6 +186,6 @@ namespace kernel
}

using CentralDifference = details::CentralDerivative<order>;
template <StaggeredGridDirection StaggerDirection>
using StaggeredDifference = details::StaggeredDerivative<StaggerDirection, order>;
template <StaggerGrid Stagger>
using StaggeredDifference = details::StaggeredDerivative<Stagger, order>;
}
2 changes: 1 addition & 1 deletion src/math_primitives.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,4 +21,4 @@ using DeviceCurrent = DeviceGradient;

// Direction of index offset when derivatives are taken on staggered grids.
// See `StaggeredDerivative` in finite_differences.cuh for more info.
enum StaggeredGridDirection { Left, Central, Right };
enum StaggerGrid { Left, Central, Right };
85 changes: 72 additions & 13 deletions test/unit_tests.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,8 @@
#include "foreach.h"
#include "finite_difference.h"

using namespace finite_difference;


// Numerical tolerance for equality with numerical calculations.
static constexpr Scalar tight_tol = 1e-12;
Expand Down Expand Up @@ -151,8 +153,8 @@ inline Gradient gradient(Field field, Stencil stencil)

auto grad = [&](int i, int j, auto stencil_y, auto stencil_x)
{
gradient[0](i, j) = finite_difference::first(stencil_y) / stencil.dy;
gradient[1](i, j) = finite_difference::first(stencil_x) / stencil.dx;
gradient[0](i, j) = first(stencil_y) / stencil.dy;
gradient[1](i, j) = first(stencil_x) / stencil.dx;
};
apply_operation(field, grad);

Expand All @@ -168,15 +170,37 @@ inline Field laplacian(const FieldRef& field, Stencil stencil)

auto lap = [&](int i, int j, auto stencil_y, auto stencil_x)
{
laplacian(i, j) = dyInv*dyInv * finite_difference::second(stencil_y)
+ dxInv*dxInv * finite_difference::second(stencil_x);
laplacian(i, j) = dyInv*dyInv * second(stencil_y)
+ dxInv*dxInv * second(stencil_x);
};
apply_operation(field, lap);

return laplacian;
}

template <StaggeredGridDirection StaggerDirection>
// Find divergence of vector field by central finite differences.
inline Field divergence(const Gradient& grad, Stencil stencil)
{
const Scalar dxInv{1/stencil.dx}, dyInv{1/stencil.dy};
const auto Ny{grad[0].rows()}, Nx{grad[0].cols()};
Field divergence = Field::Zero(Ny, Nx);

auto div_y = [&](int i, int j, auto stencil_y, auto stencil_x)
{
divergence(i, j) += dyInv * first(stencil_y);
};
auto div_x = [&](int i, int j, auto stencil_y, auto stencil_x)
{
divergence(i, j) += dxInv * first(stencil_x);
};

apply_operation(grad[0], div_y);
apply_operation(grad[1], div_x);

return divergence;
}

template <StaggerGrid Stagger>
inline Gradient staggered_gradient(Field field, Stencil stencil)
{
const auto Ny{field.rows()}, Nx{field.cols()};
Expand All @@ -186,7 +210,7 @@ inline Gradient staggered_gradient(Field field, Stencil stencil)
{
// Nearest neighbours in y-direction w/ periodic boundaries:
int ip{i};
if constexpr (StaggerDirection == Right) ip++;
if constexpr (Stagger == Right) ip++;
int im{ip-1};
if (im < 0) im += Ny;
if (ip >= Ny) ip -= Ny;
Expand All @@ -195,7 +219,7 @@ inline Gradient staggered_gradient(Field field, Stencil stencil)
{
// Nearest neighbours in x-direction w/ periodic boundaries:
int jp{j};
if constexpr (StaggerDirection == Right) jp++;
if constexpr (Stagger == Right) jp++;
int jm{jp-1};
if (jm < 0) jm += Nx;
if (jp >= Nx) jp -= Nx;
Expand All @@ -208,7 +232,7 @@ inline Gradient staggered_gradient(Field field, Stencil stencil)
return grad;
}

template <StaggeredGridDirection StaggerDirection>
template <StaggerGrid Stagger>
inline Field staggered_laplacian(Field field, Stencil stencil)
{
const auto Ny{field.rows()}, Nx{field.cols()};
Expand All @@ -218,7 +242,7 @@ inline Field staggered_laplacian(Field field, Stencil stencil)
{
// Nearest neighbours in y-direction w/ periodic boundaries:
int i3{i}, i4{i+1};
if constexpr (StaggerDirection == Right)
if constexpr (Stagger == Right)
{
i3++;
i4++;
Expand All @@ -233,7 +257,7 @@ inline Field staggered_laplacian(Field field, Stencil stencil)
{
// Nearest neighbours in x-direction w/ periodic boundaries:
int j3{j}, j4{j+1};
if constexpr (StaggerDirection == Right)
if constexpr (Stagger == Right)
{
j3++;
j4++;
Expand Down Expand Up @@ -263,7 +287,7 @@ inline Field staggered_laplacian(Field field, Stencil stencil)
return lap;
}

template <StaggeredGridDirection StaggerDirection>
template <StaggerGrid Stagger>
inline Field staggered_divergence(Gradient grad, Stencil stencil)
{
const auto Ny{grad[0].rows()}, Nx{grad[0].cols()};
Expand All @@ -273,7 +297,7 @@ inline Field staggered_divergence(Gradient grad, Stencil stencil)
{
// Nearest neighbours in y-direction w/ periodic boundaries:
int ip{i};
if constexpr (StaggerDirection == Right) ip++;
if constexpr (Stagger == Right) ip++;
int im{ip-1};
if (im < 0) im += Ny;
if (ip >= Ny) ip -= Ny;
Expand All @@ -282,7 +306,7 @@ inline Field staggered_divergence(Gradient grad, Stencil stencil)
{
// Nearest neighbours in x-direction w/ periodic boundaries:
int jp{j};
if constexpr (StaggerDirection == Right) jp++;
if constexpr (Stagger == Right) jp++;
int jm{jp-1};
if (jm < 0) jm += Nx;
if (jp >= Nx) jp -= Nx;
Expand All @@ -306,10 +330,45 @@ TEST_CASE("FiniteDifferencesTest")
int Nx{64}, Ny{32};
Stencil stencil{1e-2, 1, 0.75};

// Check $\nabla \cdot (\nabla \phi) = \nabla^2 \phi.$
{
Field field = Field::Random(Ny, Nx);
Gradient grad = staggered_gradient<Right>(field, stencil);
Field lap = laplacian(field, stencil);
CHECK(is_equal<tight_tol>(lap, staggered_divergence<Left>(grad, stencil)));
}

// Check derivatives of known analytic functions

// $\phi = x$:
{
Field field(Ny, Nx);
for (int i = 0; i < Ny; ++i)
for (int j = 0; j < Nx; ++j)
field(i, j) = j;

Gradient grad = gradient(field, stencil);
Field lap = laplacian(field, stencil);

CHECK(is_equal<tight_tol>(lap(Ny/2, Nx/2), 0));
CHECK(is_equal<tight_tol>(grad[0](Ny/2, Nx/2), 0));
CHECK(is_equal<tight_tol>(grad[1](Ny/2, Nx/2), 1));
}

// $\phi = (x^2) / 2$:
{
Field field(Ny, Nx);
for (int i = 0; i < Ny; ++i)
for (int j = 0; j < Nx; ++j)
field(i, j) = 0.5 * j*j;

Gradient grad = gradient(field, stencil);
Field lap = laplacian(field, stencil);

CHECK(is_equal<tight_tol>(lap(Ny/2, Nx/2), 1));
CHECK(is_equal<tight_tol>(grad[0](Ny/2, Nx/2), 0));
CHECK(is_equal<tight_tol>(grad[1](Ny/2, Nx/2), Nx/2));
}
}

TEST_CASE("Constructor")
Expand Down

0 comments on commit 162531c

Please sign in to comment.