From 5b15fd470c9ddb4c1bb4b7ad459b51de653d47a2 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Thu, 28 Jul 2016 14:16:29 +0200
Subject: [PATCH] finite differences
---
hysop/numerics/finite_differences.py | 50 +++++++++++++++++++++++++++-
hysop/numerics/stencil.py | 13 ++++++--
2 files changed, 59 insertions(+), 4 deletions(-)
diff --git a/hysop/numerics/finite_differences.py b/hysop/numerics/finite_differences.py
index 6ae41654b..7c0d3e221 100644
--- a/hysop/numerics/finite_differences.py
+++ b/hysop/numerics/finite_differences.py
@@ -49,6 +49,7 @@ Note
from abc import ABCMeta, abstractmethod
from hysop.constants import debug
from hysop.numerics.stencil import Stencil, StencilGenerator
+from hysop.numerics.stencil import CenteredStencil, CenteredStencilGenerator
import hysop.tools.numpywrappers as npw
import numpy as np
@@ -215,6 +216,45 @@ class FiniteDifference(object):
size as result.
"""
+class FDkCp(FiniteDifference):
+ """
+ k-th derivative, centered scheme, p-th order.
+
+ """
+
+ def __init__(self,k,p,*args,**kargs):
+ super(FDkCp,self).__init__(*args,**kargs)
+
+ csg = CenteredStencilGenerator()
+ stencil = csg.generate_exact_stencil(order=p,derivative=k)
+
+ lghosts = stencil.L
+ rghosts = stencil.R
+
+ self.stencil = stencil
+ self.right_ghost_layer_size = rghosts
+ self.left_ghosts_layer_size = lghosts
+ self.ghosts_layer_size = np.maximum(lghosts, rghosts)
+
+ def __call__(self, tab, cdir, result):
+ stencil = self.stencil
+ svars = {stencil.dx:self._step[cdir]}
+
+ self.output_indices = [idx for idx in self.output_indices]
+ result[self.output_indices] = 0
+ for (offset,coeff) in self.stencil.iteritems(svars):
+ result[self.output_indices] += coeff*tab[self._indices(offset)]
+ return result
+
+ def _compute_indices(self, step):
+ pass
+
+ def compute_and_add(self):
+ pass
+ def compute_and_add_abs(self):
+ pass
+
+
class FDC2(FiniteDifference):
"""
@@ -225,7 +265,6 @@ class FDC2(FiniteDifference):
ghosts_layer_size = 1
def _compute_indices(self, step):
- print 'lol'
self._coeff = npw.asarray(1. / (2. * step))
self._m1 = []
@@ -509,3 +548,12 @@ class FDC4(FiniteDifference):
np.add(wk[self.output_indices], result[self.output_indices],
result[self.output_indices])
return result
+
+if __name__ == '__main__':
+ N = np.asarray((16,16,16))
+ A = np.random.rand(*N)
+ dx = 1.0/(N-1)
+ indices = np.ndindex(*N)
+
+ diff = FDkCp(2,2,step=dx,indices=indices)
+ diff(A,0,A)
diff --git a/hysop/numerics/stencil.py b/hysop/numerics/stencil.py
index c28ab1ec1..a0514b1be 100644
--- a/hysop/numerics/stencil.py
+++ b/hysop/numerics/stencil.py
@@ -114,7 +114,7 @@ class Stencil(object):
def _update_attributes(self):
self.dim = self.coeffs.ndim
self.shape = self.coeffs.shape
- self.L = self.origin
+ self.L = np.asarray([self.origin] if np.isscalar(self.origin) else self.origin)
self.R = self.shape-self.origin-1
def _delete_zeros(self, coeffs):
@@ -200,7 +200,7 @@ class Stencil(object):
raise RuntimeError('Stencil is not 2d !')
return sp.sparse.coo_matrix(self.coeffs,shape=self.shape,dtype=self.dtype)
- def iteritems(self):
+ def iteritems(self,svars={},include_factor=True):
"""
Return an (offset,coefficient) iterator iterating on all **non zero** coefficients.
Offset is taken from origin.
@@ -210,8 +210,15 @@ class Stencil(object):
it : itertools.iterator
Zipped offset and coefficient iterator.
"""
+ factor = self.factor if include_factor else 1
+ def mapfun(x):
+ offset = x-self.origin
+ value = factor*self.coeffs[x]
+ if isinstance(value, sm.Expr):
+ value = value.xreplace(svars)
+ return (offset,value)
iterator = np.ndindex(self.shape)
- iterator = it.imap(lambda x: (x-self.origin,self.coeffs[x]), iterator)
+ iterator = it.imap(mapfun, iterator)
iterator = it.ifilter(lambda x: x[1]!=0, iterator)
return iterator
--
GitLab