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

FiniteDifferenceStaggeredStencilDerivation must be applied to field access 

otherwise the index gets lost
parent ed864c0e
Branches
Tags
1 merge request!109FiniteDifferenceStaggeredStencilDerivation must be applied to field access
Pipeline #20568 passed with stage
Stage:
    in 2 minutes and 18 seconds
    Hide whitespace changes
    Inline Side-by-side
    pystencils/fd/derivation.py
    + 2
    − 8
    View file @ 37aff58b
    ...@@ -340,11 +340,5 @@ class FiniteDifferenceStaggeredStencilDerivation: ...@@ -340,11 +340,5 @@ class FiniteDifferenceStaggeredStencilDerivation:
    from pystencils.stencil import plot from pystencils.stencil import plot
    plot(pts, data=ws) plot(pts, data=ws)
    def apply(self, field): def apply(self, access: Field.Access):
    if field.index_dimensions == 0: return sum([access.get_shifted(*point) * weight for point, weight in zip(self.points, self.weights)])
    return sum([field.__getitem__(point) * weight for point, weight in zip(self.points, self.weights)])
    else:
    total = field.neighbor_vector(self.points[0]) * self.weights[0]
    for point, weight in zip(self.points[1:], self.weights[1:]):
    total += field.neighbor_vector(point) * weight
    return total
    pystencils_tests/test_fd_derivation.ipynb
    + 19
    − 18
    View file @ 37aff58b
    %% 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.center) == 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.center) == 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.center) == 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.center) == 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.center) == 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.center) == 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.center) == 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.center) == 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.center) == 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.center) == c3[0, 0, 0] - c3[0, 0, -1]
    assert FiniteDifferenceStaggeredStencilDerivation("S", 2, (0,)).apply(c) == \ assert FiniteDifferenceStaggeredStencilDerivation("S", 2, (0,)).apply(c.center) == \
    (c[1, 0] + c[1, -1] - c[-1, 0] - c[-1, -1])/4 (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.center) + \
    FiniteDifferenceStaggeredStencilDerivation("NE", 2, (1,)).apply(c) == c[1, 1] - c[0, 0] FiniteDifferenceStaggeredStencilDerivation("NE", 2, (1,)).apply(c.center) == c[1, 1] - c[0, 0]
    assert FiniteDifferenceStaggeredStencilDerivation("NE", 3, (0,)).apply(c3) + \ assert FiniteDifferenceStaggeredStencilDerivation("NE", 3, (0,)).apply(c3.center) + \
    FiniteDifferenceStaggeredStencilDerivation("NE", 3, (1,)).apply(c3) == c3[1, 1, 0] - c3[0, 0, 0] FiniteDifferenceStaggeredStencilDerivation("NE", 3, (1,)).apply(c3.center) == 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.center) == 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) == \ 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)]) assert FiniteDifferenceStaggeredStencilDerivation("E", 3, (0,)).apply(v3.center_vector[i]) == \
    v3[1,0,0](i) - v3[0,0,0](i)
    ``` ```
    %% 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.
    Finish editing this message first!
    Please register or to comment