Skip to content
Snippets Groups Projects
Commit 248a5e0d authored by Markus Holzer's avatar Markus Holzer Committed by Markus Holzer
Browse files

implemented derivation of gradient weights via rotation 

derive gradient weights of other direction with
already calculated weights of one direction
via rotation and apply them to a field.

as_matrix gives now a deprecated warning and as_array is no newly implemented to replace it. as_array works with MutableDenseNDimArray now instead of Matrix
parent 3d00d87f
No related branches found
No related tags found
1 merge request!15implemented derivation of gradient weights via rotation
Hide whitespace changes
Inline Side-by-side
pystencils/fd/derivation.py
+ 49
− 0
View file @ 248a5e0d
from collections import defaultdict from collections import defaultdict
import sympy as sp import sympy as sp
import numpy as np
from pystencils.field import Field from pystencils.field import Field
from pystencils.sympyextensions import multidimensional_sum, prod from pystencils.sympyextensions import multidimensional_sum, prod
from pystencils.utils import LinearEquationSystem, fully_contains from pystencils.utils import LinearEquationSystem, fully_contains
import warnings
class FiniteDifferenceStencilDerivation: class FiniteDifferenceStencilDerivation:
"""Derives finite difference stencils. """Derives finite difference stencils.
...@@ -169,6 +172,9 @@ class FiniteDifferenceStencilDerivation: ...@@ -169,6 +172,9 @@ class FiniteDifferenceStencilDerivation:
return sum(f.get_shifted(*offset) * weight for offset, weight in zip(self.stencil, self.weights)) return sum(f.get_shifted(*offset) * weight for offset, weight in zip(self.stencil, self.weights))
def as_matrix(self): def as_matrix(self):
warnings.warn("as_matrix is deprecated and may be removed in the near future."
"Please use as_array instead which will return an MutableDenseNDimArray."
"as_array therefore can also work in 3 dimensions", category=DeprecationWarning)
dim = len(self.stencil[0]) dim = len(self.stencil[0])
assert dim == 2 assert dim == 2
max_offset = max(max(abs(e) for e in direction) for direction in self.stencil) max_offset = max(max(abs(e) for e in direction) for direction in self.stencil)
...@@ -177,6 +183,49 @@ class FiniteDifferenceStencilDerivation: ...@@ -177,6 +183,49 @@ class FiniteDifferenceStencilDerivation:
result[max_offset - direction[1], max_offset + direction[0]] = weight result[max_offset - direction[1], max_offset + direction[0]] = weight
return result return result
def as_array(self):
dim = len(self.stencil[0])
assert (dim == 2 or dim == 3), "Only 2D or 3D matrix representations are available"
max_offset = max(max(abs(e) for e in direction) for direction in self.stencil)
shape_list = []
for i in range(dim):
shape_list.append(2 * max_offset + 1)
number_of_elements = np.prod(shape_list)
shape = tuple(shape_list)
result = sp.MutableDenseNDimArray([0] * number_of_elements, shape)
if dim == 2:
for direction, weight in zip(self.stencil, self.weights):
result[max_offset - direction[1], max_offset + direction[0]] = weight
if dim == 3:
for direction, weight in zip(self.stencil, self.weights):
result[max_offset - direction[1], max_offset + direction[0], max_offset + direction[2]] = weight
return result
def rotate_weights_and_apply(self, field_access: Field.Access, axis):
"""derive gradient weights of other direction with already calculated weights of one direction
via rotation and apply them to a field."""
dim = len(self.stencil[0])
assert (dim == 2 or dim == 3), "This function is only for 2D or 3D stencils available"
rotated_weights = np.rot90(np.array(self.as_array()).reshape(self.as_array().shape), 1, axis)
result = []
max_offset = max(max(abs(e) for e in direction) for direction in self.stencil)
if dim == 2:
for direction in self.stencil:
result.append(rotated_weights[max_offset - direction[1],
max_offset + direction[0]])
if dim == 3:
for direction in self.stencil:
result.append(rotated_weights[max_offset - direction[1],
max_offset + direction[0],
max_offset + direction[2]])
f = field_access
return sum(f.get_shifted(*offset) * weight for offset, weight in zip(self.stencil, result))
def __repr__(self): def __repr__(self):
return "Finite difference stencil of accuracy {}, isotropic error: {}".format(self.accuracy, return "Finite difference stencil of accuracy {}, isotropic error: {}".format(self.accuracy,
self.is_isotropic) self.is_isotropic)
pystencils_tests/test_fd_derivation.ipynb
+ 45
− 94
View file @ 248a5e0d
%% 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]$$ AttributeError Traceback (most recent call last)
⎡0 0 0⎤ <ipython-input-4-ea41cd8e50a0> in <module>
⎢ ⎥ ----> 1 standard_2d_00.get_stencil().as_matrix()
⎢1 -2 1⎥
⎢ ⎥ ~/pystencils/pystencils/pystencils/fd/derivation.py in as_matrix(self)
⎣0 0 0⎦ 185 for direction, weight in zip(self.stencil, self.weights):
186 result[max_offset - direction[1], max_offset + direction[0]] = weight
--> 187 result.tomatrix().tomatrix()
188 print("test")
189 if dim == 3:
~/anaconda3/envs/pystencils/lib/python3.7/site-packages/sympy/matrices/matrices.py in __getattr__(self, attr)
2139 else:
2140 raise AttributeError(
-> 2141 "%s has no attribute %s." % (self.__class__.__name__, attr))
2142
2143 def __len__(self):
AttributeError: MutableDenseMatrix has no attribute tomatrix.
%% 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
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
$$\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⎤
⎢ ⎥
⎢5/6 -5/3 5/6 ⎥
⎢ ⎥
⎣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
%% 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:
``` python
type(isotropic_2d_00_res.as_matrix())
```
%% Cell type:code id: tags:
``` python
type(expected_result)
```
%% 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
$$\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⎤
⎢ ⎥
⎢2/3 -10/3 2/3⎥
⎢ ⎥
⎣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
``` ```
......
pystencils_tests/test_fd_derivation_via_rotation.ipynb 0 → 100644
+ 195
− 0
View file @ 248a5e0d
%% Cell type:code id: tags:
``` python
from pystencils.session import *
from pystencils.fd.derivation import *
```
%% Cell type:markdown id: tags:
# 2D isotropic stencils
%% Cell type:code id: tags:
``` python
stencil_2D = ((0, 0),
(0, 1), (0, -1), (-1, 0), (1, 0),
(-1, 1), (1, 1), (-1, -1), (1, -1))
f = ps.fields("f: [2D]")
```
%% Cell type:code id: tags:
``` python
Isotropic_Gradient = list()
deriv_x = FiniteDifferenceStencilDerivation((0, ), stencil_2D)
deriv_x.assume_symmetric(0, anti_symmetric=True)
deriv_x.assume_symmetric(1, anti_symmetric=False)
deriv_x.set_weight((0, 0), sp.Integer(0))
deriv_x.set_weight((1, 0), sp.Rational(1, 3))
deriv_x.set_weight((1, 1), sp.Rational(1, 12))
res_x = deriv_x.get_stencil(isotropic=True)
```
%% Cell type:code id: tags:
``` python
Isotropic_Gradient = list()
deriv_y = FiniteDifferenceStencilDerivation((1, ), stencil_2D)
deriv_y.assume_symmetric(0, anti_symmetric=False)
deriv_y.assume_symmetric(1, anti_symmetric=True)
deriv_y.set_weight((0, 0), sp.Integer(0))
deriv_y.set_weight((0, 1), sp.Rational(1, 3))
deriv_y.set_weight((1, 1), sp.Rational(1, 12))
res_y = deriv_y.get_stencil(isotropic=True)
```
%% Cell type:code id: tags:
``` python
assert res_x.apply(f.center) - res_y.rotate_weights_and_apply(f.center, (1, 0)) == 0
assert res_y.apply(f.center) - res_x.rotate_weights_and_apply(f.center, (0, 1)) == 0
```
%% Cell type:markdown id: tags:
# 3D isotropic stencils
%% Cell type:code id: tags:
``` python
stencil_3D = ((0, 0, 0),
(0, 1, 0), (0, -1, 0), (-1, 0, 0), (1, 0, 0), (0, 0, 1), (0, 0, -1),
(-1, 1, 0), (1, 1, 0), (-1, -1, 0), (1, -1, 0),
(0, 1, 1), (0, -1, 1), (-1, 0, 1), (1, 0, 1),
(0, 1, -1), (0, -1, -1), (-1, 0, -1), (1, 0, -1),
(1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1), (1, 1, -1), (-1, 1, -1), (1, -1, -1),
(-1, -1, -1))
f = ps.fields("f: [3D]")
```
%% Cell type:code id: tags:
``` python
deriv_x = FiniteDifferenceStencilDerivation((0, ), stencil_3D)
deriv_x.assume_symmetric(0, anti_symmetric=True)
deriv_x.assume_symmetric(1, anti_symmetric=False)
deriv_x.assume_symmetric(2, anti_symmetric=False)
deriv_x.set_weight((0, 0, 0), sp.Integer(0))
deriv_x.set_weight((1, 1, 1), sp.Rational(1, 3360))
res_x = deriv_x.get_stencil(isotropic=True)
```
%% Cell type:code id: tags:
``` python
deriv_y = FiniteDifferenceStencilDerivation((1, ), stencil_3D)
deriv_y.assume_symmetric(0, anti_symmetric=False)
deriv_y.assume_symmetric(1, anti_symmetric=True)
deriv_y.assume_symmetric(2, anti_symmetric=False)
deriv_y.set_weight((0, 0, 0), sp.Integer(0))
deriv_y.set_weight((1, 1, 1), sp.Rational(1, 3360))
res_y = deriv_y.get_stencil(isotropic=True)
```
%% Cell type:code id: tags:
``` python
deriv_z = FiniteDifferenceStencilDerivation((2, ), stencil_3D)
deriv_z.assume_symmetric(0, anti_symmetric=False)
deriv_z.assume_symmetric(1, anti_symmetric=False)
deriv_z.assume_symmetric(2, anti_symmetric=True)
deriv_z.set_weight((0, 0, 0), sp.Integer(0))
deriv_z.set_weight((1, 1, 1), sp.Rational(1, 3360))
res_z = deriv_z.get_stencil(isotropic=True)
```
%% Cell type:code id: tags:
``` python
assert res_x.apply(f.center) - res_y.rotate_weights_and_apply(f.center, (1, 0)) == 0
assert res_x.apply(f.center) - res_z.rotate_weights_and_apply(f.center, (2, 1)) == 0
assert res_y.apply(f.center) - res_x.rotate_weights_and_apply(f.center, (0, 1)) == 0
assert res_y.apply(f.center) - res_z.rotate_weights_and_apply(f.center, (0, 2)) == 0
assert res_z.apply(f.center) - res_x.rotate_weights_and_apply(f.center, (1, 2)) == 0
assert res_z.apply(f.center) - res_y.rotate_weights_and_apply(f.center, (2, 0)) == 0
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment