From db66c16bba4b774ea600bef3d94f7b5af8bbdba6 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Fri, 11 Jan 2019 14:49:13 +0100
Subject: [PATCH] Added forth order FD stencil for 2D, first derivatives
---
fd/spatial.py | 42 ++++++++++++++++++++++++++++++++++++++++++
1 file changed, 42 insertions(+)
diff --git a/fd/spatial.py b/fd/spatial.py
index 43e88eaea..610756f12 100644
--- a/fd/spatial.py
+++ b/fd/spatial.py
@@ -1,8 +1,11 @@
+from typing import Tuple
import sympy as sp
from functools import partial
+from pystencils.cache import memorycache
from pystencils import AssignmentCollection, Field
from pystencils.fd import Diff
from .derivative import diff_args
+from .derivation import FiniteDifferenceStencilDerivation
def fd_stencils_standard(indices, dx, fa):
@@ -51,6 +54,19 @@ def fd_stencils_isotropic(indices, dx, fa):
raise NotImplementedError("Supports only derivatives up to order 2 for 1D and 2D setups")
+def fd_stencils_forth_order_isotropic(indices, dx, fa):
+ order = len(indices)
+ if order != 1:
+ raise NotImplementedError("Forth order finite difference discretization is "
+ "currently only supported for first derivatives")
+ dim = indices[0]
+ if dim not in (0, 1):
+ raise NotImplementedError("Forth order finite difference discretization is only implemented for 2D")
+
+ stencils = forth_order_2d_derivation()
+ return stencils[dim].apply(fa) / dx
+
+
def fd_stencils_isotropic_high_density_code(indices, dx, fa):
dim = fa.field.spatial_dimensions
if dim == 1:
@@ -106,3 +122,29 @@ def discretize_spatial(expr, dx, stencil=fd_stencils_standard):
else:
new_args = [discretize_spatial(a, dx, stencil) for a in expr.args]
return expr.func(*new_args) if new_args else expr
+
+
+# -------------------------------------- special stencils --------------------------------------------------------------
+
+@memorycache(maxsize=1)
+def forth_order_2d_derivation() -> Tuple[FiniteDifferenceStencilDerivation.Result, ...]:
+ # Symmetry, isotropy and 4th order conditions are not enough to fully specify the stencil
+ # one weight has to be specifically set to a somewhat arbitrary value
+ second_neighbor_weight = sp.Rational(1, 10)
+ second_neighbor_stencil = [(i, j)
+ for i in (-2, -1, 0, 1, 2)
+ for j in (-2, -1, 0, 1, 2)
+ ]
+ x_diff = FiniteDifferenceStencilDerivation((0,), second_neighbor_stencil)
+ x_diff.set_weight((2, 0), second_neighbor_weight)
+ x_diff.assume_symmetric(0, anti_symmetric=True)
+ x_diff.assume_symmetric(1)
+ x_diff_stencil = x_diff.get_stencil(isotropic=True)
+
+ y_diff = FiniteDifferenceStencilDerivation((1,), second_neighbor_stencil)
+ y_diff.set_weight((0, 2), second_neighbor_weight)
+ y_diff.assume_symmetric(1, anti_symmetric=True)
+ y_diff.assume_symmetric(0)
+ y_diff_stencil = y_diff.get_stencil(isotropic=True)
+
+ return x_diff_stencil, y_diff_stencil
--
GitLab