From 71a1c19b87c78c932cb4c807e09498ed2229f995 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Mon, 3 Dec 2018 09:56:35 +0100
Subject: [PATCH] lbmpy phase-field: high density entropic model works :)
---
fd/spatial.py | 29 +++++++++++++++++++++++++++++
1 file changed, 29 insertions(+)
diff --git a/fd/spatial.py b/fd/spatial.py
index ff0a28bb5..6b147fe96 100644
--- a/fd/spatial.py
+++ b/fd/spatial.py
@@ -51,6 +51,35 @@ def fd_stencils_isotropic(indices, dx, fa):
raise NotImplementedError("Supports only derivatives up to order 2 for 1D and 2D setups")
+def fd_stencils_isotropic_high_density_code(indices, dx, fa):
+ dim = fa.field.spatial_dimensions
+ if dim == 1:
+ return fd_stencils_standard(indices, dx, fa)
+ elif dim == 2:
+ order = len(indices)
+
+ if order == 1:
+ idx = indices[0]
+ assert 0 <= idx < 2
+ other_idx = 1 if indices[0] == 0 else 0
+ weights = {-1: sp.Rational(1, 12) / dx,
+ 0: sp.Rational(1, 3) / dx,
+ 1: sp.Rational(1, 12) / dx}
+ upper_terms = sum(fa.neighbor(idx, +1).neighbor(other_idx, off) * w for off, w in weights.items())
+ lower_terms = sum(fa.neighbor(idx, -1).neighbor(other_idx, off) * w for off, w in weights.items())
+ return upper_terms - lower_terms
+ elif order == 2:
+ if indices[0] == indices[1]:
+ idx = indices[0]
+ diagonals = sp.Rational(1, 8) * sum(fa.neighbor(0, i).neighbor(1, j) for i in (-1, 1) for j in (-1, 1))
+ div_direction = sp.Rational(1, 2) * sum(fa.neighbor(idx, i) for i in (-1, 1))
+ center = - sp.Rational(3, 2) * fa
+ return (diagonals + div_direction + center) / (dx ** 2)
+ else:
+ return fd_stencils_standard(indices, dx, fa)
+ raise NotImplementedError("Supports only derivatives up to order 2 for 1D and 2D setups")
+
+
def discretize_spatial(expr, dx, stencil=fd_stencils_standard):
if isinstance(stencil, str):
if stencil == 'standard':
--
GitLab