Skip to content

Commit

Permalink
api: fix staggered FD corner case and enforce satggered case to be va…
Browse files Browse the repository at this point in the history 
…lid staggering
  • Loading branch information
@mloubout
mloubout committed Mar 8, 2024
1 parent ae85ddf commit e13a83a
Show file tree
Hide file tree
Showing 5 changed files with 68 additions and 59 deletions.
9 changes: 4 additions & 5 deletions devito/finite_differences/finite_difference.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@ def cross_derivative(expr, dims, fd_order, deriv_order, x0=None, **kwargs):
finite difference. Defaults to `direct`.
x0 : dict, optional
Origin of the finite-difference scheme as a map dim: origin_dim.
coefficients : strong, optional
coefficients : string, optional
Use taylor or custom coefficients (weights). Defaults to taylor.
expand : bool, optional
If True, the derivative is fully expanded as a sum of products,
Expand Down Expand Up @@ -189,7 +189,7 @@ def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, x0=None,
finite difference. Defaults to `direct`.
x0 : dict, optional
Origin of the finite-difference scheme as a map dim: origin_dim.
coefficients : strong, optional
coefficients : string, optional
Use taylor or custom coefficients (weights). Defaults to taylor.
expand : bool, optional
If True, the derivative is fully expanded as a sum of products,
Expand All @@ -201,12 +201,12 @@ def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, x0=None,
``deriv-order`` derivative of ``expr``.
"""
side = None
# First order derivative with 2nd order FD is highly non-recommended so taking
# First order derivative with 2nd order FD is strongly discouraged so taking
# first order fd that is a lot better
if deriv_order == 1 and fd_order == 2 and coefficients != 'symbolic':
fd_order = 1

# Enforce sable time coefficients
# Enforce stable time coefficients
if dim.is_Time and coefficients != 'symbolic':
coefficients = 'taylor'

Expand All @@ -219,7 +219,6 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici
# The stencil indices
indices, x0 = generate_indices(expr, dim, fd_order, side=side, matvec=matvec,
x0=x0)

# Finite difference weights corresponding to the indices. Computed via the
# `coefficients` method (`taylor` or `symbolic`)
weights = fd_weights_registry[coefficients](expr, deriv_order, indices, x0)
Expand Down
4 changes: 2 additions & 2 deletions devito/finite_differences/rsfd.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ def drot(expr, dim, dir=1, x0=None):
Rotated finite differences based on:
https://www.sciencedirect.com/science/article/pii/S0165212599000232
The rotated axis (the four diagonals of a cube) are:
The rotated axes (the four diagonals of a cube) are:
d1 = dx/dr x + dz/dl y + dy/dl z
d2 = dx/dl x + dz/dl y - dy/dr z
d3 = dx/dr x - dz/dl y + dy/dl z
Expand Down Expand Up @@ -131,7 +131,7 @@ def d45(expr, dim, x0=None, expand=True):
Expand the expression. Defaults to True.
"""
# Make sure the grid supports RSFD
if expr.grid.dim == 1:
if expr.grid.dim not in [2, 3]:
raise ValueError('RSFD only supported in 2D and 3D')

# Diagonals weights
Expand Down
14 changes: 10 additions & 4 deletions devito/finite_differences/tools.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -363,7 +363,7 @@ def generate_indices_staggered(expr, dim, order, side=None, x0=None):
indices = [start - diff/2, start + diff/2]
indices = IndexSet(dim, indices)
else:
o_min = -order//2+1
o_min = -order//2 + 1
o_max = order//2

d = make_stencil_dimension(expr, o_min, o_max)
Expand All @@ -374,11 +374,17 @@ def generate_indices_staggered(expr, dim, order, side=None, x0=None):
indices = [start, start - diff]
indices = IndexSet(dim, indices)
else:
o_min = -order//2
o_max = order//2
if start == dim + diff/2:
start = dim
o_min = -order//2 + 1
o_max = order//2
else:
start = dim + diff/2
o_min = -order//2
o_max = order//2 - 1

d = make_stencil_dimension(expr, o_min, o_max)
iexpr = start + d * diff
iexpr = ind0 + d * diff
indices = IndexSet(dim, expr=iexpr)

