Skip to content
Snippets Groups Projects
Commit efadec2e authored by Michael Kuron's avatar Michael Kuron :mortar_board: Committed by Michael Kuron
Browse files

Minor improvements to FiniteDifferenceStaggeredStencilDerivation

parent 3d9e0c6d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
pystencils/fd/derivation.py
+ 5
− 4
View file @ efadec2e
...@@ -303,11 +303,12 @@ class FiniteDifferenceStaggeredStencilDerivation: ...@@ -303,11 +303,12 @@ class FiniteDifferenceStaggeredStencilDerivation:
zero_counts[zero_count].append(weights) zero_counts[zero_count].append(weights)
best = zero_counts[max(zero_counts.keys())] best = zero_counts[max(zero_counts.keys())]
if len(best) > 1: # if there are multiple, pick the one that contains a nonzero center weight if len(best) > 1: # if there are multiple, pick the one that contains a nonzero center weight
center = [tuple(p + pos) for p in points].index((0, 0, 0)) center = [tuple(p + pos) for p in points].index((0, 0, 0)[:dim])
best = [b for b in best if b[center] != 0] best = [b for b in best if b[center] != 0]
if len(best) > 1: if len(best) > 1: # if there are still multiple, they are equivalent, so we average
raise NotImplementedError("more than one suitable set of weights found, don't know how to proceed") weights = sp.Add(*[sp.Matrix(b) for b in best]) / len(best)
weights = best[0] else:
weights = best[0]
assert weights assert weights
points_tuple = tuple([tuple(p + pos) for p in points]) points_tuple = tuple([tuple(p + pos) for p in points])
......
pystencils_tests/test_fd_derivation.ipynb
+ 3
− 0
View file @ efadec2e
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from pystencils.session import * from pystencils.session import *
from pystencils.fd.derivation import * from pystencils.fd.derivation import *
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# 2D standard stencils # 2D standard stencils
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
stencil = [(-1, 0), (1, 0), (0, -1), (0, 1), (0, 0)] stencil = [(-1, 0), (1, 0), (0, -1), (0, 1), (0, 0)]
standard_2d_00 = FiniteDifferenceStencilDerivation((0,0), stencil) standard_2d_00 = FiniteDifferenceStencilDerivation((0,0), stencil)
f = ps.fields("f: [2D]") f = ps.fields("f: [2D]")
standard_2d_00_res = standard_2d_00.get_stencil() standard_2d_00_res = standard_2d_00.get_stencil()
res = standard_2d_00_res.apply(f.center) res = standard_2d_00_res.apply(f.center)
expected = f[-1, 0] - 2 * f[0, 0] + f[1, 0] expected = f[-1, 0] - 2 * f[0, 0] + f[1, 0]
assert res == expected assert res == expected
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
assert standard_2d_00_res.accuracy == 2 assert standard_2d_00_res.accuracy == 2
assert not standard_2d_00_res.is_isotropic assert not standard_2d_00_res.is_isotropic
standard_2d_00_res standard_2d_00_res
``` ```
%% Output %% Output
Finite difference stencil of accuracy 2, isotropic error: False Finite difference stencil of accuracy 2, isotropic error: False
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
standard_2d_00.get_stencil().as_matrix() standard_2d_00.get_stencil().as_matrix()
``` ```
%% Output %% Output
$$\left[\begin{matrix}0 & 0 & 0\\1 & -2 & 1\\0 & 0 & 0\end{matrix}\right]$$ $$\left[\begin{matrix}0 & 0 & 0\\1 & -2 & 1\\0 & 0 & 0\end{matrix}\right]$$
⎡0 0 0⎤ ⎡0 0 0⎤
⎢ ⎥ ⎢ ⎥
⎢1 -2 1⎥ ⎢1 -2 1⎥
⎢ ⎥ ⎢ ⎥
⎣0 0 0⎦ ⎣0 0 0⎦
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# 2D isotropic stencils # 2D isotropic stencils
## second x-derivative ## second x-derivative
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
stencil = [(i, j) for i in (-1, 0, 1) for j in (-1, 0, 1)] stencil = [(i, j) for i in (-1, 0, 1) for j in (-1, 0, 1)]
isotropic_2d_00 = FiniteDifferenceStencilDerivation((0,0), stencil) isotropic_2d_00 = FiniteDifferenceStencilDerivation((0,0), stencil)
isotropic_2d_00_res = isotropic_2d_00.get_stencil(isotropic=True) isotropic_2d_00_res = isotropic_2d_00.get_stencil(isotropic=True)
assert isotropic_2d_00_res.is_isotropic assert isotropic_2d_00_res.is_isotropic
assert isotropic_2d_00_res.accuracy == 2 assert isotropic_2d_00_res.accuracy == 2
isotropic_2d_00_res isotropic_2d_00_res
``` ```
%% Output %% Output
Finite difference stencil of accuracy 2, isotropic error: True Finite difference stencil of accuracy 2, isotropic error: True
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
isotropic_2d_00_res.as_matrix() isotropic_2d_00_res.as_matrix()
``` ```
%% Output %% Output
$$\left[\begin{matrix}\frac{1}{12} & - \frac{1}{6} & \frac{1}{12}\\\frac{5}{6} & - \frac{5}{3} & \frac{5}{6}\\\frac{1}{12} & - \frac{1}{6} & \frac{1}{12}\end{matrix}\right]$$ $$\left[\begin{matrix}\frac{1}{12} & - \frac{1}{6} & \frac{1}{12}\\\frac{5}{6} & - \frac{5}{3} & \frac{5}{6}\\\frac{1}{12} & - \frac{1}{6} & \frac{1}{12}\end{matrix}\right]$$
⎡1/12 -1/6 1/12⎤ ⎡1/12 -1/6 1/12⎤
⎢ ⎥ ⎢ ⎥
⎢5/6 -5/3 5/6 ⎥ ⎢5/6 -5/3 5/6 ⎥
⎢ ⎥ ⎢ ⎥
⎣1/12 -1/6 1/12⎦ ⎣1/12 -1/6 1/12⎦
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
plt.figure(figsize=(2,2)) plt.figure(figsize=(2,2))
isotropic_2d_00_res.visualize() isotropic_2d_00_res.visualize()
``` ```
%% Output %% Output
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
expected_result = sp.Matrix([[1, -2, 1], [10, -20, 10], [1, -2, 1]]) / 12 expected_result = sp.Matrix([[1, -2, 1], [10, -20, 10], [1, -2, 1]]) / 12
assert expected_result == isotropic_2d_00_res.as_matrix() assert expected_result == isotropic_2d_00_res.as_matrix()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
type(isotropic_2d_00_res.as_matrix()) type(isotropic_2d_00_res.as_matrix())
``` ```
%% Output %% Output
sympy.matrices.dense.MutableDenseMatrix sympy.matrices.dense.MutableDenseMatrix
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
type(expected_result) type(expected_result)
``` ```
%% Output %% Output
sympy.matrices.dense.MutableDenseMatrix sympy.matrices.dense.MutableDenseMatrix
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Isotropic laplacian ## Isotropic laplacian
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
isotropic_2d_11 = FiniteDifferenceStencilDerivation((1,1), stencil) isotropic_2d_11 = FiniteDifferenceStencilDerivation((1,1), stencil)
isotropic_2d_11_res = isotropic_2d_11.get_stencil(isotropic=True) isotropic_2d_11_res = isotropic_2d_11.get_stencil(isotropic=True)
iso_laplacian = isotropic_2d_00_res.as_matrix() + isotropic_2d_11_res.as_matrix() iso_laplacian = isotropic_2d_00_res.as_matrix() + isotropic_2d_11_res.as_matrix()
iso_laplacian iso_laplacian
``` ```
%% Output %% Output
$$\left[\begin{matrix}\frac{1}{6} & \frac{2}{3} & \frac{1}{6}\\\frac{2}{3} & - \frac{10}{3} & \frac{2}{3}\\\frac{1}{6} & \frac{2}{3} & \frac{1}{6}\end{matrix}\right]$$ $$\left[\begin{matrix}\frac{1}{6} & \frac{2}{3} & \frac{1}{6}\\\frac{2}{3} & - \frac{10}{3} & \frac{2}{3}\\\frac{1}{6} & \frac{2}{3} & \frac{1}{6}\end{matrix}\right]$$
⎡1/6 2/3 1/6⎤ ⎡1/6 2/3 1/6⎤
⎢ ⎥ ⎢ ⎥
⎢2/3 -10/3 2/3⎥ ⎢2/3 -10/3 2/3⎥
⎢ ⎥ ⎢ ⎥
⎣1/6 2/3 1/6⎦ ⎣1/6 2/3 1/6⎦
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
expected_result = sp.Matrix([[1, 4, 1], [4, -20, 4], [1, 4, 1]]) / 6 expected_result = sp.Matrix([[1, 4, 1], [4, -20, 4], [1, 4, 1]]) / 6
assert iso_laplacian == expected_result assert iso_laplacian == expected_result
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# stencils for staggered fields # stencils for staggered fields
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
half = sp.Rational(1, 2) half = sp.Rational(1, 2)
fd_points_ex = ( fd_points_ex = (
(half, 0), (half, 0),
(-half, 0), (-half, 0),
(half, 1), (half, 1),
(half, -1), (half, -1),
(-half, 1), (-half, 1),
(-half, -1) (-half, -1)
) )
assert set(fd_points_ex) == set(FiniteDifferenceStaggeredStencilDerivation("E", 2).stencil) assert set(fd_points_ex) == set(FiniteDifferenceStaggeredStencilDerivation("E", 2).stencil)
fd_points_ey = ( fd_points_ey = (
(0, half), (0, half),
(0, -half), (0, -half),
(-1,-half), (-1,-half),
(-1, half), (-1, half),
(1, -half), (1, -half),
(1, half) (1, half)
) )
assert set(fd_points_ey) == set(FiniteDifferenceStaggeredStencilDerivation("N",2).stencil) assert set(fd_points_ey) == set(FiniteDifferenceStaggeredStencilDerivation("N",2).stencil)
fd_points_c = ( fd_points_c = (
(half, half), (half, half),
(-half, -half), (-half, -half),
(half, -half), (half, -half),
(-half, half) (-half, half)
) )
assert set(fd_points_c) == set(FiniteDifferenceStaggeredStencilDerivation("NE",2).stencil) assert set(fd_points_c) == set(FiniteDifferenceStaggeredStencilDerivation("NE",2).stencil)
assert len(FiniteDifferenceStaggeredStencilDerivation("E",3).points) == 10 assert len(FiniteDifferenceStaggeredStencilDerivation("E",3).points) == 10
assert len(FiniteDifferenceStaggeredStencilDerivation("NE",3).points) == 12 assert len(FiniteDifferenceStaggeredStencilDerivation("NE",3).points) == 12
assert len(FiniteDifferenceStaggeredStencilDerivation("TNE",3).points) == 8 assert len(FiniteDifferenceStaggeredStencilDerivation("TNE",3).points) == 8
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
c = ps.fields("c: [2D]") c = ps.fields("c: [2D]")
c3 = ps.fields("c3: [3D]") c3 = ps.fields("c3: [3D]")
assert FiniteDifferenceStaggeredStencilDerivation("E", 2, (0,)).apply(c) == c[1, 0] - c[0, 0] assert FiniteDifferenceStaggeredStencilDerivation("E", 2, (0,)).apply(c) == c[1, 0] - c[0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("W", 2, (0,)).apply(c) == c[0, 0] - c[-1, 0] assert FiniteDifferenceStaggeredStencilDerivation("W", 2, (0,)).apply(c) == c[0, 0] - c[-1, 0]
assert FiniteDifferenceStaggeredStencilDerivation("N", 2, (1,)).apply(c) == c[0, 1] - c[0, 0] assert FiniteDifferenceStaggeredStencilDerivation("N", 2, (1,)).apply(c) == c[0, 1] - c[0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("S", 2, (1,)).apply(c) == c[0, 0] - c[0, -1] assert FiniteDifferenceStaggeredStencilDerivation("S", 2, (1,)).apply(c) == c[0, 0] - c[0, -1]
assert FiniteDifferenceStaggeredStencilDerivation("E", 3, (0,)).apply(c3) == c3[1, 0, 0] - c3[0, 0, 0] assert FiniteDifferenceStaggeredStencilDerivation("E", 3, (0,)).apply(c3) == c3[1, 0, 0] - c3[0, 0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("W", 3, (0,)).apply(c3) == c3[0, 0, 0] - c3[-1, 0, 0] assert FiniteDifferenceStaggeredStencilDerivation("W", 3, (0,)).apply(c3) == c3[0, 0, 0] - c3[-1, 0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("N", 3, (1,)).apply(c3) == c3[0, 1, 0] - c3[0, 0, 0] assert FiniteDifferenceStaggeredStencilDerivation("N", 3, (1,)).apply(c3) == c3[0, 1, 0] - c3[0, 0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("S", 3, (1,)).apply(c3) == c3[0, 0, 0] - c3[0, -1, 0] assert FiniteDifferenceStaggeredStencilDerivation("S", 3, (1,)).apply(c3) == c3[0, 0, 0] - c3[0, -1, 0]
assert FiniteDifferenceStaggeredStencilDerivation("T", 3, (2,)).apply(c3) == c3[0, 0, 1] - c3[0, 0, 0] assert FiniteDifferenceStaggeredStencilDerivation("T", 3, (2,)).apply(c3) == c3[0, 0, 1] - c3[0, 0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("B", 3, (2,)).apply(c3) == c3[0, 0, 0] - c3[0, 0, -1] assert FiniteDifferenceStaggeredStencilDerivation("B", 3, (2,)).apply(c3) == c3[0, 0, 0] - c3[0, 0, -1]
assert FiniteDifferenceStaggeredStencilDerivation("S", 2, (0,)).apply(c) == \
(c[1, 0] + c[1, -1] - c[-1, 0] - c[-1, -1])/4
assert FiniteDifferenceStaggeredStencilDerivation("NE", 2, (0,)).apply(c) + \ assert FiniteDifferenceStaggeredStencilDerivation("NE", 2, (0,)).apply(c) + \
FiniteDifferenceStaggeredStencilDerivation("NE", 2, (1,)).apply(c) == c[1, 1] - c[0, 0] FiniteDifferenceStaggeredStencilDerivation("NE", 2, (1,)).apply(c) == c[1, 1] - c[0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("NE", 3, (0,)).apply(c3) + \ assert FiniteDifferenceStaggeredStencilDerivation("NE", 3, (0,)).apply(c3) + \
FiniteDifferenceStaggeredStencilDerivation("NE", 3, (1,)).apply(c3) == c3[1, 1, 0] - c3[0, 0, 0] FiniteDifferenceStaggeredStencilDerivation("NE", 3, (1,)).apply(c3) == c3[1, 1, 0] - c3[0, 0, 0]
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
d = FiniteDifferenceStaggeredStencilDerivation("NE", 2, (0, 1)) d = FiniteDifferenceStaggeredStencilDerivation("NE", 2, (0, 1))
assert d.apply(c) == c[0,0] + c[1,1] - c[1,0] - c[0,1] assert d.apply(c) == c[0,0] + c[1,1] - c[1,0] - c[0,1]
d.visualize() d.visualize()
``` ```
%% Output %% Output
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
v3 = ps.fields("v(3): [3D]") v3 = ps.fields("v(3): [3D]")
assert FiniteDifferenceStaggeredStencilDerivation("E", 3, (0,)).apply(v3) == \ assert FiniteDifferenceStaggeredStencilDerivation("E", 3, (0,)).apply(v3) == \
sp.Matrix([v3[1,0,0](i) - v3[0,0,0](i) for i in range(*v3.index_shape)]) sp.Matrix([v3[1,0,0](i) - v3[0,0,0](i) for i in range(*v3.index_shape)])
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
``` ```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment