/
test_symbolic_coefficients.py
96 lines (74 loc) · 3.24 KB
/
test_symbolic_coefficients.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
import numpy as np
import pytest
from conftest import skipif
from devito import (Grid, Function, TimeFunction, Eq, Coefficient, Substitutions,
Dimension, solve, Operator)
from devito.finite_differences import Differentiable
from devito.tools import as_tuple
_PRECISION = 9
pytestmark = skipif(['yask', 'ops'])
class TestSC(object):
"""
Class for testing symbolic coefficients functionality
"""
@pytest.mark.parametrize('order', [1, 2, 6])
def test_default_rules(self, order):
"""
Test that the default replacement rules return the same
as standard FD.
"""
grid = Grid(shape=(20, 20))
u0 = TimeFunction(name='u', grid=grid, time_order=order, space_order=order)
u1 = TimeFunction(name='u', grid=grid, time_order=order, space_order=order,
coefficients='symbolic')
eq0 = Eq(-u0.dx+u0.dt)
eq1 = Eq(u1.dt-u1.dx)
assert(eq0.evalf(_PRECISION).__repr__() == eq1.evalf(_PRECISION).__repr__())
@pytest.mark.parametrize('expr, sorder, dorder, dim, weights, expected', [
('u.dx', 2, 1, 0, (-0.6, 0.1, 0.6),
'0.1*u(x, y) - 0.6*u(x - h_x, y) + 0.6*u(x + h_x, y)'),
('u.dy2', 3, 2, 1, (0.121, -0.223, 1.648, -2.904),
'1.648*u(x, y) + 0.121*u(x, y - 2*h_y) - 0.223*u(x, y - h_y) \
- 2.904*u(x, y + h_y)')])
def test_coefficients(self, expr, sorder, dorder, dim, weights, expected):
"""Test that custom coefficients return the expected result"""
grid = Grid(shape=(10, 10))
u = Function(name='u', grid=grid, space_order=sorder, coefficients='symbolic')
x = grid.dimensions
order = dorder
dim = x[dim]
weights = np.array(weights)
coeffs = Coefficient(order, u, dim, weights)
eq = Eq(eval(expr), coefficients=Substitutions(coeffs))
assert isinstance(eq.lhs, Differentiable)
assert expected == str(eq.evaluate.lhs)
def test_function_coefficients(self):
"""Test that custom function coefficients return the expected result"""
so = 2
grid = Grid(shape=(4, 4))
f0 = TimeFunction(name='f0', grid=grid, space_order=so, coefficients='symbolic')
f1 = TimeFunction(name='f1', grid=grid, space_order=so)
x, y = grid.dimensions
s = Dimension(name='s')
ncoeffs = so+1
wshape = list(grid.shape)
wshape.append(ncoeffs)
wshape = as_tuple(wshape)
wdims = list(grid.dimensions)
wdims.append(s)
wdims = as_tuple(wdims)
w = Function(name='w', dimensions=wdims, shape=wshape)
w.data[:, :, 0] = -0.5/grid.spacing[0]
w.data[:, :, 1] = 0.0
w.data[:, :, 2] = 0.5/grid.spacing[0]
f_x_coeffs = Coefficient(1, f0, x, w)
subs = Substitutions(f_x_coeffs)
eq0 = Eq(f0.dt + f0.dx, 1, coefficients=subs)
eq1 = Eq(f1.dt + f1.dx, 1)
stencil0 = solve(eq0.evaluate, f0.forward)
stencil1 = solve(eq1.evaluate, f1.forward)
op0 = Operator(Eq(f0.forward, stencil0))
op1 = Operator(Eq(f1.forward, stencil1))
op0(time_m=0, time_M=5, dt=1.0)
op1(time_m=0, time_M=5, dt=1.0)
assert np.all(np.isclose(f0.data[:] - f1.data[:], 0.0, atol=1e-5, rtol=0))