return start, indices
Expand Down
66 changes: 31 additions & 35 deletions examples/seismic/tutorials/06_elastic_varying_parameters.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -392,22 +392,22 @@
"name": "stderr",
"output_type": "stream",
"text": [
"Operator `Kernel` ran in 0.23 s\n"
"Operator `Kernel` ran in 0.24 s\n"
]
},
{
"data": {
"text/plain": [
"PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
" PerfEntry(time=0.20520499999999994, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.2058980000000001, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section1', rank=None),\n",
" PerfEntry(time=0.005932999999999998, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.006354999999999981, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section2', rank=None),\n",
" PerfEntry(time=0.006593000000000015, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.007110000000000014, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section3', rank=None),\n",
" PerfEntry(time=0.005976000000000029, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.006563000000000019, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section4', rank=None),\n",
" PerfEntry(time=0.00602000000000003, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
" PerfEntry(time=0.006876000000000027, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
]
},
"execution_count": 14,
Expand Down Expand Up @@ -511,22 +511,22 @@
"name": "stderr",
"output_type": "stream",
"text": [
"Operator `Kernel` ran in 0.27 s\n"
"Operator `Kernel` ran in 0.31 s\n"
]
},
{
"data": {
"text/plain": [
"PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
" PerfEntry(time=0.23217499999999996, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.24031900000000017, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section1', rank=None),\n",
" PerfEntry(time=0.008015000000000006, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.015194000000000017, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section2', rank=None),\n",
" PerfEntry(time=0.012265999999999938, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.015348, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section3', rank=None),\n",
" PerfEntry(time=0.007797000000000036, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.015376999999999998, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section4', rank=None),\n",
" PerfEntry(time=0.009437999999999978, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
" PerfEntry(time=0.015624999999999997, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
]
},
"execution_count": 17,
Expand Down Expand Up @@ -699,20 +699,20 @@
"name": "stderr",
"output_type": "stream",
"text": [
"Operator `Kernel` ran in 0.24 s\n"
"Operator `Kernel` ran in 0.23 s\n"
]
},
{
"data": {
"text/plain": [
"PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
" PerfEntry(time=0.21262299999999978, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.20202699999999996, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section1', rank=None),\n",
" PerfEntry(time=0.008176999999999999, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.007293999999999994, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section2', rank=None),\n",
" PerfEntry(time=0.008674000000000015, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.007614000000000006, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section3', rank=None),\n",
" PerfEntry(time=0.008485000000000003, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
" PerfEntry(time=0.007164000000000016, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
]
},
"execution_count": 24,
Expand Down Expand Up @@ -796,20 +796,20 @@
"name": "stderr",
"output_type": "stream",
"text": [
"Operator `Kernel` ran in 0.25 s\n"
"Operator `Kernel` ran in 0.23 s\n"
]
},
{
"data": {
"text/plain": [
"PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
" PerfEntry(time=0.2165199999999999, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.19974000000000014, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section1', rank=None),\n",
" PerfEntry(time=0.00951100000000004, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.007187999999999981, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section2', rank=None),\n",
" PerfEntry(time=0.01492799999999995, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.007523000000000016, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section3', rank=None),\n",
" PerfEntry(time=0.00861800000000002, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
" PerfEntry(time=0.007113000000000016, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
]
},
"execution_count": 26,
Expand Down Expand Up @@ -946,16 +946,12 @@
{
"data": {
"text/latex": [
"$\\displaystyle \\left(- \\frac{f(x - h_x, y - h_y)}{4 h_{x}} + \\frac{f(x + h_x, y + h_y)}{4 h_{x}}\\right) + \\left(- \\frac{f(x - h_x, y + h_y)}{4 h_{x}} + \\frac{f(x + h_x, y - h_y)}{4 h_{x}}\\right)$"
"$\\displaystyle \\left(\\frac{f(x, y)}{2 h_{x}} - \\frac{f(x - h_x, y - h_y)}{2 h_{x}}\\right) + \\left(\\frac{f(x, y)}{2 h_{x}} - \\frac{f(x - h_x, y + h_y)}{2 h_{x}}\\right)$"
],
"text/plain": [
" f(x - hₓ, y - h_y) f(x + hₓ, y + h_y) f(x - hₓ, y + h_y) f(x + hₓ, y\n",
"- ────────────────── + ────────────────── + - ────────────────── + ───────────\n",
" 4⋅hₓ 4⋅hₓ 4⋅hₓ 4⋅hₓ\n",
"\n",
" - h_y)\n",
"───────\n",
" "
"f(x, y) f(x - hₓ, y - h_y) f(x, y) f(x - hₓ, y + h_y)\n",
"─────── - ────────────────── + ─────── - ──────────────────\n",
" 2⋅hₓ 2⋅hₓ 2⋅hₓ 2⋅hₓ "
]
},
"execution_count": 32,
Expand Down Expand Up @@ -1023,15 +1019,15 @@
"data": {
"text/plain": [
"PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
" PerfEntry(time=0.5155960000000004, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.5144760000000006, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section1', rank=None),\n",
" PerfEntry(time=0.011958000000000087, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.012137000000000096, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section2', rank=None),\n",
" PerfEntry(time=0.013548000000000008, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.013851000000000006, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section3', rank=None),\n",
" PerfEntry(time=0.012433999999999917, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.012857999999999913, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section4', rank=None),\n",
" PerfEntry(time=0.012537999999999919, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
" PerfEntry(time=0.012548999999999916, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
]
},
"execution_count": 34,
Expand Down
34 changes: 21 additions & 13 deletions tests/test_derivatives.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -265,14 +265,14 @@ def test_fd_space(self, derivative, space_order):
@pytest.mark.parametrize('staggered', [(True, True), (False, False),
(True, False), (False, True)])
@pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
def test_fd_space_45(self, staggered, space_order):
@pytest.mark.parametrize('ndim', [2, 3])
def test_fd_space_45(self, staggered, space_order, ndim):
"""
Rotated derivatives require at least 2D to get access to diagonal points.
We create a simple 1D gradient over a 2D grid to check against a polynomial
"""
# dummy axis dimension
nx = 100
ny = nx
xx = np.linspace(-1, 1, nx)
dx = xx[1] - xx[0]
if staggered[0] and not staggered[1]:
Expand All @@ -282,7 +282,7 @@ def test_fd_space_45(self, staggered, space_order):
else:
xx_s = xx
# Symbolic data
grid = Grid(shape=(nx, ny), dtype=np.float32)
grid = Grid(shape=tuple([nx]*ndim), dtype=np.float32)
x = grid.dimensions[0]
u = Function(name="u", grid=grid, space_order=space_order,
staggered=None if staggered[0] else grid.dimensions)
Expand All @@ -293,8 +293,7 @@ def test_fd_space_45(self, staggered, space_order):
polynome = sum([coeffs[i]*x**i for i in range(0, space_order)])
polyvalues = np.array([polynome.subs(x, xi) for xi in xx], np.float32)
# Fill original data with the polynomial values
for i in range(ny):
u.data[:, i] = polyvalues
u.data[:] = polyvalues.reshape(nx, *[1]*(ndim-1))
# True derivative of the polynome
Dpolynome = diff(polynome)
Dpolyvalues = np.array([Dpolynome.subs(x, xi) for xi in xx_s], np.float32)
Expand All @@ -307,13 +306,15 @@ def test_fd_space_45(self, staggered, space_order):

# Check exactness of the numerical derivative except inside space_brd
space_border = space_order
error = abs(du.data[space_border:-space_border, ny//2] -
mid = tuple([slice(space_border, -space_border, 1)] +
[nx//2 for i in range(ndim-1)])
error = abs(du.data[mid] -
Dpolyvalues[space_border:-space_border])

assert np.isclose(np.mean(error), 0., atol=1e-3)

@pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
@pytest.mark.parametrize('stagger', [centered, left, right])
@pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
# Only test x and t as y and z are the same as x
def test_fd_space_staggered(self, space_order, stagger):
"""
Expand All @@ -333,26 +334,32 @@ def test_fd_space_staggered(self, space_order, stagger):
if stagger == left:
off = -.5
side = -x
xx2 = xx + off * dx
xx_u = xx + off * dx
duside = NODE
xx_du = xx
elif stagger == right:
off = .5
side = x
xx2 = xx + off * dx
xx_u = xx + off * dx
duside = NODE
xx_du = xx
else:
side = NODE
xx2 = xx
xx_u = xx
xx_du = xx + .5 * dx
duside = x

u = Function(name="u", grid=grid, space_order=space_order, staggered=side)
du = Function(name="du", grid=grid, space_order=space_order, staggered=side)
du = Function(name="du", grid=grid, space_order=space_order, staggered=duside)
# Define polynomial with exact fd
coeffs = np.ones((space_order-1,), dtype=np.float32)
polynome = sum([coeffs[i]*x**i for i in range(0, space_order-1)])
polyvalues = np.array([polynome.subs(x, xi) for xi in xx2], np.float32)
polyvalues = np.array([polynome.subs(x, xi) for xi in xx_u], np.float32)
# Fill original data with the polynomial values
u.data[:] = polyvalues
# True derivative of the polynome
Dpolynome = diff(polynome)
Dpolyvalues = np.array([Dpolynome.subs(x, xi) for xi in xx2], np.float32)
Dpolyvalues = np.array([Dpolynome.subs(x, xi) for xi in xx_du], np.float32)
# Compute numerical FD
stencil = Eq(du, u.dx)
op = Operator(stencil, subs={x.spacing: dx})
Expand Down Expand Up @@ -788,6 +795,7 @@ def test_dx2(self):
assert isinstance(term0, EvalDerivative)

term1 = f.dx2._evaluate(expand=False)
print(type(term1))
assert isinstance(term1, DiffDerivative)
assert term1.depth == 1
term1 = term1.evaluate
Expand Down

0 comments on commit e13a83a

Please sign in to comment.