Commit dc8453b3 authored by Martin Bauer's avatar Martin Bauer
Browse files

Merge branch 'derivation_gradient' into 'master'

implemented derivation of gradient weights via rotation

See merge request pycodegen/pystencils!15
parents ddb86435 248a5e0d
from collections import defaultdict
import sympy as sp
import numpy as np
from pystencils.field import Field
from pystencils.sympyextensions import multidimensional_sum, prod
from pystencils.utils import LinearEquationSystem, fully_contains
import warnings
class FiniteDifferenceStencilDerivation:
"""Derives finite difference stencils.
......@@ -169,6 +172,9 @@ class FiniteDifferenceStencilDerivation:
return sum(f.get_shifted(*offset) * weight for offset, weight in zip(self.stencil, self.weights))
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])
assert dim == 2
max_offset = max(max(abs(e) for e in direction) for direction in self.stencil)
......@@ -177,6 +183,49 @@ class FiniteDifferenceStencilDerivation:
result[max_offset - direction[1], max_offset + direction[0]] = weight
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):
return "Finite difference stencil of accuracy {}, isotropic error: {}".format(self.accuracy,
self.is_isotropic)
%% Cell type:code id: tags:
``` python
from pystencils.session import *
from pystencils.fd.derivation import *
```
%% Cell type:markdown id: tags:
# 2D standard stencils
%% Cell type:code id: tags:
``` python
stencil = [(-1, 0), (1, 0), (0, -1), (0, 1), (0, 0)]
standard_2d_00 = FiniteDifferenceStencilDerivation((0,0), stencil)
f = ps.fields("f: [2D]")
standard_2d_00_res = standard_2d_00.get_stencil()
res = standard_2d_00_res.apply(f.center)
expected = f[-1, 0] - 2 * f[0, 0] + f[1, 0]
assert res == expected
```
%% Cell type:code id: tags:
``` python
assert standard_2d_00_res.accuracy == 2
assert not standard_2d_00_res.is_isotropic
standard_2d_00_res
```
%%%% Output: execute_result
Finite difference stencil of accuracy 2, isotropic error: False
%% Cell type:code id: tags:
``` python
standard_2d_00.get_stencil().as_matrix()
```
%%%% Output: execute_result
%%%% Output: error
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAGkAAABLCAMAAABEOFk1AAAAPFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAo1xBWAAAAE3RSTlMAMquZdlQiEEDpMESJZs3d77ts6LWnRAAAAAlwSFlzAAAOxAAADsQBlSsOGwAAAj9JREFUWAntmV17gyAMRlFQ1olf4///1wl0kpc+Qlba3QxvBnKas6TUaio6649evOtYg0CIzkp1HMO7RGJ24XvrTOPbJDHwgCajZtPF1csRk6MYmpajgnotZ8jkAAPTvLkkzHKZyn2BySEGpt1vv8HqgorJIQYm602TLe1CJocYNWkrXTKTVfmcmFyCUdNojVN0JROTSzA0+Zy6ILzOawy5l7gEo6Yk3UsVk0swahLhLRyYO6LIYTgw7bvLYy7vch6H4cCEH7XL6gkmhxiYxOquRlv5wsfkAEOTNsfVvSwSTA4wNF0XrH6lmWpq2KrXqocVaDsC6yG0lPvOuEz6l6XVG7bSLVi0aXcrZewUzzyMSDgw6d70xa/BGE36f2rzt6PxbBxhODAdUPkLN0baVjfuc88PJFyNafH59zbzTr3IFLJbc/V+pWnK3h2+0rRmH1qfNOmFHPc3R/rb3rhLklFi+rh9ngBZOs9lBiovolv564bPub8zDU40Zp4hSbiaXS46/2yinjOp3J5NyjiuxhjZ+89vsvQzJeEwJ9lvdpHzD1f4u4SGxrUJwqGpELpquZlqyteq16qHFfgHO4I2GTF5nDE5imH1oMmIsWHG5AADE7YqIDZMmBxiYMImIwSHCZNDDEzYZITgMGFyiFFT0iaD4HTC5BKMmpImIw0OYyaXYGhqPctSb7P1LGHTnRNoRp5nHwaA0b13dEZbz/KhXLkTWL0cWbvWTDUVDNX7q1+OtfthV6lcp6kmmaNFFATiG+4bKVe2WSLVAAAAAElFTkSuQmCC)
$$\left[\begin{matrix}0 & 0 & 0\\1 & -2 & 1\\0 & 0 & 0\end{matrix}\right]$$
⎡0 0 0⎤
⎢ ⎥
⎢1 -2 1⎥
⎢ ⎥
⎣0 0 0⎦
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
<ipython-input-4-ea41cd8e50a0> in <module>
----> 1 standard_2d_00.get_stencil().as_matrix()
~/pystencils/pystencils/pystencils/fd/derivation.py in as_matrix(self)
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:
# 2D isotropic stencils
## second x-derivative
%% Cell type:code id: tags:
``` python
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_res = isotropic_2d_00.get_stencil(isotropic=True)
assert isotropic_2d_00_res.is_isotropic
assert isotropic_2d_00_res.accuracy == 2
isotropic_2d_00_res
```
%%%% Output: execute_result
Finite difference stencil of accuracy 2, isotropic error: True
%% Cell type:code id: tags:
``` python
isotropic_2d_00_res.as_matrix()
```
%%%% Output: execute_result
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAIIAAABNCAMAAABdY0yuAAAAYFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACACGXEAAAAH3RSTlMAMquZdlQiEEDpMES73YnNZu+F48P7apvf1/Fc0cls39hR/gAAAAlwSFlzAAAOxAAADsQBlSsOGwAAA8dJREFUaAXtmumW2yAMhcFrU6/pdJ0ufv+3LNjGBnQBkcn0NOfEPxIbpMsXGQMyEcI/ZO+XsK55bhwrWZYDq0nXiOcWsJLLehSHZHcLghA8N9tq2FoWQi5lpY763yPMut1i0QjN0fp6YoO6NdErnhuxqp8IKqzPKOi+BaNQ3/ZQ8tyIFUBoymEpqmjfR5U8N2AFEJD+e5Y9LkI7WYd8W4zcKGxDtvXJFz+cgi6HhTnZLV2EoPt7VhCEfTofy75sMxsex2qOu0BxD8FM56OauIoprufXTp1orn6hfR0Q9xDMpD8oBDKB2nL0vFLrjTYxnGyzpCceQuiEaBb1wT+u53oj6LQjuOIBBC1SL1mdYamrKtV9zrWCJR5BGBJRdX9rs4zmLroV9tWJYImHEcbSdk6et+ttS9y7A8EWDyLMeQRCbAjxp9IgOOIhhE6HNWvknXTrvCi44gShHmQhm08vdV33Wd3x60U2317iNwyKewh6Ov+sxvAv6zge1/NqO+XyPQodEPcQtGxXy6iS17S55LkBK4SQ1QUMAbPngP6FEFRulRpjjobPk47lBqwAQqPGzzk/t+a5ASuAoH9Y5hRlYsFz86wAgr4J21hnlFnfPDdgBRCuaqJu8qYozchzA1YAYVRRsMdwVgyUC8sNWAEEMY6lGp6zD54btUII2Y2/zeFxEe6dyny4fLRCaRIN9fol5+C5Ias/l5vHoRzAmC3sC+mcBEmOVVkwplgijhCSOQkCEKVaNZXpqYWKAwROTgIg9FgyJzsQEAcInJwEIOiidBSAOEDg5CQYob0mVztAnCKwchKIUE3J5AeJUwRWTgIRGDcCiVMEVk4SQGiWVG4LEh6AwMlJKMLaD5olnk0JAcQBQtfLvs58vaGQhkb28zU1NgFxgCD7H683rKB//vr9Wni7CjRYQBwgpHN0qryWmKQ1UL0XE6sngo7MMwrPKOxPCO4LZOPEWMe/eW7ECnRHsHESb3qr5bkBK4DAae+eNo+LcO9U5nzlcGQa5oQfcOMRXsAeFuZkF/8fbwTcOOEFgyQpxA2Ke1EIbJwQLVRAkxTPKiDuIZi1grdx4mnBS5CkELttreCJhxDcjROiBQpAkkKsdgRXPICgna2NE6IFCkCSQqzOFZMlHkGwNk6IFi1ASQqxOhEs8TBC5ks3lKSEEWzxIIKzcUK0QAFIUoiViYIjHkJwN06IFigASQqx2hFc8Q1hHTK3v1Zt07mcsndlOvV6I5UB+eLHX6ta/QenqlKPijDT+XWFIr8iWlCXYzwDouLrX6uqSvwFDr5NBGUlpz0AAAAASUVORK5CYII=)
$$\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:
``` python
plt.figure(figsize=(2,2))
isotropic_2d_00_res.visualize()
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAKQAAACaCAYAAAAnxeOcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAE+FJREFUeJztnXtwXNV5wH9nd7UPrWTJ8soWsh6uATPBIQ/KKzEloSngdcOAE9M0E6cUApmUuJPAJI1pnLSJA2RqbEKGkjKEAs0wDiUYOoAX8iyEDEljYGyS1ji2gm1JfuhaWkkrrR6rPf1jtbZkr/Z1X+de3d+MZ6x9nPt9+n66j71nvyOklHh4qILP7gA8PGYSsGpDIqHVARcBUUBU8FYJpIBdMh4bMSM2IxAJbQHwp0AtleVXiCxwBNgt47Gs3tjMQiS0ZcB5QE2Fb53iVH6zDtHCikO2SGg3AF8HgjqGGQe+IeOxp42JyjhEQvs08BUqL0wpeoCbZTz2jsHj6kIktBBwP3ClzqEOAjfJeKwn/4Dph2yR0JqBb6JPRoAQsFkktJj+qIxDJLR2YBPGywiwFNhswrh6+TT6ZQToBP5p5gNWnEP+mYHb8QOrDBrLKK4wefyLREKrNXkblfIhA8daJRKaP/+DFUI2GDxeo8Hj6cXo/E7HB9SbvI1KMTLnALnrCsAaIfWe4Js9nhNQLWfT4rH/Y5/hAR+f/3AHazvOZf9uveeZzsKtuevIy34hw9Esm5/s4ZKrh+0OxXLcmruOvOwXsiYITUum7A7DFtyau4687BfSw2MGnpAeSuEJ6aEUlt3LLsrGtUs5uDdMb1eQa9Yn+ejNQ3aHZBluzb3KvNQQ8tvP9JR+kUtxa+5V5uUdsj2UQo09pBu5Y3U7e1+PFHzu3Peluf8nhy2OyBHYI+Sa5hUVvX5n3z6TIjGPbS+WFq7Y78FpORtUU3uEfHLffjaubaO3K8SW5w5yznsn2PNqmIe/vphAjWThkgx3fv8INe65m1aQnX372PWzCD/ctohsVnDtLQNc+fGU3WFVhUE1teccstCtpZZlk2x5/jD3vXSYJR2TvPJsnS2xWcl4WrDjwSbu3tHNtsRhx8oIhtXUHiEL3Vpa3DZFuDY3fd3vlwgXXm9JCb1dZ5Psy00y3vNqhGA4y6Yb2vjaJ1rRev0lRlAXg2qqXtW7D9Tw5stRLr82xehwHccOdTIyVIcbvh2ZTkWR0sfIcCNSQv9xP8cOBfnWU92sXj/IY3cpNRu+JJmMn6S2iKOHOslOze3SzJqWQK2r7FTSx723tfDlB48QDEF2KsPkRIhk31kMapJowwDjo/XgrLqdJDXYhJQCgSSdqqOuIct5F6YJhuDiq0Z56oEmu0Msi6TWSGailfF0FCkFPn8G4Sv8ZbTTa1oCdfaQmUm46+az+NSXTrDs/EkAQpEx/IFJpBRksz5SySae2HKfEMLsWdrGk5kMMDkeBkBKH6nkQlZemqb7QBCZhbffCNHSMWlzlCURQnyefW+sYmy0LvfHJST1jf2IAnN2C9W0BPbtIU+/teTzw4E9YbZvXcT2rYuI35jkqk8OU7/wBIPaEqQUSCloWvJbwHnzB1ODDUh5qmoT4xGaWnxctjrF7fF2fALueOCojRGWywvAd5k5azy6YBAov6ZFsE/IQreW1tx45v3O2vphBrXFCAFSCk4cvcyK8Aynti4FUpAabEIISbShH38gw7oNSdZtSNodXgUsRspTR9ZI3RA+f+4Ev9yaFkGdQ/Zc+HySBU0akbpBWjq6ph+dEsJhl+HB8DiNzRoA4egwjbETBQ9zCiOEuAT4Df7ABC0dXQTDo9Qv7DdyG1YUNa17hPqFSZqWHCcQzHD95zZMP6qKlGMWbEP/71AnJ2UEjYs+8l8EghkWt3VTEzTivHc8/x8rCrrL0NGWr3yR3BfMQQ0pXzd5/P0yHhs0eRtFmSmjlLIZY3N+S8Zj1gkp47E/ADsMGu5pGY8dkFIeQhEpZTy2G0iYNHwG2GrS2GVRQEaAHwC9Bgw/Adw3a3sW9fbxAR+c/ldHdc2mfgW8NrP5khCig1x/GAC/lNKWxkzTnRcuBz5AqWZT27duonnpbv7ir58rMuQUcBR4ScZjfzQy1kqYQ8bccwmtAVhN9c2meoEXZTx2aNa4Tu8PqYqU5SKEkMBjUsqb7I6lGMVkNBO7z790o9Lh2y3YJSO4QEjwpDQSO2UElwgJnpRGYLeM4CIhwZNSDyrICC4TEjwpq0EVGcGFQoInZSWoJCO4VEjwpCwH1WQEFwsJnpTFUFFGcLmQ4ElZCFVlhHkgJHhSzkRlGWGeCAmelKC+jDCPhIT5LaUTZIR5JiTMTymdIiPMQyFhfknpJBlhngoJ80NKp8kI81hIcLeUTpQR5rmQ4E4pnSojeEIC7pLSyTKCJ+RJ3CCl02UET8hZOFlKN8gInpBn4EQp3SIjeEIWxElSuklG8IScEydI6TYZwROyKCpL6UYZwROyJCpK6VYZwROyLFSS0s0ygidk2aggpdtlBBs76E43aCrWpGhSxmNVrUpvFlLKQ0KITnK9hKaEELN6CYmEJoAgxZpN1YTAH/CLhBYusqmsjMcmZj7gBBlFQgsy905OAhMyHivaTMryZlMioV0L3Equa1Yp3gYelvFYsU5hlnN6gyt29q0A/gG4lFJ/5N37VxCpG2JRS6l+4v3A88C/sKb5/Sgq4/Qf1kYgDjSWePkE8Cpwt4zHCi69Z+lhRyS0VcAWypOR6ddtEQntA+ZFVTmzDt/h6BRSPg6swtgjThPwNxzc+z0UlXGazcAnKS0j5I4efw48Pn2EPAOrz4Ouo7LekEy//joTYtHFSSnf80HoOXCJKQs7jafDTE58Br9fSRlFQosA11Tx1qXARYWesFrI9irf12ZoFAYhpTzErZv/EYCeAysMlXI8Haavp4PauimeO3q+cQMbSjNQejWkwhSsqdVCVrs9tVYcm8nSswdp6cytDmGUlHkZfb4pWpcfAFRdA1GPPwVzsr/QPV0Bbr+mk7azc1eVmx7vPWMRR9UJ1GRo6ezi6MHl9BxYwdKz95W15Eeh3KMLak6T0VnorKf9QgK866I039huRBN1+6hWypm5n7lndCY66qnGB+P73ozwxavbeeirMdRuEV6cvJRQ/uE7n/v37mzheLfzZQRd9bRfyObWKR7Z1cV9Lx5mUPPzi6edvXB7JVLmc//2M30MHFvAmy9nHS+jznraL2QwLKmtyy3uveqjKbreqvaqTR3KlTIYlvj9IbQjHVx4ZZbjhwYsjNIcdNbTfiFTyVMxvPVahNblyi/RWxblSNl/LHLynPGd/0u6Ined9bT/ouaN/46wfWuMUCTL4rZJbvmGZndIhlHoQifPeDrMazvbee4RCNdOsLgt4IrcddbTfiGvuH6EK64fmfN5KaH3j21CNP+K3L33D1oXnAGcLiVAdspPX08H7//QFH95k+POGYUQFwA/5ewLfsrWRJBg6NREkFL1LIH9Qs7F5EQNqcFGRocb6PrdcnKxjkzPtlGHh3/dxJL20r/HWOshtN4OAMbTUXy+LIvbDpKZKP3ex+9pE2seUOnc+kIgSir5V/R1+/EHJqhrHKCuoaK1sQuhrpADx1qYGI8AkJ3Kn5dEgXfsCqkgzz4E1332zMfvvQ32vVn4PWe/B77ykI+jh5aXtY2f/+dvdURoHpLcESwzGSLZ1+IOIdc0r5jzuSf/cIxUsgmfT5C7ABuUUjZZFlsZiIR2G/CFM574zo9n/5z/0Hsmny0yiWln36nzzSd+/yEZj5WarmYZQoirgWcRQiBEiGBkhPrG3CcExeo5M6c5sF/InX372PWzCD/ctohsVnDtLQNc+fHUyefrGwcJRd4htzpqrV1h6mLmHZhs1k+4dpix0XrufQH+7c4x/AGJzy/Z+PARmpc64bZpF7CXWGuCls5P4A+cirlUPUtgv5DjacGOB5u4e0c3wTlOk5qWDEgpb7Q2MIM4/XZg9/4VCJ+kpbOL7NRyvvS9MO3n7uOFRxfwwqMN/O2mfrtDLoWUcj9woUhoy4B1s54sp55FsP9zyD2vRgiGs2y6oY2vfaIVrVfVmS2VU+zedKAmQ+vyLny+3OeUo8M+Ot81McdIzkFnPe0Xsv+4n2OHgnzrqW5Wrx/ksbtidodkCOVMlAjUZBjq7+buz0DiPxaz4v1jFkdpPDrrab+QdQ1ZzrswTTAEF181Svf+oN0h6aaSWTvnXzLKd3/axXW3wg/uWW7KzHMr0VlP+4VceWma7gNBZBbefiNES4ezb59VIuPEWG5+WqAmQ3PbEYJh42eeW43Oetp/UdPYnOWy1Sluj7fjE3DHA8p8vFExlc5n3LsrxL9vbsbng5pQltvvP0hmsrOi+ZSqobOe9gsJsG5DknUbknaHoYtqJte+5/IxvvPS7K+DZiYrn+SrGjrqaf8h2w0YOdO7mkm+LsJqITNVvk/d88rfvdZu+NcOTpcyqan6YbmeuAq6YLWQXVW+T8kZMUKIS9i+daMpXzvIS5lKwvqVx1XoulaAY0C1M3sKumB1kj+i8r3dJPC0CbHo4mSvnd//WqN1+cumbCRQk6Gx+Ztks6BIK8CZTPcfqqY2/wvsLvSEpRc1Mh7bLRLaLcCNwEpKNJsCfg88KuOxt6yIr1xmNX4aH2sWCa0V+DvgA+Tut899JTI8ADBOrnfPXGSBI8ALtHQ+BjzCHA2uFOAe4Di53j5LKN5sahj4JfCvczWdsrzZlNPR24VMCCGBx6SUN1X4vlkNrhST0jCUOgSojp0t8VToT2kFrkzKDFTozzgfpHRdQmaggox53C6lq5IxA5VkzONmKV2TiBmoKGMet0rpiiTMQGUZ87hRSscnYAZOkDGP26R0dPBm4CQZ87hJSscGbgZOlDGPW6R0ZNBm4GQZ87hBSscFbAZukDGP06V0VLBm4CYZ8zhZSscEagZulDGPU6V0RJBm4GYZ8zhRSuUDNIP5IGMep0mpdHBmMJ9kzOMkKZUNzAzmo4x5nCKlkkGZwXyWMY8TpFQuIDPwZDyF6lIqFYwZeDKeicpSKhOIGXgyzo2qUioRhBl4MpZGRSltD8AMPBnLRzUpXSekJ2PlqCSlq4T0ZKweVaR0jZCejPpRQUpXCOnJaBx2S+l4IT0ZjcdOKS1pNiUSWgi4DlgFVLTC/DQp4FfAs9Mt4HLjKiKjSGgR4Hpy3c+iRV/80hM3s6jlD1z0kV8WeZUEeoEXZDz2G8MCrZBiDa5EQusEPgacR/EudoXIAj3A8zIe2zVrmxYJ+QhwuQFDvSLjsVtBKRn9wOPAxWW9oXv/CiJ1QyxqKbcZ/FdlPPajauPTSyEpRUI7B9gOLNA5vAS+LOOx5/IPmL4rFgntvRgjI8AVIqG9WxUZp7mYcmWsjttMHLskcxy+16NfRsj10fzczAesODe4wNDRXv/5WtSREeDdJo+/VCQ0W1fAPV1KpDSypudMn/IA1gipf2WuZF+Mvp42xkYjvPzMN1FHRjAiPzW2UZRZUv7Pjz/GeDrEkXf+hIkxIxaWP3kOqv5VdnbKx8jQQibGImi97fj9owrJWD7pkVr6eltz/08t4MSRs8hmHbUIzSwp+3o6mcrUMDRg6N7bfiGHB3x8/sMdrO04l/27z9wTpIZy5ypS5ooXqT8shDDir9JaJsdDjI+e+oQhPVLPUH+gaO5q0oQQp9pJj4/WMZWZveJrqZoWwX4hw9Esm5/s4ZKrh894TkpIJZtOyiiEJJ1qQ5UVyCohumCIXH/xHKHICNEFk3Pmri4rAWblkhpsmPWKYjUtgf1C1gShaUnhBXjGRqJkpwIIIQnUTNDYfJTP3f0FKWW1a6PYhz8wRTA8CuSKWdc4UDR3RZFSPsEFq35OdMEAQmSRwMhg06wVx3TkpfaeRogs4doU9Qv7CUXya0k7qoCzqG/sZzwdRfiyhGtH7Q6naiLRMSJRjQWLTjA6XM/YSB1Sill7zSpRW8hwNE04mrY7DMMI1abxBzLU1icduajm6fh8krqGIeoahowaUm0h3YYQcNayapfXmxeoIeTGtUs5uDdMb1eQa9Yn+ejNhv3FKY9bc68yLzWE/PYzPXaHYBtuzb3KvOy/yvbwmIEae0g3csfqdva+Hin43LnvS3P/Tw5bHJEjsEfINc0rKnr9zr59JkViHtteLC1csd+D03I2qKb2CPnkvv1sXNtGb1eILc8d5Jz3TrDn1TAPf30xgRrJwiUZ7vz+EWqccjetSnb27WPXzyL8cNsislnBtbcMcOXHU3aHVRUG1dSec8hCt5Zalk2y5fnD3PfSYZZ0TPLKs9XMLHcW42nBjgebuHtHN9sShx0rIxhWU3v2kIVuLS1uO/Wz3y+xv4mC+ex5NUIwnGXTDW2EIlm+cN8xYq3OvBNlUE3Vq3r3gRrefDnK5dc6d29RLv3H/Rw7FORbT3Wzev0gj90VszskU6igpmoJmUr6uPe2Fr784BGCzpthVjF1DVnOuzBNMAQXXzVK9373nTRXWFN1hMxMwl03n8WnvnSCZedP2h2OJay8NE33gSAyC2+/EaKlw115V1FT+z6HPP3Wks8PB/aE2b51Edu3LiJ+Y5KrPumkeYKV09ic5bLVKW6Pt+MTcMcD5X4TUU0MqKl9Qha6tbTmRnfcx62EdRuSrNuQtDsMQzCgplYcsrOlX2LreHox/4vt6uVsdDwnf4dWCHlC8fH0YnY8GUC1I4eROaeBk5OVrRDyZWDcoLHGgWItSOzgF5g7i/0VGY+NlX6ZpfzEwLF+IeOxk78/04WU8dgQ8PfAcZ1DHQM2TI+nDDIe6wO+CGgmDL8L2GTCuHp5EngCmCj1wiJI4DXgn2c+aElvHwCR0ARwDtU3m9ov4zFrgq0CkdB85PIr3myqPLLAURmPHTNgLNMQCa0OWEZ1zaZ6p/+YZ49plZAeHuWgzgfjHh54Qnoohiekh1J4Qnooxf8DQobOZydRBWoAAAAASUVORK5CYII=)
%% Cell type:code id: tags:
``` python
expected_result = sp.Matrix([[1, -2, 1], [10, -20, 10], [1, -2, 1]]) / 12
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:
## Isotropic laplacian
%% Cell type:code id: tags:
``` python
isotropic_2d_11 = FiniteDifferenceStencilDerivation((1,1), stencil)
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
```
%%%% Output: execute_result
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAHoAAABNCAMAAABe+IDDAAAAVFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACwjqETAAAAG3RSTlMAMquZdlQiEEDpMES73Ylmze/h0WzfvadcyfFyvev/AAAACXBIWXMAAA7EAAAOxAGVKw4bAAADhklEQVRoBe2a2XqlIBCE8YCajHuSWX3/9xyWg9LlQvfEzPLN8cKNpn5B0AJVKl2KNj3K7fe61XUuaE0/Ey+07tbI7F5fKXUbs2H3gB3xYvbLzUc0EnRn0cVsV8wlineB6DJrY5fyF9CNUtVsV8wlogfHu81w3TGVKaZUOfNvNhEv343uDPsq1bXoXvPJ16IHCflSdNMr1RTscp9XeCnpXMVYlmXLb2ZEHJpZpbv5xm83k++f3EKDOKC5KlfEPdBKhUerXfMqVBaO0X9Thfe9GXhldlG90Td+31JEHEs9Nqqa2Ghtr1LzvQUVB7SxL+2a36+1fZgNzIahFIgDegovbXaxbSC/1CAO6Lk0RmK2bB1N7Ec4iFN0NbvXgeQpbkb27UFxiq692RE4HkmFozhFq4AW9C7nzbjtA8QBPToqu9T+Plcz90pBHNDNi27bV3YDd2Z4mJ5Czza9sS3lZAFxQPdvuv/MRxeW9uVrGDUMllsGL3+AB3FASz29h4Qu0Tk7fmqMQRzRQk+/ouvZ9W94avjkZeUuLhkwANqFSTy9l/WlDuOfLtfkEvEdtMTTp2j3BsvmTQK2aJGnT9Fu1Ddlnm2p+AYt8/QrOtzrzCOBiCNa6OlXtMq3cEXFAS319B4djL2xQ6Dzfg3igJZ5eg9ejL3p0zvp0+gKxAFNYz/26N9B12OysO3Jfu25Uj89f7onokvfz0POLlniDklND2JAtJFvzzCXkgZ/7D7ea6GnjxfHy0ajAC309JHMywZRiBZ5+khWvKEARAHaifE9/YLmZ0vEt2iBp0/RvGxp1AYt8PQpmZeNRG3Qf7LCBZ4+LTYzWzJgoKUWevqI5mXDKIpWXVW0w8SfJ7h/BQhDgVw2EAf0t9fvP278qfX4FcANBfLZQBzQwiEumPp4Aw63ZPz8TjQ19YfImHAh2kkmpj4SDrcXoxNTf4iMCdeiM1YwQsP2UjQx9ZSzc3Qlmpr6HRg9dY4mE/U04/YITP02AM4Qcehci5+HPAeHYOoPouJpEAd0jPod2wf68RWA2c6oqc9menwFsK88Mtw7nXLbr8/E1O8HLGdBHNAwUb/kOt5JTf1xlE8BcYrGifqMlk0mpv48HMUpGifqz7VCKrvCUZyi/6+vAH6aw02jN3ZCv2T/XmNHDO4rQM76x9u2ii+/xNTuBxVj3FS6KnUv+M7Fs/4RvYr7X2KMUT8Bqn5A97SBuwoAAAAASUVORK5CYII=)
$$\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:
``` python
expected_result = sp.Matrix([[1, 4, 1], [4, -20, 4], [1, 4, 1]]) / 6
assert iso_laplacian == expected_result
```
......
%% 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
```
Markdown is supported
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