From eaf9f7d44ea5f346afcad1549a4eab1b90963f5b Mon Sep 17 00:00:00 2001 From: maroba Date: Mon, 9 Dec 2024 08:30:54 +0100 Subject: [PATCH 1/6] Moved old code to legacy subpackage --- findiff/__init__.py | 10 +-- findiff/legacy/__init__.py | 0 findiff/{ => legacy}/coefs.py | 0 findiff/{ => legacy}/diff.py | 120 +++++++++++++------------ findiff/{ => legacy}/grids.py | 0 findiff/{ => legacy}/operators.py | 0 findiff/{ => legacy}/pde.py | 0 findiff/{ => legacy}/stencils.py | 0 findiff/{ => legacy}/symbolic.py | 0 findiff/{ => legacy}/utils.py | 0 findiff/{ => legacy}/vector.py | 0 test/legacy/__init__.py | 0 test/{ => legacy}/test_bugs.py | 0 test/{ => legacy}/test_coefs.py | 70 ++++++++++----- test/{ => legacy}/test_diff.py | 47 +++++----- test/{ => legacy}/test_findiff.py | 137 +++++++++++++++-------------- test/{ => legacy}/test_grids.py | 4 +- test/{ => legacy}/test_pde.py | 93 ++++++++++++-------- test/{ => legacy}/test_scaling.py | 26 +++--- test/{ => legacy}/test_stencils.py | 80 +++++++++-------- test/{ => legacy}/test_symbolic.py | 4 +- test/{ => legacy}/test_utils.py | 8 +- test/{ => legacy}/test_vector.py | 73 ++++++++------- 23 files changed, 382 insertions(+), 290 deletions(-) create mode 100644 findiff/legacy/__init__.py rename findiff/{ => legacy}/coefs.py (100%) rename findiff/{ => legacy}/diff.py (79%) rename findiff/{ => legacy}/grids.py (100%) rename findiff/{ => legacy}/operators.py (100%) rename findiff/{ => legacy}/pde.py (100%) rename findiff/{ => legacy}/stencils.py (100%) rename findiff/{ => legacy}/symbolic.py (100%) rename findiff/{ => legacy}/utils.py (100%) rename findiff/{ => legacy}/vector.py (100%) create mode 100644 test/legacy/__init__.py rename test/{ => legacy}/test_bugs.py (100%) rename test/{ => legacy}/test_coefs.py (74%) rename test/{ => legacy}/test_diff.py (80%) rename test/{ => legacy}/test_findiff.py (77%) rename test/{ => legacy}/test_grids.py (81%) rename test/{ => legacy}/test_pde.py (76%) rename test/{ => legacy}/test_scaling.py (72%) rename test/{ => legacy}/test_stencils.py (76%) rename test/{ => legacy}/test_symbolic.py (95%) rename test/{ => legacy}/test_utils.py (89%) rename test/{ => legacy}/test_vector.py (68%) diff --git a/findiff/__init__.py b/findiff/__init__.py index 4ad4311..d15a5f5 100644 --- a/findiff/__init__.py +++ b/findiff/__init__.py @@ -19,8 +19,8 @@ __version__ = "0.10.2" -from .coefs import coefficients -from .operators import FinDiff, Coef, Identity, Coefficient -from .pde import PDE, BoundaryConditions -from .symbolic import SymbolicMesh, SymbolicDiff -from .vector import Gradient, Divergence, Curl, Laplacian +from .legacy.coefs import coefficients +from .legacy.operators import FinDiff, Coef, Identity, Coefficient +from .legacy.pde import PDE, BoundaryConditions +from .legacy.symbolic import SymbolicMesh, SymbolicDiff +from .legacy.vector import Gradient, Divergence, Curl, Laplacian diff --git a/findiff/legacy/__init__.py b/findiff/legacy/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/findiff/coefs.py b/findiff/legacy/coefs.py similarity index 100% rename from findiff/coefs.py rename to findiff/legacy/coefs.py diff --git a/findiff/diff.py b/findiff/legacy/diff.py similarity index 79% rename from findiff/diff.py rename to findiff/legacy/diff.py index eeacee4..79d5415 100644 --- a/findiff/diff.py +++ b/findiff/legacy/diff.py @@ -4,7 +4,7 @@ import scipy.sparse as sparse -from findiff.coefs import coefficients, coefficients_non_uni +from .coefs import coefficients, coefficients_non_uni from .stencils import StencilSet from .utils import * @@ -12,12 +12,13 @@ class Operator(object): - """ Base class for all operator classes """ + """Base class for all operator classes""" + pass class BinaryOperator(Operator): - """ Base class for all binary operators (like addition) that allow to combine or chain differential operators """ + """Base class for all binary operators (like addition) that allow to combine or chain differential operators""" def __init__(self, left, right): self.left = left @@ -34,7 +35,7 @@ def stencil(self, shape): class Plus(BinaryOperator): - """ A class to add two differential operators """ + """A class to add two differential operators""" def __init__(self, left, right): super().__init__(left, right) @@ -81,7 +82,7 @@ def matrix(self, shape, *args, **kwargs): class Minus(BinaryOperator): - """ A class to subtract two differential operators from each other """ + """A class to subtract two differential operators from each other""" def __init__(self, left, right): super().__init__(left, right) @@ -128,7 +129,7 @@ def matrix(self, shape, *args, **kwargs): class Mul(BinaryOperator): - """ A class to multiply (chain) two differential operators """ + """A class to multiply (chain) two differential operators""" def __init__(self, left, right): super().__init__(left, right) @@ -167,7 +168,7 @@ def apply(self, rhs, *args, **kwargs): return result def matrix(self, shape, *args, **kwargs): - """ Matrix representation of given operator product on an equidistant grid of given shape. + """Matrix representation of given operator product on an equidistant grid of given shape. :param shape: tuple with the shape of the grid :return: scipy sparse matrix representing the operator product @@ -182,7 +183,9 @@ def matrix(self, shape, *args, **kwargs): if isinstance(self.right, np.ndarray): right = sparse.diags(self.right.reshape(-1), 0) - elif isinstance(self.right, LinearMap) or isinstance(self.right, BinaryOperator): + elif isinstance(self.right, LinearMap) or isinstance( + self.right, BinaryOperator + ): right = self.right.matrix(shape, *args, **kwargs) else: right = self.right * sparse.diags(np.ones(shape).reshape(-1), 0) @@ -218,21 +221,21 @@ def __call__(self, rhs, *args, **kwargs): class Diff(LinearMap): - """ Representation of a single partial derivative based on finite differences. + """Representation of a single partial derivative based on finite differences. - This class is usually not used directly by the user, but is wrapped in - a FinDiff object. + This class is usually not used directly by the user, but is wrapped in + a FinDiff object. - :param axis: the numpy axis along which to apply the derivative + :param axis: the numpy axis along which to apply the derivative - :param order: the order of the derivative + :param order: the order of the derivative - :param kwargs: optional keyword arguments + :param kwargs: optional keyword arguments - Allowed keywords: + Allowed keywords: - `acc`: even integer - The desired accuracy order. + `acc`: even integer + The desired accuracy order. """ def __init__(self, axis, order, **kwargs): @@ -243,11 +246,11 @@ def __init__(self, axis, order, **kwargs): self.axis = axis self.order = order self.acc = None - if 'acc' in kwargs: - self.acc = kwargs['acc'] + if "acc" in kwargs: + self.acc = kwargs["acc"] def apply(self, u, *args, **kwargs): - """ Applies the partial derivative to a numpy array.""" + """Applies the partial derivative to a numpy array.""" h = None acc = DEFAULT_ACC @@ -260,15 +263,15 @@ def get_h(a): return h for key, value in kwargs.items(): - if key == 'h' or key == 'grid': + if key == "h" or key == "grid": h = get_h(value) break if h is None: h = get_h(args[0]) - if 'acc' in kwargs: - acc = kwargs['acc'] + if "acc" in kwargs: + acc = kwargs["acc"] if isinstance(h, np.ndarray): return self.diff_non_uni(u, h, **kwargs) @@ -278,10 +281,10 @@ def get_h(a): def diff(self, y, h, acc): """The core function to take a partial derivative on a uniform grid. - Central coefficients will be used whenever possible. Backward or forward - coefficients will be used if not enough points are available on either side, - i.e. forward coefficients for the low index boundary and backward coefficients - for the high index boundary. + Central coefficients will be used whenever possible. Backward or forward + coefficients will be used if not enough points are available on either side, + i.e. forward coefficients for the low index boundary and backward coefficients + for the high index boundary. """ dim = self.axis @@ -292,7 +295,8 @@ def diff(self, y, h, acc): npts = y.shape[dim] except AttributeError as err: raise ValueError( - "FinDiff objects can only be applied to arrays or evaluated(!) functions returning arrays") from err + "FinDiff objects can only be applied to arrays or evaluated(!) functions returning arrays" + ) from err scheme = "center" weights = coefs[scheme]["coefficients"] @@ -300,7 +304,9 @@ def diff(self, y, h, acc): num_bndry_points = len(weights) // 2 ref_slice = slice(num_bndry_points, npts - num_bndry_points, 1) - off_slices = [self._shift_slice(ref_slice, offsets[k], npts) for k in range(len(offsets))] + off_slices = [ + self._shift_slice(ref_slice, offsets[k], npts) for k in range(len(offsets)) + ] yd = np.zeros_like(y) @@ -311,7 +317,9 @@ def diff(self, y, h, acc): offsets = coefs[scheme]["offsets"] ref_slice = slice(0, num_bndry_points, 1) - off_slices = [self._shift_slice(ref_slice, offsets[k], npts) for k in range(len(offsets))] + off_slices = [ + self._shift_slice(ref_slice, offsets[k], npts) for k in range(len(offsets)) + ] self._apply_to_array(yd, y, weights, off_slices, ref_slice, dim) @@ -320,11 +328,13 @@ def diff(self, y, h, acc): offsets = coefs[scheme]["offsets"] ref_slice = slice(npts - num_bndry_points, npts, 1) - off_slices = [self._shift_slice(ref_slice, offsets[k], npts) for k in range(len(offsets))] + off_slices = [ + self._shift_slice(ref_slice, offsets[k], npts) for k in range(len(offsets)) + ] self._apply_to_array(yd, y, weights, off_slices, ref_slice, dim) - h_inv = 1. / h ** deriv + h_inv = 1.0 / h**deriv return yd * h_inv def diff_non_uni(self, y, coords, **kwargs): @@ -362,22 +372,24 @@ def diff_non_uni(self, y, coords, **kwargs): return yd - def matrix(self, shape, h=None, acc=None, coords=None, sparse_type=sparse.csr_matrix): - """ Matrix representation of the partial derivative. + def matrix( + self, shape, h=None, acc=None, coords=None, sparse_type=sparse.csr_matrix + ): + """Matrix representation of the partial derivative. - :param shape: Tuple with the shape of the grid (number of grid points in each dimension) + :param shape: Tuple with the shape of the grid (number of grid points in each dimension) - :param h: The grid spacing for the axis of the partial derivative - (only used for uniform grids) + :param h: The grid spacing for the axis of the partial derivative + (only used for uniform grids) - :param coords: The coordinate values of the grid on the axis of the partial derivative - (only used for non-uniform grids) + :param coords: The coordinate values of the grid on the axis of the partial derivative + (only used for non-uniform grids) - :param acc: The accuracy order of the derivative (even int) + :param acc: The accuracy order of the derivative (even int) - :param sparse_type: The scipy sparse matrix type used for the matrix representation. + :param sparse_type: The scipy sparse matrix type used for the matrix representation. - :returns matrix representation (scipy sparse matrix) + :returns matrix representation (scipy sparse matrix) """ if h is not None: @@ -385,7 +397,7 @@ def matrix(self, shape, h=None, acc=None, coords=None, sparse_type=sparse.csr_ma elif coords is not None: return sparse_type(self._matrix_nonuniform(shape, coords, acc)) else: - raise ValueError('Neither spacing nor coordinates given.') + raise ValueError("Neither spacing nor coordinates given.") def _matrix_nonuniform(self, shape, coords, acc): @@ -404,7 +416,7 @@ def _matrix_nonuniform(self, shape, coords, acc): for base_ind_long, base_ind_short in enumerate(short_inds): cd = coef_dicts[base_ind_short[self.axis]] - cs, os = cd['coefficients'], cd['offsets'] + cs, os = cd["coefficients"], cd["offsets"] for c, o in zip(cs, os): off_short = np.zeros(len(shape), dtype=int) off_short[self.axis] = int(o) @@ -430,10 +442,10 @@ def _matrix_uniform(self, shape, h=None, acc=None): mat = sparse.lil_matrix((siz, siz)) coeff_dict = coefficients(order, acc) - for scheme in ['center', 'forward', 'backward']: + for scheme in ["center", "forward", "backward"]: - offsets_1d = coeff_dict[scheme]['offsets'] - coeffs = coeff_dict[scheme]['coefficients'] + offsets_1d = coeff_dict[scheme]["offsets"] + coeffs = coeff_dict[scheme]["coefficients"] # translate offsets of given scheme to long format offsets_long = [] @@ -444,12 +456,12 @@ def _matrix_uniform(self, shape, h=None, acc=None): offsets_long.append(o_long) # determine points where to evaluate current scheme in long format - nside = len(coeff_dict['center']['coefficients']) // 2 - if scheme == 'center': + nside = len(coeff_dict["center"]["coefficients"]) // 2 + if scheme == "center": multi_slice = [slice(None, None)] * ndims multi_slice[axis] = slice(nside, -nside) Is = long_indices_nd[tuple(multi_slice)].reshape(-1) - elif scheme == 'forward': + elif scheme == "forward": multi_slice = [slice(None, None)] * ndims multi_slice[axis] = slice(0, nside) Is = long_indices_nd[tuple(multi_slice)].reshape(-1) @@ -459,7 +471,7 @@ def _matrix_uniform(self, shape, h=None, acc=None): Is = long_indices_nd[tuple(multi_slice)].reshape(-1) for o, c in zip(offsets_long, coeffs): - v = c / h ** order + v = c / h**order mat[Is, Is + o] = v return mat @@ -489,7 +501,7 @@ def _apply_to_array(self, yd, y, weights, off_slices, ref_slice, dim): for w, s in zip(weights, off_slices): off_multi_slice = [all] * ndims off_multi_slice[dim] = s - if abs(1 - w) < 1.E-14: + if abs(1 - w) < 1.0e-14: yd[tuple(ref_multi_slice)] += y[tuple(off_multi_slice)] else: yd[tuple(ref_multi_slice)] += w * y[tuple(off_multi_slice)] @@ -503,7 +515,7 @@ def _shift_slice(self, sl, off, max_index): class Id(LinearMap): - """ The identity operator. When applied to an array, returns the same array (not a copy) """ + """The identity operator. When applied to an array, returns the same array (not a copy)""" def __init__(self): self.value = 1 @@ -512,7 +524,7 @@ def apply(self, rhs, *args, **kwargs): return rhs def matrix(self, shape): - """ Matrix representation of the identity operator, i.e. identity matrix of given shape. + """Matrix representation of the identity operator, i.e. identity matrix of given shape. :param shape: Shape of the arrays to which Id shall be applied :type shape: tuple of ints diff --git a/findiff/grids.py b/findiff/legacy/grids.py similarity index 100% rename from findiff/grids.py rename to findiff/legacy/grids.py diff --git a/findiff/operators.py b/findiff/legacy/operators.py similarity index 100% rename from findiff/operators.py rename to findiff/legacy/operators.py diff --git a/findiff/pde.py b/findiff/legacy/pde.py similarity index 100% rename from findiff/pde.py rename to findiff/legacy/pde.py diff --git a/findiff/stencils.py b/findiff/legacy/stencils.py similarity index 100% rename from findiff/stencils.py rename to findiff/legacy/stencils.py diff --git a/findiff/symbolic.py b/findiff/legacy/symbolic.py similarity index 100% rename from findiff/symbolic.py rename to findiff/legacy/symbolic.py diff --git a/findiff/utils.py b/findiff/legacy/utils.py similarity index 100% rename from findiff/utils.py rename to findiff/legacy/utils.py diff --git a/findiff/vector.py b/findiff/legacy/vector.py similarity index 100% rename from findiff/vector.py rename to findiff/legacy/vector.py diff --git a/test/legacy/__init__.py b/test/legacy/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/test/test_bugs.py b/test/legacy/test_bugs.py similarity index 100% rename from test/test_bugs.py rename to test/legacy/test_bugs.py diff --git a/test/test_coefs.py b/test/legacy/test_coefs.py similarity index 74% rename from test/test_coefs.py rename to test/legacy/test_coefs.py index 9054dd6..667a907 100644 --- a/test/test_coefs.py +++ b/test/legacy/test_coefs.py @@ -1,10 +1,11 @@ import sys -sys.path.insert(1, '..') + +sys.path.insert(1, "..") import unittest from findiff import coefficients -from findiff.coefs import coefficients_non_uni -from findiff.coefs import calc_coefs +from findiff.legacy.coefs import coefficients_non_uni +from findiff.legacy.coefs import calc_coefs import numpy as np from sympy import Rational @@ -19,7 +20,7 @@ def test_order2_acc2(self): with self.subTest(): coefs = c["center"]["coefficients"] - np.testing.assert_array_almost_equal(np.array([1., -2., 1.]), coefs) + np.testing.assert_array_almost_equal(np.array([1.0, -2.0, 1.0]), coefs) offs = c["center"]["offsets"] np.testing.assert_array_almost_equal(np.array([-1, 0, 1]), offs) @@ -29,7 +30,9 @@ def test_order2_acc2(self): np.testing.assert_array_almost_equal(np.array([0, 1, 2, 3]), offs) coefs = c["backward"]["coefficients"] - np.testing.assert_array_almost_equal(np.array(([2, -5, 4, -1])[::-1]), coefs) + np.testing.assert_array_almost_equal( + np.array(([2, -5, 4, -1])[::-1]), coefs + ) offs = c["backward"]["offsets"] np.testing.assert_array_almost_equal(np.array([-3, -2, -1, 0]), offs) @@ -50,7 +53,9 @@ def test_order1_acc2(self): np.testing.assert_array_almost_equal(np.array([0, 1, 2]), offs) coefs = c["backward"]["coefficients"] - np.testing.assert_array_almost_equal(-np.array([-1.5, 2, -0.5])[::-1], coefs) + np.testing.assert_array_almost_equal( + -np.array([-1.5, 2, -0.5])[::-1], coefs + ) offs = c["backward"]["offsets"] np.testing.assert_array_almost_equal(np.array([-2, -1, 0]), offs) @@ -60,19 +65,27 @@ def test_order1_acc4(self): c = coefficients(deriv=1, acc=4, analytic_inv=analytic_inv) with self.subTest(): coefs = c["center"]["coefficients"] - np.testing.assert_array_almost_equal(np.array([1/12, -2/3, 0, 2/3, -1/12]), coefs) + np.testing.assert_array_almost_equal( + np.array([1 / 12, -2 / 3, 0, 2 / 3, -1 / 12]), coefs + ) offs = c["center"]["offsets"] np.testing.assert_array_almost_equal(np.array([-2, -1, 0, 1, 2]), offs) coefs = c["forward"]["coefficients"] - np.testing.assert_array_almost_equal(np.array([-25/12, 4, -3, 4/3, -1/4]), coefs) + np.testing.assert_array_almost_equal( + np.array([-25 / 12, 4, -3, 4 / 3, -1 / 4]), coefs + ) offs = c["forward"]["offsets"] np.testing.assert_array_almost_equal(np.array([0, 1, 2, 3, 4]), offs) coefs = c["backward"]["coefficients"] - np.testing.assert_array_almost_equal(-np.array([-25/12, 4, -3, 4/3, -1/4])[::-1], coefs) + np.testing.assert_array_almost_equal( + -np.array([-25 / 12, 4, -3, 4 / 3, -1 / 4])[::-1], coefs + ) offs = c["backward"]["offsets"] - np.testing.assert_array_almost_equal(-np.array([0, 1, 2, 3, 4])[::-1], offs) + np.testing.assert_array_almost_equal( + -np.array([0, 1, 2, 3, 4])[::-1], offs + ) def test_order2_acc4(self): @@ -80,19 +93,28 @@ def test_order2_acc4(self): c = coefficients(deriv=2, acc=4, analytic_inv=analytic_inv) with self.subTest(): coefs = c["center"]["coefficients"] - np.testing.assert_array_almost_equal(np.array([-1/12, 4/3, -2.5, 4/3, -1/12]), coefs) + np.testing.assert_array_almost_equal( + np.array([-1 / 12, 4 / 3, -2.5, 4 / 3, -1 / 12]), coefs + ) offs = c["center"]["offsets"] np.testing.assert_array_almost_equal(np.array([-2, -1, 0, 1, 2]), offs) coefs = c["forward"]["coefficients"] - np.testing.assert_array_almost_equal(np.array([15/4, -77/6, 107/6, -13, 61/12, -5/6]), coefs) + np.testing.assert_array_almost_equal( + np.array([15 / 4, -77 / 6, 107 / 6, -13, 61 / 12, -5 / 6]), coefs + ) offs = c["forward"]["offsets"] np.testing.assert_array_almost_equal(np.array([0, 1, 2, 3, 4, 5]), offs) coefs = c["backward"]["coefficients"] - np.testing.assert_array_almost_equal(np.array([15/4, -77/6, 107/6, -13, 61/12, -5/6])[::-1], coefs) + np.testing.assert_array_almost_equal( + np.array([15 / 4, -77 / 6, 107 / 6, -13, 61 / 12, -5 / 6])[::-1], + coefs, + ) offs = c["backward"]["offsets"] - np.testing.assert_array_almost_equal(-np.array([0, 1, 2, 3, 4, 5])[::-1], offs) + np.testing.assert_array_almost_equal( + -np.array([0, 1, 2, 3, 4, 5])[::-1], offs + ) def test_calc_accuracy_central_deriv2_acc2(self): for analytic_inv in [True, False]: @@ -120,13 +142,17 @@ def test_calc_accuracy_left0_right3_deriv1_acc3(self): def test_calc_accuracy_from_offsets_symbolic(self): for analytic_inv in [True, False]: - coefs = calc_coefs(2, [0, 1, 2, 3], symbolic=True, analytic_inv=analytic_inv) + coefs = calc_coefs( + 2, [0, 1, 2, 3], symbolic=True, analytic_inv=analytic_inv + ) with self.subTest(): self.assertEqual(2, coefs["accuracy"]) def test_calc_accuracy_from_offsets_symbolic(self): for analytic_inv in [True, False]: - coefs = calc_coefs(2, [-4, -2, 0, 2, 4], symbolic=True, analytic_inv=analytic_inv) + coefs = calc_coefs( + 2, [-4, -2, 0, 2, 4], symbolic=True, analytic_inv=analytic_inv + ) with self.subTest(): self.assertEqual(4, coefs["accuracy"]) @@ -134,13 +160,17 @@ def test_calc_coefs_from_offsets(self): for analytic_inv in [True, False]: coefs = calc_coefs(1, [-2, 0, 1], analytic_inv=analytic_inv) with self.subTest(): - np.testing.assert_array_almost_equal(coefs["coefficients"], [-1./6, -0.5, 2./3]) + np.testing.assert_array_almost_equal( + coefs["coefficients"], [-1.0 / 6, -0.5, 2.0 / 3] + ) def test_calc_coefs_from_offsets_no_central_point(self): for analytic_inv in [True, False]: coefs = calc_coefs(1, [-2, 1, 2], analytic_inv=analytic_inv) with self.subTest(): - np.testing.assert_array_almost_equal(coefs["coefficients"], [-1./4, 0, 1./4]) + np.testing.assert_array_almost_equal( + coefs["coefficients"], [-1.0 / 4, 0, 1.0 / 4] + ) def test_calc_coefs_symbolic(self): for analytic_inv in [True, False]: @@ -156,7 +186,7 @@ def test_non_uniform(self): for analytic_inv in [True, False]: c_uni = coefficients(deriv=2, acc=2, analytic_inv=analytic_inv) - coefs_uni = c_uni["center"]["coefficients"]/dx**2 + coefs_uni = c_uni["center"]["coefficients"] / dx**2 c_non_uni = coefficients_non_uni(deriv=2, acc=2, coords=x, idx=5) coefs_non_uni = c_non_uni["coefficients"] @@ -181,5 +211,5 @@ def test_invalid_deriv_raises_exception(self): coefficients_non_uni(-1, 2, None, None) -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/test/test_diff.py b/test/legacy/test_diff.py similarity index 80% rename from test/test_diff.py rename to test/legacy/test_diff.py index 7e6ac08..ffcd105 100644 --- a/test/test_diff.py +++ b/test/legacy/test_diff.py @@ -1,20 +1,19 @@ import sys -sys.path.insert(1, '..') + +sys.path.insert(1, "..") import unittest -import numpy as np -from numpy.testing import assert_array_almost_equal, assert_array_equal -from findiff.diff import * -#import matplotlib.pyplot as plt -import findiff +from numpy.testing import assert_array_almost_equal +from findiff.legacy.diff import * + +# import matplotlib.pyplot as plt from numpy import cos, sin -from findiff.grids import UniformGrid class DiffTest(unittest.TestCase): def test_partial_d_dx(self): - shape = 101, + shape = (101,) x, dx = make_grid(shape, (0, 1)) u = x**2 @@ -26,10 +25,10 @@ def test_partial_d_dx(self): assert_array_almost_equal(expected, actual) def test_partial_d2_dx2(self): - shape = 101, + shape = (101,) x, dx = make_grid(shape, (0, 1)) - u = x ** 2 + u = x**2 expected = 2 fd = Diff(0, 2) @@ -38,15 +37,17 @@ def test_partial_d2_dx2(self): assert_array_almost_equal(expected, actual) def test_partial_d_dx_acc(self): - shape = 11, + shape = (11,) x, dx = make_grid(shape, (0, 1)) - u = x ** 3 - expected = 3*x**2 + u = x**3 + expected = 3 * x**2 fd = Diff(0, 1) actual = fd(u, dx) - np.testing.assert_raises(AssertionError, assert_array_almost_equal, expected, actual) + np.testing.assert_raises( + AssertionError, assert_array_almost_equal, expected, actual + ) actual = fd(u, dx, acc=4) assert_array_almost_equal(expected, actual) @@ -124,7 +125,7 @@ def test_multiply_with_variable_from_left(self): fd = Coef(3 * X * Y) * Diff(0, 2) actual = fd(u, dy, acc=4) - expected = - 3*X *Y* np.sin(X) * np.sin(Y) + expected = -3 * X * Y * np.sin(X) * np.sin(Y) assert_array_almost_equal(expected, actual) def test_identity(self): @@ -144,10 +145,10 @@ def test_linear_comb(self): (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) u = np.sin(X) * np.sin(Y) - fd = Coef(3*X)*Diff(0, 2) - Coef(Y) * (Diff(1, 2) + Diff(0, 1)) - actual = fd(u, h={0: dx, 1:dy}, acc=4) + fd = Coef(3 * X) * Diff(0, 2) - Coef(Y) * (Diff(1, 2) + Diff(0, 1)) + actual = fd(u, h={0: dx, 1: dy}, acc=4) - expected = -3*X*u - Y * (-sin(X)*sin(Y) + cos(X)*sin(Y)) + expected = -3 * X * u - Y * (-sin(X) * sin(Y) + cos(X) * sin(Y)) assert_array_almost_equal(expected, actual) def test_apply_base_binary_operator_raises_exception(self): @@ -161,10 +162,12 @@ def make_grid(shape, edges): if len(shape) == 1: x = np.linspace(*edges, shape[0]) - dx = x[1]-x[0] + dx = x[1] - x[0] return x, dx - axes = tuple([np.linspace(edges[k][0], edges[k][1], shape[k]) for k in range(len(shape))]) - coords = np.meshgrid(*axes, indexing='ij') - spacings = [axes[k][1]-axes[k][0] for k in range(len(shape))] + axes = tuple( + [np.linspace(edges[k][0], edges[k][1], shape[k]) for k in range(len(shape))] + ) + coords = np.meshgrid(*axes, indexing="ij") + spacings = [axes[k][1] - axes[k][0] for k in range(len(shape))] return axes, spacings, coords diff --git a/test/test_findiff.py b/test/legacy/test_findiff.py similarity index 77% rename from test/test_findiff.py rename to test/legacy/test_findiff.py index 529cfb0..5bf51c3 100644 --- a/test/test_findiff.py +++ b/test/legacy/test_findiff.py @@ -1,10 +1,11 @@ import sys -sys.path.insert(1, '..') + +sys.path.insert(1, "..") import unittest import numpy as np from numpy.testing import assert_array_almost_equal, assert_array_equal -from findiff.operators import FinDiff, Coef, Identity +from findiff.legacy.operators import FinDiff, Coef, Identity class FinDiffTest(unittest.TestCase): @@ -13,7 +14,7 @@ def test_partial_diff_1d(self): nx = 11 x = np.linspace(0, 1, nx) u = x**3 - ux_ex = 3*x**2 + ux_ex = 3 * x**2 fd = FinDiff(0, x[1] - x[0], 1) @@ -21,7 +22,6 @@ def test_partial_diff_1d(self): assert_array_almost_equal(ux, ux_ex, decimal=5) - def test_partial_diff(self): nx = 100 x = np.linspace(0, np.pi, nx) @@ -36,7 +36,7 @@ def test_partial_diff(self): ny = 100 y = np.linspace(0, np.pi, ny) - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") u = np.sin(X) * np.sin(Y) uxy_ex = np.cos(X) * np.cos(Y) @@ -58,7 +58,7 @@ def test_plus(self): d = d_dx + d_dy u1 = d(u) - u1_ex = 2*X + 2*Y + u1_ex = 2 * X + 2 * Y assert_array_almost_equal(u1, u1_ex) @@ -72,7 +72,7 @@ def test_multiply(self): d = Coef(X) * d2_dx2 u1 = d(u) - assert_array_almost_equal(u1, 2*X) + assert_array_almost_equal(u1, 2 * X) def test_multiply_operators(self): @@ -85,7 +85,7 @@ def test_multiply_operators(self): uxx = d2_dx2_test(u) - assert_array_almost_equal(uxx, np.ones_like(X)*2) + assert_array_almost_equal(uxx, np.ones_like(X) * 2) def test_laplace(self): @@ -98,7 +98,7 @@ def test_laplace(self): laplace = d2_dx2 + d2_dy2 + d2_dz2 lap_u = laplace(u) - assert_array_almost_equal(lap_u, 6*X + 6*Y + 6*Z) + assert_array_almost_equal(lap_u, 6 * X + 6 * Y + 6 * Z) d_dx, d_dy, d_dz = [FinDiff(i, h[i]) for i in range(3)] @@ -130,7 +130,7 @@ def test_identity_2d(self): (X, Y), (x, y), _ = grid(2, 100, -1, 1) - u = X ** 2 + Y ** 2 + u = X**2 + Y**2 identity = Identity() assert_array_equal(u, identity(u)) @@ -144,7 +144,7 @@ def test_identity_2d(self): dx = x[1] - x[0] d_dx = FinDiff(0, dx) linop = d_dx + 2 * Identity() - assert_array_almost_equal(2 * X + 2*u, linop(u)) + assert_array_almost_equal(2 * X + 2 * u, linop(u)) def test_spac(self): @@ -155,13 +155,13 @@ def test_spac(self): d_dx = FinDiff(0, dx) d_dy = FinDiff(1, dy) - assert_array_almost_equal(2*X, d_dx(u)) - assert_array_almost_equal(2*Y, d_dy(u)) + assert_array_almost_equal(2 * X, d_dx(u)) + assert_array_almost_equal(2 * Y, d_dy(u)) d_dx = FinDiff(0, dx) d_dy = FinDiff(1, dy) - u = X*Y + u = X * Y d2_dxdy = d_dx * d_dy assert_array_almost_equal(np.ones_like(u), d2_dxdy(u)) @@ -174,14 +174,14 @@ def test_mixed_partials(self): d3_dxdydz = FinDiff((0, dx), (1, dy), (2, dz)) diffed = d3_dxdydz(u) - assert_array_almost_equal(8*X*Y*Z, diffed) + assert_array_almost_equal(8 * X * Y * Z, diffed) def test_linear_combinations(self): (X, Y, Z), _, (dx, dy, dz) = grid(3, 30, 0, 1) - u = X ** 2 + Y ** 2 + Z ** 2 + u = X**2 + Y**2 + Z**2 d = Coef(X) * FinDiff(0, dx) + Coef(Y**2) * FinDiff(1, dy, 2) - assert_array_almost_equal(d(u), 2*X**2 + 2*Y**2) + assert_array_almost_equal(d(u), 2 * X**2 + 2 * Y**2) def test_minus(self): @@ -195,7 +195,7 @@ def test_minus(self): d = d_dx - d_dy u1 = d(u) - u1_ex = 2*X - 2*Y + u1_ex = 2 * X - 2 * Y assert_array_almost_equal(u1, u1_ex) @@ -224,7 +224,7 @@ def test_local_stencil_single_axis_center_2d_compared_with_findiff(self): n = 70 (X, Y), _, (dx, dy) = grid(2, n, -1, 1) - u = X ** 3 * Y ** 3 + u = X**3 * Y**3 d4_dx2dy2 = FinDiff(1, dx, 2) expected = d4_dx2dy2(u) @@ -253,7 +253,7 @@ def test_local_stencil_operator_addition(self): def test_local_stencil_operator_mixed_partials(self): x = np.linspace(0, 10, 101) y = np.linspace(0, 10, 101) - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") u = X * Y dx = x[1] - x[0] dy = y[1] - y[0] @@ -266,7 +266,7 @@ def test_local_stencil_operator_mixed_partials(self): def test_local_stencil_operator_multiplication(self): x = np.linspace(0, 10, 101) y = np.linspace(0, 10, 101) - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") u = X * Y dx = x[1] - x[0] dy = y[1] - y[0] @@ -279,7 +279,7 @@ def test_local_stencil_operator_multiplication(self): def test_local_stencil_operator_with_coef(self): x = np.linspace(0, 10, 101) y = np.linspace(0, 10, 101) - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") u = X * Y dx = x[1] - x[0] dy = y[1] - y[0] @@ -287,7 +287,7 @@ def test_local_stencil_operator_with_coef(self): stencil1 = d1x.stencil(u.shape) du_dx = stencil1.apply_all(u) - np.testing.assert_array_almost_equal(2*np.ones_like(X), du_dx) + np.testing.assert_array_almost_equal(2 * np.ones_like(X), du_dx) def dict_almost_equal(self, d1, d2): @@ -299,12 +299,14 @@ def dict_almost_equal(self, d1, d2): def test_matrix_1d(self): x = np.linspace(0, 6, 7) - d2_dx2 = FinDiff(0, x[1]-x[0], 2) + d2_dx2 = FinDiff(0, x[1] - x[0], 2) u = x**2 mat = d2_dx2.matrix(u.shape) - np.testing.assert_array_almost_equal(2*np.ones_like(x), mat.dot(u.reshape(-1))) + np.testing.assert_array_almost_equal( + 2 * np.ones_like(x), mat.dot(u.reshape(-1)) + ) def test_matrix_2d(self): thr = np.get_printoptions()["threshold"] @@ -312,21 +314,23 @@ def test_matrix_2d(self): np.set_printoptions(threshold=np.inf) np.set_printoptions(linewidth=500) x, y = [np.linspace(0, 4, 5)] * 2 - X, Y = np.meshgrid(x, y, indexing='ij') - laplace = FinDiff(0, x[1]-x[0], 2) + FinDiff(0, y[1]-y[0], 2) - #d = FinDiff(1, y[1]-y[0], 2) + X, Y = np.meshgrid(x, y, indexing="ij") + laplace = FinDiff(0, x[1] - x[0], 2) + FinDiff(0, y[1] - y[0], 2) + # d = FinDiff(1, y[1]-y[0], 2) u = X**2 + Y**2 mat = laplace.matrix(u.shape) - np.testing.assert_array_almost_equal(4 * np.ones_like(X).reshape(-1), mat.dot(u.reshape(-1))) + np.testing.assert_array_almost_equal( + 4 * np.ones_like(X).reshape(-1), mat.dot(u.reshape(-1)) + ) np.set_printoptions(threshold=thr) np.set_printoptions(linewidth=lw) def test_matrix_2d_mixed(self): x, y = [np.linspace(0, 5, 6), np.linspace(0, 6, 7)] - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") d2_dxdy = FinDiff((0, x[1] - x[0]), (1, y[1] - y[0])) u = X**2 * Y**2 @@ -337,7 +341,7 @@ def test_matrix_2d_mixed(self): np.testing.assert_array_almost_equal(expected, actual) def test_matrix_1d_coeffs(self): - shape = 11, + shape = (11,) x = np.linspace(0, 10, 11) dx = x[1] - x[0] @@ -352,19 +356,19 @@ def test_matrix_1d_coeffs(self): @unittest.skip def test_eigs(self): - shape = 8, + shape = (8,) x = np.linspace(-10, 10, shape[0]) dx = x[1] - x[0] hbar = 1 m = 1 omega = 1 - T = Coef(-hbar**2/(2*m)) * FinDiff(0, dx, 2) - V = Coef(0.5*omega*x**2) * Identity() + T = Coef(-(hbar**2) / (2 * m)) * FinDiff(0, dx, 2) + V = Coef(0.5 * omega * x**2) * Identity() H = T + V - print("\n",T.matrix(shape).toarray()) - print("\n",V.matrix(shape).toarray()) - print("\n",H.matrix(shape).toarray()) + print("\n", T.matrix(shape).toarray()) + print("\n", V.matrix(shape).toarray()) + print("\n", H.matrix(shape).toarray()) print(H.matrix(shape).toarray()) vals, vecs = H.eigs(shape) @@ -376,30 +380,29 @@ class TestFinDiffNonUniform(unittest.TestCase): def test_1d_different_accs(self): x = np.r_[np.arange(0, 4, 0.05), np.arange(4, 10, 1)] - f = np.exp(-x**2) + f = np.exp(-(x**2)) d_dx = FinDiff(0, x, acc=4) f_x = d_dx(f) - assert_array_almost_equal(- 2*x * np.exp(-x**2), f_x, decimal=4) + assert_array_almost_equal(-2 * x * np.exp(-(x**2)), f_x, decimal=4) # same, but this time with default acc x = np.linspace(0, 1, 100) - f = np.exp(-x ** 2) + f = np.exp(-(x**2)) d_dx = FinDiff(0, x) f_x = d_dx(f) - assert_array_almost_equal(- 2*x * np.exp(-x**2), f_x, decimal=4) - + assert_array_almost_equal(-2 * x * np.exp(-(x**2)), f_x, decimal=4) def test_3d_different_accs(self): x = np.r_[np.arange(0, 4, 0.05), np.arange(4, 10, 1)] y = np.r_[np.arange(0, 4, 0.05), np.arange(4, 10, 1)] z = np.r_[np.arange(0, 4.5, 0.05), np.arange(4.5, 10, 1)] - X, Y, Z = np.meshgrid(x, y, z, indexing='ij') - f = np.exp(-X**2-Y**2-Z**2) + X, Y, Z = np.meshgrid(x, y, z, indexing="ij") + f = np.exp(-(X**2) - Y**2 - Z**2) d_dy = FinDiff(1, y, acc=4) fy = d_dy(f) - fye = - 2 * Y * np.exp(-X**2-Y**2-Z**2) + fye = -2 * Y * np.exp(-(X**2) - Y**2 - Z**2) assert_array_almost_equal(fy, fye, decimal=4) d_dy = FinDiff(1, y, acc=6) @@ -411,28 +414,28 @@ def test_1d_linear_combination(self): d1 = Coef(x) * FinDiff(0, x, acc=4) d2 = Coef(x**2) * FinDiff(0, x, 2, acc=4) - assert_array_almost_equal(3*x**3, d1(x**3)) - assert_array_almost_equal(6*x**3, d2(x**3)) - assert_array_almost_equal(9*x**3, (d1 + d2)(x**3)) + assert_array_almost_equal(3 * x**3, d1(x**3)) + assert_array_almost_equal(6 * x**3, d2(x**3)) + assert_array_almost_equal(9 * x**3, (d1 + d2)(x**3)) def test_3d_linear_combination(self): x, y, z = [np.sqrt(np.linspace(0, 1, 10))] * 3 - X, Y, Z = np.meshgrid(x, y, z, indexing='ij') + X, Y, Z = np.meshgrid(x, y, z, indexing="ij") d1 = Coef(X) * FinDiff(0, x, acc=4) d2 = Coef(Y) * FinDiff(1, y, acc=4) - assert_array_almost_equal(3 * X ** 3, d1(X**3 + Y**3 + Z**3)) - assert_array_almost_equal(3 * Y ** 3, d2(X**3 + Y**3 + Z**3)) + assert_array_almost_equal(3 * X**3, d1(X**3 + Y**3 + Z**3)) + assert_array_almost_equal(3 * Y**3, d2(X**3 + Y**3 + Z**3)) assert_array_almost_equal(3 * (X**3 + Y**3), (d1 + d2)(X**3 + Y**3 + Z**3)) def test_several_tuples_as_args(self): x, y, z = [np.sqrt(np.linspace(0, 1, 10))] * 3 - X, Y, Z = np.meshgrid(x, y, z, indexing='ij') + X, Y, Z = np.meshgrid(x, y, z, indexing="ij") d1 = FinDiff((0, x), (1, y), acc=4) - assert_array_almost_equal(d1(X*Y), np.ones_like(X)) + assert_array_almost_equal(d1(X * Y), np.ones_like(X)) def test_accepts_not_more_than_three_args(self): x = np.linspace(0, 1, 10) @@ -446,7 +449,9 @@ def test_matrix_1d_nonuni(self): mat = d2_dx2.matrix(u.shape) - np.testing.assert_array_almost_equal(2*np.ones_like(x), mat.dot(u.reshape(-1))) + np.testing.assert_array_almost_equal( + 2 * np.ones_like(x), mat.dot(u.reshape(-1)) + ) def test_matrix_2d_nonuni(self): thr = np.get_printoptions()["threshold"] @@ -454,21 +459,23 @@ def test_matrix_2d_nonuni(self): np.set_printoptions(threshold=np.inf) np.set_printoptions(linewidth=500) x, y = [np.linspace(0, 4, 5)] * 2 - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") laplace = FinDiff(0, x, 2) + FinDiff(1, y, 2) - #d = FinDiff(1, y[1]-y[0], 2) + # d = FinDiff(1, y[1]-y[0], 2) u = X**2 + Y**2 mat = laplace.matrix(u.shape) - np.testing.assert_array_almost_equal(4 * np.ones_like(X).reshape(-1), mat.dot(u.reshape(-1))) + np.testing.assert_array_almost_equal( + 4 * np.ones_like(X).reshape(-1), mat.dot(u.reshape(-1)) + ) np.set_printoptions(threshold=thr) np.set_printoptions(linewidth=lw) def test_matrix_2d_mixed_nonuni(self): x, y = [np.linspace(0, 5, 6), np.linspace(0, 6, 7)] - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") d2_dxdy = FinDiff((0, x), (1, y)) u = X**2 * Y**2 @@ -479,7 +486,7 @@ def test_matrix_2d_mixed_nonuni(self): np.testing.assert_array_almost_equal(expected, actual) def test_matrix_1d_coeffs_nonuni(self): - shape = 11, + shape = (11,) x = np.linspace(0, 10, 11) L = Coef(x) * FinDiff(0, x, 2) @@ -494,13 +501,15 @@ def test_matrix_1d_coeffs_nonuni(self): def test_matrix_3d_nonuni_performance(self): x = y = z = np.linspace(0, 4, 30) - X, Y, Z = np.meshgrid(x, y, z, indexing='ij') + X, Y, Z = np.meshgrid(x, y, z, indexing="ij") laplace = FinDiff(0, x, 2) + FinDiff(1, y, 2) + FinDiff(2, z, 2) u = X**2 + Y**2 + Z**2 mat = laplace.matrix(u.shape) - np.testing.assert_array_almost_equal(6 * np.ones_like(X).reshape(-1), mat.dot(u.reshape(-1))) + np.testing.assert_array_almost_equal( + 6 * np.ones_like(X).reshape(-1), mat.dot(u.reshape(-1)) + ) def grid(ndim, npts, a, b): @@ -513,11 +522,11 @@ def grid(ndim, npts, a, b): npts = [npts] * ndim coords = [np.linspace(a[i], b[i], npts[i]) for i in range(ndim)] - mesh = np.meshgrid(*coords, indexing='ij') + mesh = np.meshgrid(*coords, indexing="ij") spac = [coords[i][1] - coords[i][0] for i in range(ndim)] return mesh, coords, spac -if __name__ == '__main__': - unittest.main() \ No newline at end of file +if __name__ == "__main__": + unittest.main() diff --git a/test/test_grids.py b/test/legacy/test_grids.py similarity index 81% rename from test/test_grids.py rename to test/legacy/test_grids.py index b613534..de07530 100644 --- a/test/test_grids.py +++ b/test/legacy/test_grids.py @@ -2,7 +2,7 @@ import numpy as np -from findiff.grids import UniformGrid +from findiff.legacy.grids import UniformGrid class TestUniformGrid(unittest.TestCase): @@ -18,4 +18,4 @@ def test_init_2d(self): def test_init_1d(self): grid = UniformGrid(30, 0.1) self.assertEqual(0.1, grid.spacing(0)) - np.testing.assert_array_equal((0,), grid.center) \ No newline at end of file + np.testing.assert_array_equal((0,), grid.center) diff --git a/test/test_pde.py b/test/legacy/test_pde.py similarity index 76% rename from test/test_pde.py rename to test/legacy/test_pde.py index 762b0ef..0f8392d 100644 --- a/test/test_pde.py +++ b/test/legacy/test_pde.py @@ -1,13 +1,14 @@ import sys -sys.path.insert(1, '..') + +sys.path.insert(1, "..") import unittest -import numpy as np -from numpy.testing import assert_array_almost_equal, assert_array_equal -from findiff.operators import FinDiff, Identity, Coef -from findiff.pde import * -#import matplotlib.pyplot as plt -#from mpl_toolkits import mplot3d +from findiff.legacy.operators import FinDiff, Identity, Coef +from findiff.legacy.pde import * + +# import matplotlib.pyplot as plt +# from mpl_toolkits import mplot3d + class TestPDE(unittest.TestCase): @@ -43,7 +44,7 @@ def test_1d_dirichlet_inhom(self): bc[0] = 1 bc[-1] = 2 - pde = PDE(L, 6*x, bc) + pde = PDE(L, 6 * x, bc) u = pde.solve() expected = x**3 + 1 @@ -65,7 +66,7 @@ def test_1d_neumann_hom(self): pde = PDE(L, np.zeros_like(x), bc) u = pde.solve() - expected = 2*x + 1 + expected = 2 * x + 1 np.testing.assert_array_almost_equal(expected, u) def test_2d_dirichlet_hom(self): @@ -74,7 +75,7 @@ def test_2d_dirichlet_hom(self): x, y = np.linspace(0, 1, shape[0]), np.linspace(0, 1, shape[1]) dx, dy = x[1] - x[0], y[1] - y[0] - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") L = FinDiff(0, dx, 2) + FinDiff(1, dy, 2) expected = X + 1 @@ -96,11 +97,11 @@ def test_2d_dirichlet_inhom(self): x, y = np.linspace(0, 1, shape[0]), np.linspace(0, 1, shape[1]) dx, dy = x[1] - x[0], y[1] - y[0] - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") L = FinDiff(0, dx, 2) + FinDiff(1, dy, 2) expected = X**3 + Y**3 + 1 - f = 6*X + 6*Y + f = 6 * X + 6 * Y bc = BoundaryConditions(shape) @@ -118,7 +119,7 @@ def test_2d_neumann_hom(self): x, y = np.linspace(0, 1, shape[0]), np.linspace(0, 1, shape[1]) dx, dy = x[1] - x[0], y[1] - y[0] - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") L = FinDiff(0, dx, 2) + FinDiff(1, dy, 2) expected = X**2 + Y**2 + 1 @@ -130,8 +131,8 @@ def test_2d_neumann_hom(self): bc[0, :] = expected bc[-1, :] = expected - bc[:, 0] = d_dy, 2*Y - bc[:, -1] = d_dy, 2*Y + bc[:, 0] = d_dy, 2 * Y + bc[:, -1] = d_dy, 2 * Y pde = PDE(L, f, bc) u = pde.solve() @@ -140,7 +141,7 @@ def test_2d_neumann_hom(self): def test_1d_oscillator_free_dirichlet(self): n = 300 - shape = n, + shape = (n,) t = np.linspace(0, 5, n) dt = t[1] - t[0] L = FinDiff(0, dt, 2) + Identity() @@ -152,18 +153,18 @@ def test_1d_oscillator_free_dirichlet(self): eq = PDE(L, np.zeros_like(t), bc) u = eq.solve() - expected = np.cos(t)-(np.cos(5)-2)*np.sin(t)/np.sin(5) + expected = np.cos(t) - (np.cos(5) - 2) * np.sin(t) / np.sin(5) np.testing.assert_array_almost_equal(expected, u, decimal=4) def test_1d_damped_osc_driv_dirichlet(self): n = 100 - shape = n, + shape = (n,) t = np.linspace(0, 1, n) dt = t[1] - t[0] L = FinDiff(0, dt, 2) - FinDiff(0, dt) + Identity() - f = -3*np.exp(-t)*np.cos(t) + 2*np.exp(-t)*np.sin(t) + f = -3 * np.exp(-t) * np.cos(t) + 2 * np.exp(-t) * np.sin(t) - expected = np.exp(-t)*np.sin(t) + expected = np.exp(-t) * np.sin(t) bc = BoundaryConditions(shape) @@ -177,7 +178,7 @@ def test_1d_damped_osc_driv_dirichlet(self): def test_1d_oscillator_driv_neumann(self): n = 200 - shape = n, + shape = (n,) t = np.linspace(0, 1, n) dt = t[1] - t[0] L = FinDiff(0, dt, 2) - FinDiff(0, dt) + Identity() @@ -197,11 +198,11 @@ def test_1d_oscillator_driv_neumann(self): def test_1d_with_coeffs(self): n = 200 - shape = n, + shape = (n,) t = np.linspace(0, 1, n) dt = t[1] - t[0] L = Coef(t) * FinDiff(0, dt, 2) - f = 6*t**2 + f = 6 * t**2 bc = BoundaryConditions(shape) @@ -219,11 +220,15 @@ def test_mixed_equation__with_coeffs_2d(self): x, y = np.linspace(0, 1, shape[0]), np.linspace(0, 1, shape[1]) dx, dy = x[1] - x[0], y[1] - y[0] - X, Y = np.meshgrid(x, y, indexing='ij') - L = FinDiff(0, dx, 2) + Coef(X*Y) * FinDiff((0, dx, 1), (1, dy, 1)) + FinDiff(1, dy, 2) + X, Y = np.meshgrid(x, y, indexing="ij") + L = ( + FinDiff(0, dx, 2) + + Coef(X * Y) * FinDiff((0, dx, 1), (1, dy, 1)) + + FinDiff(1, dy, 2) + ) - expected = X ** 3 + Y ** 3 + 1 - f = 6*(X + Y) + expected = X**3 + Y**3 + 1 + f = 6 * (X + Y) bc = BoundaryConditions(shape) @@ -241,9 +246,13 @@ def test_2d_inhom_const_coefs_dirichlet_all(self): shape = (41, 50) (x, y), (dx, dy), (X, Y) = make_grid(shape, edges=[(-1, 1), (-1, 1)]) - expected = X**3 + Y**3 + X*Y + 1 + expected = X**3 + Y**3 + X * Y + 1 - L = Coef(3) * FinDiff(0, dx, 2) + Coef(2) * FinDiff((0, dx, 1), (1, dy, 1)) + FinDiff(1, dy, 2) + L = ( + Coef(3) * FinDiff(0, dx, 2) + + Coef(2) * FinDiff((0, dx, 1), (1, dy, 1)) + + FinDiff(1, dy, 2) + ) f = 2 + 18 * X + 6 * Y bc = BoundaryConditions(shape) @@ -261,10 +270,14 @@ def test_2d_inhom_var_coefs_dirichlet_all(self): shape = (41, 50) (x, y), (dx, dy), (X, Y) = make_grid(shape, edges=[(-1, 1), (-1, 1)]) - expected = X**3 + Y**3 + X*Y + 1 + expected = X**3 + Y**3 + X * Y + 1 - L = Coef(3*X) * FinDiff(0, dx, 2) + Coef(2*Y) * FinDiff((0, dx, 1), (1, dy, 1)) + FinDiff(1, dy, 2) - f = 18 * X**2 + 8*Y + L = ( + Coef(3 * X) * FinDiff(0, dx, 2) + + Coef(2 * Y) * FinDiff((0, dx, 1), (1, dy, 1)) + + FinDiff(1, dy, 2) + ) + f = 18 * X**2 + 8 * Y bc = BoundaryConditions(shape) bc[0, :] = expected @@ -283,11 +296,11 @@ def test_2d_inhom_var_coefs_with_identity_all_dirichlet(self): shape = (5, 5) (x, y), (dx, dy), (X, Y) = make_grid(shape, edges=[(-1, 1), (-1, 1)]) - expected = X**3 + Y**3 + X*Y + 1 + expected = X**3 + Y**3 + X * Y + 1 - #L = Coef(3*X) * FinDiff(0, dx, 2) + Coef(2*Y) * FinDiff((0, dx, 1), (1, dy, 1)) + FinDiff(1, dy, 2) + Coef(5*X*Y) * Identity() - L = Coef(5*X*Y) * FinDiff(0, dx, 2) #Identity() - #f = 18 * X**2 + 8*Y + 5*X*Y*expected + # L = Coef(3*X) * FinDiff(0, dx, 2) + Coef(2*Y) * FinDiff((0, dx, 1), (1, dy, 1)) + FinDiff(1, dy, 2) + Coef(5*X*Y) * Identity() + L = Coef(5 * X * Y) * FinDiff(0, dx, 2) # Identity() + # f = 18 * X**2 + 8*Y + 5*X*Y*expected mat = L.matrix(shape) print(mat) @@ -305,7 +318,9 @@ def test_2d_inhom_var_coefs_with_identity_all_dirichlet(self): def make_grid(shape, edges): - axes = tuple([np.linspace(edges[k][0], edges[k][1], shape[k]) for k in range(len(shape))]) - coords = np.meshgrid(*axes, indexing='ij') - spacings = [axes[k][1]-axes[k][0] for k in range(len(shape))] + axes = tuple( + [np.linspace(edges[k][0], edges[k][1], shape[k]) for k in range(len(shape))] + ) + coords = np.meshgrid(*axes, indexing="ij") + spacings = [axes[k][1] - axes[k][0] for k in range(len(shape))] return axes, spacings, coords diff --git a/test/test_scaling.py b/test/legacy/test_scaling.py similarity index 72% rename from test/test_scaling.py rename to test/legacy/test_scaling.py index 9275a2e..6605e4e 100644 --- a/test/test_scaling.py +++ b/test/legacy/test_scaling.py @@ -1,15 +1,17 @@ import sys -sys.path.insert(1, '..') + +sys.path.insert(1, "..") from math import log import unittest import numpy as np -from findiff.operators import FinDiff +from findiff.legacy.operators import FinDiff + class TestScaling(unittest.TestCase): def fit_1d(self, acc): - nx_list = 10 ** np.linspace(1.75, 2., 10) + nx_list = 10 ** np.linspace(1.75, 2.0, 10) Lx = np.pi log_err_list = [] @@ -17,7 +19,7 @@ def fit_1d(self, acc): for nx in nx_list: nx = int(nx) - x = np.linspace(0., Lx, nx) + x = np.linspace(0.0, Lx, nx) dx = x[1] - x[0] f = np.sin(x) d_dx = FinDiff(0, dx) @@ -42,7 +44,7 @@ def fit_2d(self, acc): x = np.linspace(0, Lx, nx) y = np.linspace(0, Ly, ny) dx, dy = x[1] - x[0], y[1] - y[0] - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") f = np.sin(X) * np.sin(Y) d_dx = FinDiff(0, dx) fx = d_dx(f, acc=acc) @@ -55,20 +57,20 @@ def fit_2d(self, acc): return fit[0] def test_1d_acc2(self): - self.assertAlmostEqual(2., self.fit_1d(acc=2), 1) + self.assertAlmostEqual(2.0, self.fit_1d(acc=2), 1) def test_1d_acc4(self): - self.assertAlmostEqual(4., self.fit_1d(acc=4), 1) + self.assertAlmostEqual(4.0, self.fit_1d(acc=4), 1) def test_1d_acc6(self): - self.assertAlmostEqual(6., self.fit_1d(acc=6), 1) + self.assertAlmostEqual(6.0, self.fit_1d(acc=6), 1) def test_2d_acc2(self): - self.assertAlmostEqual(2., self.fit_2d(acc=2), 1) + self.assertAlmostEqual(2.0, self.fit_2d(acc=2), 1) def test_2d_acc4(self): - self.assertAlmostEqual(4., self.fit_2d(acc=4), 1) + self.assertAlmostEqual(4.0, self.fit_2d(acc=4), 1) -if __name__ == '__main__': - unittest.main() \ No newline at end of file +if __name__ == "__main__": + unittest.main() diff --git a/test/test_stencils.py b/test/legacy/test_stencils.py similarity index 76% rename from test/test_stencils.py rename to test/legacy/test_stencils.py index 5230d5a..6472b91 100644 --- a/test/test_stencils.py +++ b/test/legacy/test_stencils.py @@ -3,7 +3,7 @@ import numpy as np from findiff import Identity, FinDiff -from findiff.stencils import Stencil +from findiff.legacy.stencils import Stencil class TestStencilOperations(unittest.TestCase): @@ -13,23 +13,36 @@ def test_solve_laplace_2d_with_5_points(self): stencil = Stencil(offsets, {(2, 0): 1, (0, 2): 1}) - expected = { - (0, 0): -4, - (-1, 0): 1, (1, 0): 1, (0, 1): 1, (0, -1): 1 - } + expected = {(0, 0): -4, (-1, 0): 1, (1, 0): 1, (0, 1): 1, (0, -1): 1} self.assertEqual(expected, stencil.values) self.assertEqual(2, stencil.accuracy) def test_solve_laplace_2d_with_9_points(self): - offsets = [(-1, 0), (0, 0), (1, 0), (0, 1), (0, -1), (-2, 0), (2, 0), (0, -2), (0, 2)] + offsets = [ + (-1, 0), + (0, 0), + (1, 0), + (0, 1), + (0, -1), + (-2, 0), + (2, 0), + (0, -2), + (0, 2), + ] stencil = Stencil(offsets, {(2, 0): 1, (0, 2): 1}) expected = { (0, 0): -5, - (-1, 0): 4 / 3., (1, 0): 4 / 3., (0, 1): 4 / 3., (0, -1): 4 / 3., - (-2, 0): -1 / 12., (2, 0): -1 / 12., (0, -2): -1 / 12., (0, 2): -1 / 12. + (-1, 0): 4 / 3.0, + (1, 0): 4 / 3.0, + (0, 1): 4 / 3.0, + (0, -1): 4 / 3.0, + (-2, 0): -1 / 12.0, + (2, 0): -1 / 12.0, + (0, -2): -1 / 12.0, + (0, 2): -1 / 12.0, } self.assertEqual(4, stencil.accuracy) @@ -42,10 +55,7 @@ def test_solve_laplace_2d_with_5_points_times_2(self): stencil = Stencil(offsets, {(2, 0): 2, (0, 2): 2}) - expected = { - (0, 0): -8, - (-1, 0): 2, (1, 0): 2, (0, 1): 2, (0, -1): 2 - } + expected = {(0, 0): -8, (-1, 0): 2, (1, 0): 2, (0, 1): 2, (0, -1): 2} self.assertEqual(expected, stencil.values) self.assertEqual(2, stencil.accuracy) @@ -55,10 +65,7 @@ def test_solve_laplace_2d_with_5_points_times_2_and_spacing(self): stencil = Stencil(offsets, {(2, 0): 2, (0, 2): 2}, spacings=(0.1, 0.1)) - expected = { - (0, 0): -800, - (-1, 0): 200, (1, 0): 200, (0, 1): 200, (0, -1): 200 - } + expected = {(0, 0): -800, (-1, 0): 200, (1, 0): 200, (0, 1): 200, (0, -1): 200} self.assertEqual(len(expected), len(stencil.values)) for off, coeff in stencil.values.items(): @@ -69,10 +76,10 @@ def test_apply_laplacian_laplacian(self): x = y = np.linspace(0, 1, 21) dx = dy = x[1] - x[0] - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") f = X**3 + Y**3 - expected = 6*X + 6*Y + expected = 6 * X + 6 * Y offsets = [(-1, 0), (0, 0), (1, 0), (0, 1), (0, -1)] stencil = Stencil(offsets, {(2, 0): 1, (0, 2): 1}, spacings=(dx, dy)) @@ -84,10 +91,10 @@ def test_apply_laplacian_laplacian_stencil_x_form(self): x = y = np.linspace(0, 1, 21) dx = dy = x[1] - x[0] - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") f = X**3 + Y**3 - expected = 6*X + 6*Y + expected = 6 * X + 6 * Y offsets = [(0, 0), (1, 1), (-1, -1), (1, -1), (-1, 1)] stencil = Stencil(offsets, {(2, 0): 1, (0, 2): 1}, spacings=(dx, dy)) @@ -99,7 +106,7 @@ def test_apply_laplacian_laplacian_stencil_outside_grid(self): x = y = np.linspace(0, 1, 21) dx = dy = x[1] - x[0] - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") f = X**3 + Y**3 @@ -117,10 +124,10 @@ def test_apply_mixed_deriv(self): x = y = np.linspace(0, 1, 101) dx = dy = x[1] - x[0] - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") - f = np.exp(-X**2-Y**2) - expected = 4*X*Y*f + f = np.exp(-(X**2) - Y**2) + expected = 4 * X * Y * f offsets = [(0, 0), (1, 1), (-1, -1), (1, -1), (-1, 1)] stencil = Stencil(offsets, partials={(1, 1): 1}, spacings=(dx, dy)) @@ -131,20 +138,20 @@ def test_apply_mixed_deriv(self): def test_laplace_1d_9points(self): x = np.linspace(0, 1, 101) f = x**3 - expected = 6*x + expected = 6 * x offsets = list(range(-4, 5)) - stencil = Stencil(offsets, partials={(2,): 1}, spacings=(x[1]-x[0],)) - at = 8, + stencil = Stencil(offsets, partials={(2,): 1}, spacings=(x[1] - x[0],)) + at = (8,) actual = stencil(f, at) self.assertAlmostEqual(expected[at], actual, places=5) def tests_apply_stencil_on_multislice(self): x = y = np.linspace(0, 1, 21) dx = dy = x[1] - x[0] - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") - f = X ** 3 + Y ** 3 - expected = 6*X + 6*Y + f = X**3 + Y**3 + expected = 6 * X + 6 * Y offsets = [(-1, 0), (0, 0), (1, 0), (0, 1), (0, -1)] stencil = Stencil(offsets, {(2, 0): 1, (0, 2): 1}, spacings=(dx, dy)) @@ -156,10 +163,10 @@ def tests_apply_stencil_on_multislice(self): def tests_apply_stencil_on_mask(self): x = y = np.linspace(0, 1, 21) dx = dy = x[1] - x[0] - X, Y = np.meshgrid(x, y, indexing='ij') + X, Y = np.meshgrid(x, y, indexing="ij") - f = X ** 3 + Y ** 3 - expected = 6*X + 6*Y + f = X**3 + Y**3 + expected = 6 * X + 6 * Y offsets = [(-1, 0), (0, 0), (1, 0), (0, 1), (0, -1)] stencil = Stencil(offsets, {(2, 0): 1, (0, 2): 1}, spacings=(dx, dy)) @@ -176,8 +183,11 @@ def test_helmholtz_stencil_issue_60(self): stencil_set = H.stencil((10,)) - expected = {('L',): {(0,): -1.0, (1,): 5.0, (2,): -4.0, (3,): 1.0}, ('C',): {(-1,): -1.0, (0,): 3.0, (1,): -1.0}, - ('H',): {(-3,): 1.0, (-2,): -4.0, (-1,): 5.0, (0,): -1.0}} + expected = { + ("L",): {(0,): -1.0, (1,): 5.0, (2,): -4.0, (3,): 1.0}, + ("C",): {(-1,): -1.0, (0,): 3.0, (1,): -1.0}, + ("H",): {(-3,): 1.0, (-2,): -4.0, (-1,): 5.0, (0,): -1.0}, + } actual = stencil_set.data self.assertEqual(len(expected), len(actual)) diff --git a/test/test_symbolic.py b/test/legacy/test_symbolic.py similarity index 95% rename from test/test_symbolic.py rename to test/legacy/test_symbolic.py index 83300c2..ea4b30e 100644 --- a/test/test_symbolic.py +++ b/test/legacy/test_symbolic.py @@ -1,8 +1,8 @@ import unittest -from sympy import IndexedBase, Symbol, Expr, Eq, symbols, latex +from sympy import IndexedBase, Symbol, symbols, latex -from findiff.symbolic import SymbolicMesh, SymbolicDiff +from findiff.legacy.symbolic import SymbolicMesh, SymbolicDiff class TestSymbolicMesh(unittest.TestCase): diff --git a/test/test_utils.py b/test/legacy/test_utils.py similarity index 89% rename from test/test_utils.py rename to test/legacy/test_utils.py index 95ca2e9..018f613 100644 --- a/test/test_utils.py +++ b/test/legacy/test_utils.py @@ -1,10 +1,10 @@ import sys -sys.path.insert(1, '..') + +sys.path.insert(1, "..") import unittest -import numpy as np -from numpy.testing import assert_array_almost_equal, assert_array_equal -from findiff.utils import * +from findiff.legacy.utils import * + class TestUtils(unittest.TestCase): diff --git a/test/test_vector.py b/test/legacy/test_vector.py similarity index 68% rename from test/test_vector.py rename to test/legacy/test_vector.py index e917fc8..b135dfd 100644 --- a/test/test_vector.py +++ b/test/legacy/test_vector.py @@ -1,10 +1,11 @@ import sys -sys.path.insert(1, '..') + +sys.path.insert(1, "..") import unittest from numpy.testing import assert_array_almost_equal import numpy as np -from findiff.vector import Gradient, Divergence, Curl, Laplacian +from findiff.legacy.vector import Gradient, Divergence, Curl, Laplacian class TestGradient(unittest.TestCase): @@ -12,11 +13,13 @@ class TestGradient(unittest.TestCase): def test_3d_gradient_on_scalar_func(self): axes, h, [X, Y, Z] = init_mesh(3, (50, 50, 50)) f = np.sin(X) * np.sin(Y) * np.sin(Z) - grad_f_ex = np.array([ - np.cos(X) * np.sin(Y) * np.sin(Z), - np.sin(X) * np.cos(Y) * np.sin(Z), - np.sin(X) * np.sin(Y) * np.cos(Z), - ]) + grad_f_ex = np.array( + [ + np.cos(X) * np.sin(Y) * np.sin(Z), + np.sin(X) * np.cos(Y) * np.sin(Z), + np.sin(X) * np.sin(Y) * np.cos(Z), + ] + ) grad = Gradient(spac=h, acc=4) grad_f = grad(f) assert_array_almost_equal(grad_f, grad_f_ex) @@ -24,11 +27,13 @@ def test_3d_gradient_on_scalar_func(self): def test_spacing_with_h(self): axes, h, [X, Y, Z] = init_mesh(3, (50, 50, 50)) f = np.sin(X) * np.sin(Y) * np.sin(Z) - grad_f_ex = np.array([ - np.cos(X) * np.sin(Y) * np.sin(Z), - np.sin(X) * np.cos(Y) * np.sin(Z), - np.sin(X) * np.sin(Y) * np.cos(Z), - ]) + grad_f_ex = np.array( + [ + np.cos(X) * np.sin(Y) * np.sin(Z), + np.sin(X) * np.cos(Y) * np.sin(Z), + np.sin(X) * np.sin(Y) * np.cos(Z), + ] + ) grad = Gradient(h=h, acc=4) grad_f = grad(f) assert_array_almost_equal(grad_f, grad_f_ex) @@ -36,33 +41,37 @@ def test_spacing_with_h(self): def test_3d_gradient_on_scalar_func_non_uni(self): axes, h, [X, Y, Z] = init_mesh(3, (50, 50, 50)) f = np.sin(X) * np.sin(Y) * np.sin(Z) - grad_f_ex = np.array([ - np.cos(X) * np.sin(Y) * np.sin(Z), - np.sin(X) * np.cos(Y) * np.sin(Z), - np.sin(X) * np.sin(Y) * np.cos(Z), - ]) + grad_f_ex = np.array( + [ + np.cos(X) * np.sin(Y) * np.sin(Z), + np.sin(X) * np.cos(Y) * np.sin(Z), + np.sin(X) * np.sin(Y) * np.cos(Z), + ] + ) grad = Gradient(coords=axes, acc=4) grad_f = grad(f) assert_array_almost_equal(grad_f, grad_f_ex) def test_3d_gradient_on_vector_func_should_fail(self): axes, h, [X, Y, Z] = init_mesh(3, (50, 50, 50)) - f = np.array([np.sin(X) * np.sin(Y) * np.sin(Z), - np.sin(X) * np.sin(Y) * np.sin(Z) - ]) + f = np.array( + [np.sin(X) * np.sin(Y) * np.sin(Z), np.sin(X) * np.sin(Y) * np.sin(Z)] + ) grad = Gradient(spac=h, acc=4) self.assertRaises(ValueError, grad, f) + class TestDivergence(unittest.TestCase): def test_3d_divergence_on_vector_func(self): axes, h, [X, Y, Z] = init_mesh(3, (50, 50, 50)) f = np.array([np.sin(X) * np.sin(Y) * np.sin(Z)] * 3) assert f.shape == (3, 50, 50, 50) - div_f_ex = \ - np.cos(X) * np.sin(Y) * np.sin(Z) +\ - np.sin(X) * np.cos(Y) * np.sin(Z) +\ - np.sin(X) * np.sin(Y) * np.cos(Z) + div_f_ex = ( + np.cos(X) * np.sin(Y) * np.sin(Z) + + np.sin(X) * np.cos(Y) * np.sin(Z) + + np.sin(X) * np.sin(Y) * np.cos(Z) + ) div = Divergence(spac=h, acc=4) div_f = div(f) @@ -82,11 +91,13 @@ def test_curl_on_3d_vector_func(self): axes, h, [X, Y, Z] = init_mesh(3, (50, 50, 50)) f = np.array([np.sin(X) * np.sin(Y) * np.sin(Z)] * 3) assert f.shape == (3, 50, 50, 50) - curl_f_ex = np.array([ - np.sin(X) * np.cos(Y) * np.sin(Z) - np.sin(X) * np.sin(Y) * np.cos(Z), - np.sin(X) * np.sin(Y) * np.cos(Z) - np.cos(X) * np.sin(Y) * np.sin(Z), - np.cos(X) * np.sin(Y) * np.sin(Z) - np.sin(X) * np.cos(Y) * np.sin(Z), - ]) + curl_f_ex = np.array( + [ + np.sin(X) * np.cos(Y) * np.sin(Z) - np.sin(X) * np.sin(Y) * np.cos(Z), + np.sin(X) * np.sin(Y) * np.cos(Z) - np.cos(X) * np.sin(Y) * np.sin(Z), + np.cos(X) * np.sin(Y) * np.sin(Z) - np.sin(X) * np.cos(Y) * np.sin(Z), + ] + ) curl = Curl(spac=h, acc=4) curl_f = curl(f) assert_array_almost_equal(curl_f, curl_f_ex) @@ -117,5 +128,5 @@ def init_mesh(ndims, npoints): return axes, h, mesh -if __name__ == '__main__': - unittest.main() \ No newline at end of file +if __name__ == "__main__": + unittest.main() From b2bc5f5797f06b6a98efabf57b019874b71dec97 Mon Sep 17 00:00:00 2001 From: maroba Date: Mon, 9 Dec 2024 19:03:13 +0100 Subject: [PATCH 2/6] ... --- CHANGES.md | 8 ++ README.md | 2 +- findiff/__init__.py | 10 +- findiff/current/__init__.py | 3 + findiff/current/operators.py | 122 ++++++++++++++++ findiff/legacy/__init__.py | 23 +++ test/legacy/test_bugs.py | 101 -------------- test/legacy/test_diff.py | 173 ----------------------- {test => tests}/__init__.py | 0 {test => tests}/legacy/__init__.py | 0 tests/legacy/test_bugs.py | 143 +++++++++++++++++++ {test => tests}/legacy/test_coefs.py | 2 +- {test => tests}/legacy/test_findiff.py | 0 {test => tests}/legacy/test_grids.py | 0 {test => tests}/legacy/test_pde.py | 0 {test => tests}/legacy/test_scaling.py | 0 {test => tests}/legacy/test_stencils.py | 3 +- {test => tests}/legacy/test_symbolic.py | 0 {test => tests}/legacy/test_utils.py | 0 {test => tests}/legacy/test_vector.py | 0 tests/test_diff.py | 176 +++++++++++++++++++++++ tests/test_operators.py | 177 ++++++++++++++++++++++++ 22 files changed, 660 insertions(+), 283 deletions(-) create mode 100644 findiff/current/__init__.py create mode 100644 findiff/current/operators.py delete mode 100644 test/legacy/test_bugs.py delete mode 100644 test/legacy/test_diff.py rename {test => tests}/__init__.py (100%) rename {test => tests}/legacy/__init__.py (100%) create mode 100644 tests/legacy/test_bugs.py rename {test => tests}/legacy/test_coefs.py (99%) rename {test => tests}/legacy/test_findiff.py (100%) rename {test => tests}/legacy/test_grids.py (100%) rename {test => tests}/legacy/test_pde.py (100%) rename {test => tests}/legacy/test_scaling.py (100%) rename {test => tests}/legacy/test_stencils.py (99%) rename {test => tests}/legacy/test_symbolic.py (100%) rename {test => tests}/legacy/test_utils.py (100%) rename {test => tests}/legacy/test_vector.py (100%) create mode 100644 tests/test_diff.py create mode 100644 tests/test_operators.py diff --git a/CHANGES.md b/CHANGES.md index 32a1a22..ae00975 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,5 +1,13 @@ # Change Log +## Version 0.10.x + +- Create symbolic representations of finite difference schemes + +## Version 0.9.x + +- Generate differential operators for generic stencils + ## Version 0.8.x - Formulate and solve partial differential equations diff --git a/README.md b/README.md index 07016c6..402d776 100644 --- a/README.md +++ b/README.md @@ -387,5 +387,5 @@ python setup.py develop From the console: ``` -python -m unittest discover test +python -m unittest discover tests ``` diff --git a/findiff/__init__.py b/findiff/__init__.py index d15a5f5..6d8763e 100644 --- a/findiff/__init__.py +++ b/findiff/__init__.py @@ -15,12 +15,10 @@ - New in version 0.8: Solve partial differential equations with Dirichlet or Neumann boundary conditions - New in version 0.9: Generate differential operators for generic stencils - New in version 0.10: Create symbolic representations of finite difference schemes +- Version 1.0: Completely remodeled API (backward compatibility is maintained, though) """ -__version__ = "0.10.2" +__version__ = "1.0.0.dev" -from .legacy.coefs import coefficients -from .legacy.operators import FinDiff, Coef, Identity, Coefficient -from .legacy.pde import PDE, BoundaryConditions -from .legacy.symbolic import SymbolicMesh, SymbolicDiff -from .legacy.vector import Gradient, Divergence, Curl, Laplacian + +from .current import Diff, Identity diff --git a/findiff/current/__init__.py b/findiff/current/__init__.py new file mode 100644 index 0000000..6f8596a --- /dev/null +++ b/findiff/current/__init__.py @@ -0,0 +1,3 @@ +API_VERSION = 1 + +from .operators import * diff --git a/findiff/current/operators.py b/findiff/current/operators.py new file mode 100644 index 0000000..f0ab4e7 --- /dev/null +++ b/findiff/current/operators.py @@ -0,0 +1,122 @@ +import numbers +from abc import ABC, abstractmethod + +import numpy as np +from scipy import sparse + + +class Node(ABC): + __array_priority__ = 100 # Make sure custom multiplication is called over numpy's + + @abstractmethod + def __call__(self, f, *args, **kwargs): + pass + + def __add__(self, other): + return Add(self, other) + + def __radd__(self, other): + return Add(self, other) + + def __sub__(self, other): + return Add(FieldOperator(-1) * other, self) + + def __rsub__(self, other): + return Add(FieldOperator(-1) * other, self) + + def __mul__(self, other): + return Mul(self, other) + + def __rmul__(self, other): + return Mul(other, self) + + +class FieldOperator(Node): + """An operator that multiplies an array pointwise.""" + + def __init__(self, value): + self.value = value + + def __call__(self, f, *args, **kwargs): + if isinstance(f, (numbers.Number, np.ndarray)): + return self.value * f + return self.value * super().__call__(f, *args, **kwargs) + + +class ScalarOperator(FieldOperator): + def __init__(self, value): + if not isinstance(value, numbers.Number): + raise ValueError("Expected number, got " + str(type(value))) + super().__init__(value) + + +class Identity(ScalarOperator): + def __init__(self): + super().__init__(1) + + +class Add(Node): + def __init__(self, left, right): + if isinstance(left, (numbers.Number, np.ndarray)): + left = FieldOperator(left) + if isinstance(right, (numbers.Number, np.ndarray)): + right = FieldOperator(right) + self.left = left + self.right = right + + def __call__(self, f, *args, **kwargs): + return self.left(f, *args, **kwargs) + self.right(f, *args, **kwargs) + + +class Mul(Node): + def __init__(self, left, right): + if isinstance(left, (numbers.Number, np.ndarray)): + left = FieldOperator(left) + if isinstance(right, (numbers.Number, np.ndarray)): + right = FieldOperator(right) + self.left = left + self.right = right + + def __call__(self, f, *args, **kwargs): + return self.left(self.right(f, *args, **kwargs), *args, **kwargs) + + +class Diff(Node): + + DEFAULT_ACC = 2 + + def __init__(self, spacing, order=1, axis=0, acc=DEFAULT_ACC): + if spacing <= 0: + raise ValueError("spacing must be greater than zero") + if order <= 0: + raise ValueError("order must be greater than zero") + + self.spacing = spacing + self.axis = axis + self.order = order + self.acc = acc + + self._fd = None + + def __call__(self, f, *args, **kwargs): + + if "acc" in kwargs: + # allow to pass down new accuracy deep in expression tree + new_acc = kwargs["acc"] + + if new_acc != self.acc: + self._fd = None + + self.acc = new_acc + + if self._fd is None: + self._build_differentiator() + + if isinstance(f, Node): + f = f(*args, **kwargs) + return self._fd(f) + + def _build_differentiator(self): + from findiff.legacy import FinDiff + + self._fd = FinDiff(self.axis, self.spacing, self.order, acc=self.acc) diff --git a/findiff/legacy/__init__.py b/findiff/legacy/__init__.py index e69de29..3647c29 100644 --- a/findiff/legacy/__init__.py +++ b/findiff/legacy/__init__.py @@ -0,0 +1,23 @@ +from .coefs import coefficients +from .operators import FinDiff, Coef, Coefficient, Identity +from .pde import PDE, BoundaryConditions +from .symbolic import SymbolicMesh, SymbolicDiff +from .vector import Gradient, Divergence, Curl, Laplacian + +API_VERSION = 0 + +__all__ = [ + "coefficients", + "FinDiff", + "Identity", + "Coef", + "Coefficient", + "PDE", + "BoundaryConditions", + "SymbolicMesh", + "SymbolicDiff", + "Gradient", + "Divergence", + "Curl", + "Laplacian", +] diff --git a/test/legacy/test_bugs.py b/test/legacy/test_bugs.py deleted file mode 100644 index 9bdc82d..0000000 --- a/test/legacy/test_bugs.py +++ /dev/null @@ -1,101 +0,0 @@ -import sys - -sys.path.insert(1, '..') - -import unittest -import numpy as np -from findiff import FinDiff -import findiff - - -class TestOldBugs(unittest.TestCase): - - def test_findiff_should_raise_exception_when_applied_to_unevaluated_function(self): - def f(x, y): - return 5 * x ** 2 - 5 * x + 10 * y ** 2 - 10 * y - - d_dx = FinDiff(1, 0.01) - self.assertRaises(ValueError, lambda ff: d_dx(ff), f) - - def test_matrix_representation_doesnt_work_for_order_greater_2_issue_24(self): - x = np.zeros((10)) - d3_dx3 = FinDiff((0, 1, 3)) - mat = d3_dx3.matrix(x.shape) - self.assertAlmostEqual(-2.5, mat[0, 0]) - self.assertAlmostEqual(-2.5, mat[1, 1]) - self.assertAlmostEqual(-0.5, mat[2, 0]) - - def test_high_accuracy_results_in_type_error(self): - # in issue 25 the following line resulted in a TypeError - findiff.coefficients(deriv=1, acc=16) - - def test_matrix_repr_with_different_accs(self): - # issue 28 - shape = (11,) - d1 = findiff.FinDiff(0, 1, 2).matrix(shape) - d2 = findiff.FinDiff(0, 1, 2, acc=4).matrix(shape) - - self.assertTrue(np.max(np.abs((d1 - d2).toarray())) > 1) - - x = np.linspace(0, 10, 11) - f = x ** 2 - df = d2.dot(f) - np.testing.assert_almost_equal(2 * np.ones_like(f), df) - - def test_accuracy_should_be_passed_down_to_stencil(self): - # issue 31 - - shape = 11, 11 - dx = 1. - d1x = FinDiff(0, dx, 1, acc=4) - stencil1 = d1x.stencil(shape) - - expected = { - ('L', 'L'): {(0, 0): -2.083333333333331, (1, 0): 3.9999999999999916, (2, 0): -2.999999999999989, - (3, 0): 1.3333333333333268, (4, 0): -0.24999999999999858}, - ('L', 'C'): {(0, 0): -2.083333333333331, (1, 0): 3.9999999999999916, (2, 0): -2.999999999999989, - (3, 0): 1.3333333333333268, (4, 0): -0.24999999999999858}, - ('L', 'H'): {(0, 0): -2.083333333333331, (1, 0): 3.9999999999999916, (2, 0): -2.999999999999989, - (3, 0): 1.3333333333333268, (4, 0): -0.24999999999999858}, - ('C', 'L'): {(-2, 0): 0.08333333333333333, (-1, 0): -0.6666666666666666, - (1, 0): 0.6666666666666666, (2, 0): -0.08333333333333333}, - ('C', 'C'): {(-2, 0): 0.08333333333333333, (-1, 0): -0.6666666666666666, - (1, 0): 0.6666666666666666, (2, 0): -0.08333333333333333}, - ('C', 'H'): {(-2, 0): 0.08333333333333333, (-1, 0): -0.6666666666666666, - (1, 0): 0.6666666666666666, (2, 0): -0.08333333333333333}, - ('H', 'L'): {(-4, 0): 0.24999999999999958, (-3, 0): -1.3333333333333313, (-2, 0): 2.9999999999999956, - (-1, 0): -3.999999999999996, (0, 0): 2.0833333333333317}, - ('H', 'C'): {(-4, 0): 0.24999999999999958, (-3, 0): -1.3333333333333313, (-2, 0): 2.9999999999999956, - (-1, 0): -3.999999999999996, (0, 0): 2.0833333333333317}, - ('H', 'H'): {(-4, 0): 0.24999999999999958, (-3, 0): -1.3333333333333313, (-2, 0): 2.9999999999999956, - (-1, 0): -3.999999999999996, (0, 0): 2.0833333333333317}, - } - - for char_pt in stencil1.data: - stl = stencil1.data[char_pt] - self.assert_dict_almost_equal(expected[char_pt], stl) - - d1x = FinDiff(0, dx, 1, acc=4) - stencil1 = d1x.stencil(shape) - for char_pt in stencil1.data: - stl = stencil1.data[char_pt] - self.assert_dict_almost_equal(expected[char_pt], stl) - - def test_order_as_numpy_integer(self): - - order = np.ones(3, dtype=np.int32)[0] - d_dx = FinDiff(0, 0.1, order) # raised an AssertionError with the bug - - np.testing.assert_allclose(d_dx(np.linspace(0, 1, 11)), np.ones(11)) - - - def assert_dict_almost_equal(self, actual, expected, places=7): - if len(actual) != len(expected): - return False - - for key, value in actual.items(): - self.assertAlmostEqual(actual[key], expected[key], places=places) - - -if __name__ == '__main__': - unittest.main() diff --git a/test/legacy/test_diff.py b/test/legacy/test_diff.py deleted file mode 100644 index ffcd105..0000000 --- a/test/legacy/test_diff.py +++ /dev/null @@ -1,173 +0,0 @@ -import sys - -sys.path.insert(1, "..") - -import unittest -from numpy.testing import assert_array_almost_equal -from findiff.legacy.diff import * - -# import matplotlib.pyplot as plt -from numpy import cos, sin - - -class DiffTest(unittest.TestCase): - - def test_partial_d_dx(self): - shape = (101,) - x, dx = make_grid(shape, (0, 1)) - - u = x**2 - expected = 2 * x - - fd = Diff(0, 1) - actual = fd(u, dx) - - assert_array_almost_equal(expected, actual) - - def test_partial_d2_dx2(self): - shape = (101,) - x, dx = make_grid(shape, (0, 1)) - - u = x**2 - expected = 2 - - fd = Diff(0, 2) - actual = fd(u, dx) - - assert_array_almost_equal(expected, actual) - - def test_partial_d_dx_acc(self): - shape = (11,) - x, dx = make_grid(shape, (0, 1)) - - u = x**3 - expected = 3 * x**2 - - fd = Diff(0, 1) - actual = fd(u, dx) - np.testing.assert_raises( - AssertionError, assert_array_almost_equal, expected, actual - ) - - actual = fd(u, dx, acc=4) - assert_array_almost_equal(expected, actual) - - def test_partial_d2_dxdy(self): - - shape = 50, 50 - (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) - u = np.sin(X) * np.sin(Y) - - fd = Diff(0, 1) * Diff(1, 1) - actual = fd(u, h={0: dx, 1: dy}) - - expected = np.cos(X) * np.cos(Y) - assert_array_almost_equal(expected, actual, decimal=3) - - fd = Diff(1, 1) * Diff(0, 1) - actual = fd(u, h={0: dx, 1: dy}) - assert_array_almost_equal(expected, actual, decimal=3) - - def test_partial_d_dx_on_2d_array(self): - - shape = 50, 50 - (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) - u = np.sin(X) * np.sin(Y) - - fd = Diff(1, 1) - actual = fd(u, dy) - - expected = np.sin(X) * np.cos(Y) - assert_array_almost_equal(expected, actual, decimal=3) - - def test_add_two_diffs(self): - - shape = 50, 50 - (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) - u = np.sin(X) * np.sin(Y) - - fd = Diff(0, 2) + Diff(1, 2) - actual = fd(u, dy) - - expected = -2 * np.sin(X) * np.sin(Y) - assert_array_almost_equal(expected, actual, decimal=3) - - def test_sub_two_diffs(self): - - shape = 50, 50 - (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) - u = np.sin(X) * np.sin(Y) - - fd = Diff(0, 2) - Diff(1, 2) - actual = fd(u, dy) - - expected = np.zeros(shape) - assert_array_almost_equal(expected, actual, decimal=3) - - def test_multiply_with_scalar(self): - - shape = 50, 50 - (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) - u = np.sin(X) * np.sin(Y) - - fd = 3 * Diff(0, 2) - actual = fd(u, dy, acc=4) - - expected = -3 * np.sin(X) * np.sin(Y) - assert_array_almost_equal(expected, actual) - - def test_multiply_with_variable_from_left(self): - - shape = 50, 50 - (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) - u = np.sin(X) * np.sin(Y) - - fd = Coef(3 * X * Y) * Diff(0, 2) - actual = fd(u, dy, acc=4) - - expected = -3 * X * Y * np.sin(X) * np.sin(Y) - assert_array_almost_equal(expected, actual) - - def test_identity(self): - - shape = 50, 50 - (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) - u = np.sin(X) * np.sin(Y) - - fd = Diff(0, 2) + Id() - actual = fd(u, dy, acc=4) - - expected = np.zeros_like(u) - assert_array_almost_equal(expected, actual) - - def test_linear_comb(self): - shape = 50, 50 - (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) - u = np.sin(X) * np.sin(Y) - - fd = Coef(3 * X) * Diff(0, 2) - Coef(Y) * (Diff(1, 2) + Diff(0, 1)) - actual = fd(u, h={0: dx, 1: dy}, acc=4) - - expected = -3 * X * u - Y * (-sin(X) * sin(Y) + cos(X) * sin(Y)) - assert_array_almost_equal(expected, actual) - - def test_apply_base_binary_operator_raises_exception(self): - - op = BinaryOperator(None, None) - with self.assertRaises(NotImplementedError): - op(None) - - -def make_grid(shape, edges): - - if len(shape) == 1: - x = np.linspace(*edges, shape[0]) - dx = x[1] - x[0] - return x, dx - - axes = tuple( - [np.linspace(edges[k][0], edges[k][1], shape[k]) for k in range(len(shape))] - ) - coords = np.meshgrid(*axes, indexing="ij") - spacings = [axes[k][1] - axes[k][0] for k in range(len(shape))] - return axes, spacings, coords diff --git a/test/__init__.py b/tests/__init__.py similarity index 100% rename from test/__init__.py rename to tests/__init__.py diff --git a/test/legacy/__init__.py b/tests/legacy/__init__.py similarity index 100% rename from test/legacy/__init__.py rename to tests/legacy/__init__.py diff --git a/tests/legacy/test_bugs.py b/tests/legacy/test_bugs.py new file mode 100644 index 0000000..d4a72c3 --- /dev/null +++ b/tests/legacy/test_bugs.py @@ -0,0 +1,143 @@ +import sys + +sys.path.insert(1, "..") + +import unittest +import numpy as np + +from findiff.legacy import FinDiff +import findiff.legacy as findiff + + +class TestOldBugs(unittest.TestCase): + + def test_findiff_should_raise_exception_when_applied_to_unevaluated_function(self): + def f(x, y): + return 5 * x**2 - 5 * x + 10 * y**2 - 10 * y + + d_dx = FinDiff(1, 0.01) + self.assertRaises(ValueError, lambda ff: d_dx(ff), f) + + def test_matrix_representation_doesnt_work_for_order_greater_2_issue_24(self): + x = np.zeros((10)) + d3_dx3 = FinDiff((0, 1, 3)) + mat = d3_dx3.matrix(x.shape) + self.assertAlmostEqual(-2.5, mat[0, 0]) + self.assertAlmostEqual(-2.5, mat[1, 1]) + self.assertAlmostEqual(-0.5, mat[2, 0]) + + def test_high_accuracy_results_in_type_error(self): + # in issue 25 the following line resulted in a TypeError + findiff.coefficients(deriv=1, acc=16) + + def test_matrix_repr_with_different_accs(self): + # issue 28 + shape = (11,) + d1 = findiff.FinDiff(0, 1, 2).matrix(shape) + d2 = findiff.FinDiff(0, 1, 2, acc=4).matrix(shape) + + self.assertTrue(np.max(np.abs((d1 - d2).toarray())) > 1) + + x = np.linspace(0, 10, 11) + f = x**2 + df = d2.dot(f) + np.testing.assert_almost_equal(2 * np.ones_like(f), df) + + def test_accuracy_should_be_passed_down_to_stencil(self): + # issue 31 + + shape = 11, 11 + dx = 1.0 + d1x = FinDiff(0, dx, 1, acc=4) + stencil1 = d1x.stencil(shape) + + expected = { + ("L", "L"): { + (0, 0): -2.083333333333331, + (1, 0): 3.9999999999999916, + (2, 0): -2.999999999999989, + (3, 0): 1.3333333333333268, + (4, 0): -0.24999999999999858, + }, + ("L", "C"): { + (0, 0): -2.083333333333331, + (1, 0): 3.9999999999999916, + (2, 0): -2.999999999999989, + (3, 0): 1.3333333333333268, + (4, 0): -0.24999999999999858, + }, + ("L", "H"): { + (0, 0): -2.083333333333331, + (1, 0): 3.9999999999999916, + (2, 0): -2.999999999999989, + (3, 0): 1.3333333333333268, + (4, 0): -0.24999999999999858, + }, + ("C", "L"): { + (-2, 0): 0.08333333333333333, + (-1, 0): -0.6666666666666666, + (1, 0): 0.6666666666666666, + (2, 0): -0.08333333333333333, + }, + ("C", "C"): { + (-2, 0): 0.08333333333333333, + (-1, 0): -0.6666666666666666, + (1, 0): 0.6666666666666666, + (2, 0): -0.08333333333333333, + }, + ("C", "H"): { + (-2, 0): 0.08333333333333333, + (-1, 0): -0.6666666666666666, + (1, 0): 0.6666666666666666, + (2, 0): -0.08333333333333333, + }, + ("H", "L"): { + (-4, 0): 0.24999999999999958, + (-3, 0): -1.3333333333333313, + (-2, 0): 2.9999999999999956, + (-1, 0): -3.999999999999996, + (0, 0): 2.0833333333333317, + }, + ("H", "C"): { + (-4, 0): 0.24999999999999958, + (-3, 0): -1.3333333333333313, + (-2, 0): 2.9999999999999956, + (-1, 0): -3.999999999999996, + (0, 0): 2.0833333333333317, + }, + ("H", "H"): { + (-4, 0): 0.24999999999999958, + (-3, 0): -1.3333333333333313, + (-2, 0): 2.9999999999999956, + (-1, 0): -3.999999999999996, + (0, 0): 2.0833333333333317, + }, + } + + for char_pt in stencil1.data: + stl = stencil1.data[char_pt] + self.assert_dict_almost_equal(expected[char_pt], stl) + + d1x = FinDiff(0, dx, 1, acc=4) + stencil1 = d1x.stencil(shape) + for char_pt in stencil1.data: + stl = stencil1.data[char_pt] + self.assert_dict_almost_equal(expected[char_pt], stl) + + def test_order_as_numpy_integer(self): + + order = np.ones(3, dtype=np.int32)[0] + d_dx = FinDiff(0, 0.1, order) # raised an AssertionError with the bug + + np.testing.assert_allclose(d_dx(np.linspace(0, 1, 11)), np.ones(11)) + + def assert_dict_almost_equal(self, actual, expected, places=7): + if len(actual) != len(expected): + return False + + for key, value in actual.items(): + self.assertAlmostEqual(actual[key], expected[key], places=places) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/legacy/test_coefs.py b/tests/legacy/test_coefs.py similarity index 99% rename from test/legacy/test_coefs.py rename to tests/legacy/test_coefs.py index 667a907..12af237 100644 --- a/test/legacy/test_coefs.py +++ b/tests/legacy/test_coefs.py @@ -3,7 +3,7 @@ sys.path.insert(1, "..") import unittest -from findiff import coefficients +from findiff.legacy import coefficients from findiff.legacy.coefs import coefficients_non_uni from findiff.legacy.coefs import calc_coefs diff --git a/test/legacy/test_findiff.py b/tests/legacy/test_findiff.py similarity index 100% rename from test/legacy/test_findiff.py rename to tests/legacy/test_findiff.py diff --git a/test/legacy/test_grids.py b/tests/legacy/test_grids.py similarity index 100% rename from test/legacy/test_grids.py rename to tests/legacy/test_grids.py diff --git a/test/legacy/test_pde.py b/tests/legacy/test_pde.py similarity index 100% rename from test/legacy/test_pde.py rename to tests/legacy/test_pde.py diff --git a/test/legacy/test_scaling.py b/tests/legacy/test_scaling.py similarity index 100% rename from test/legacy/test_scaling.py rename to tests/legacy/test_scaling.py diff --git a/test/legacy/test_stencils.py b/tests/legacy/test_stencils.py similarity index 99% rename from test/legacy/test_stencils.py rename to tests/legacy/test_stencils.py index 6472b91..cf89727 100644 --- a/test/legacy/test_stencils.py +++ b/tests/legacy/test_stencils.py @@ -2,7 +2,8 @@ import numpy as np -from findiff import Identity, FinDiff +from findiff.legacy import Identity, FinDiff + from findiff.legacy.stencils import Stencil diff --git a/test/legacy/test_symbolic.py b/tests/legacy/test_symbolic.py similarity index 100% rename from test/legacy/test_symbolic.py rename to tests/legacy/test_symbolic.py diff --git a/test/legacy/test_utils.py b/tests/legacy/test_utils.py similarity index 100% rename from test/legacy/test_utils.py rename to tests/legacy/test_utils.py diff --git a/test/legacy/test_vector.py b/tests/legacy/test_vector.py similarity index 100% rename from test/legacy/test_vector.py rename to tests/legacy/test_vector.py diff --git a/tests/test_diff.py b/tests/test_diff.py new file mode 100644 index 0000000..782e00e --- /dev/null +++ b/tests/test_diff.py @@ -0,0 +1,176 @@ +import numpy as np + +# import matplotlib.pyplot as plt +from numpy import cos, sin +from numpy.testing import assert_array_almost_equal + +from findiff import Diff, Identity + + +def test_partial_d_dx(): + shape = (101,) + x, dx = make_grid(shape, (0, 1)) + + u = x**2 + expected = 2 * x + + fd = Diff(dx, 1, 0) + actual = fd(u) + + assert_array_almost_equal(expected, actual) + + +def test_partial_d2_dx2(): + shape = (101,) + x, dx = make_grid(shape, (0, 1)) + + u = x**2 + expected = 2 + + fd = Diff( + dx, + 2, + 0, + ) + actual = fd(u, dx) + + assert_array_almost_equal(expected, actual) + + +def test_partial_d_dx_acc(): + shape = (11,) + x, dx = make_grid(shape, (0, 1)) + + u = x**3 + expected = 3 * x**2 + + fd = Diff(dx, 1, 0) + actual = fd(u, dx) + np.testing.assert_raises( + AssertionError, assert_array_almost_equal, expected, actual + ) + + actual = fd(u, acc=4) + assert_array_almost_equal(expected, actual) + + +def test_partial_d2_dxdy(): + + shape = 50, 50 + (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) + u = np.sin(X) * np.sin(Y) + + fd = Diff(dx, 1, 0) * Diff(dy, 1, 1) + actual = fd(u, h={0: dx, 1: dy}) + + expected = np.cos(X) * np.cos(Y) + assert_array_almost_equal(expected, actual, decimal=3) + + fd = Diff(dy, 1, 1) * Diff(dy, 1, 0) + actual = fd(u) + assert_array_almost_equal(expected, actual, decimal=3) + + +def test_partial_d_dx_on_2d_array(): + + shape = 50, 50 + (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) + u = np.sin(X) * np.sin(Y) + + fd = Diff(dy, 1, 1) + actual = fd(u) + + expected = np.sin(X) * np.cos(Y) + assert_array_almost_equal(expected, actual, decimal=3) + + +def test_add_two_diffs(): + + shape = 50, 50 + (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) + u = np.sin(X) * np.sin(Y) + + fd = Diff(dx, 2, 0) + Diff(dy, 2, 1) + actual = fd(u) + + expected = -2 * np.sin(X) * np.sin(Y) + assert_array_almost_equal(expected, actual, decimal=3) + + +def test_sub_two_diffs(): + + shape = 50, 50 + (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) + u = np.sin(X) * np.sin(Y) + + fd = Diff(dx, 2, 0) - Diff(dy, 2, 1) + actual = fd(u) + + expected = np.zeros(shape) + assert_array_almost_equal(expected, actual, decimal=3) + + +def test_multiply_with_scalar(): + + shape = 50, 50 + (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) + u = np.sin(X) * np.sin(Y) + + fd = 3 * Diff(dx, 2, 0) + actual = fd(u, acc=4) + + expected = -3 * np.sin(X) * np.sin(Y) + assert_array_almost_equal(expected, actual) + + +def test_multiply_with_variable_from_left(): + + shape = 50, 50 + (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) + u = np.sin(X) * np.sin(Y) + + fd = 3 * X * Y * Diff(dx, 2, 0) + actual = fd(u, acc=4) + + expected = -3 * X * Y * np.sin(X) * np.sin(Y) + assert_array_almost_equal(expected, actual) + + +def test_identity(): + + shape = 50, 50 + (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) + u = np.sin(X) * np.sin(Y) + + fd = Diff(dx, 2, 0) + Identity() + actual = fd(u, acc=4) + + expected = np.zeros_like(u) + assert_array_almost_equal(expected, actual) + + +def test_linear_comb(): + shape = 50, 50 + (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) + u = np.sin(X) * np.sin(Y) + + fd = 3 * X * Diff(dx, 2, 0) - Y * (Diff(dy, 2, 1) + Diff(dx, 1, 0)) + actual = fd(u, acc=4) + + expected = -3 * X * u - Y * (-sin(X) * sin(Y) + cos(X) * sin(Y)) + assert_array_almost_equal(expected, actual) + + +def make_grid(shape, edges): + + if len(shape) == 1: + x = np.linspace(*edges, shape[0]) + dx = x[1] - x[0] + return x, dx + + axes = tuple( + [np.linspace(edges[k][0], edges[k][1], shape[k]) for k in range(len(shape))] + ) + coords = np.meshgrid(*axes, indexing="ij") + spacings = [axes[k][1] - axes[k][0] for k in range(len(shape))] + return axes, spacings, coords diff --git a/tests/test_operators.py b/tests/test_operators.py new file mode 100644 index 0000000..080b91d --- /dev/null +++ b/tests/test_operators.py @@ -0,0 +1,177 @@ +import numpy as np +from numpy.testing import assert_array_almost_equal + +from findiff.current import Add, FieldOperator +from findiff import Diff, Identity + + +def test_applying_identity(): + + f = np.array([1, 2, 3]) + I = Identity() + + assert np.array_equal(I(f), f) + + +def test_add_identities(): + + f = np.array([1, 2, 3]) + I = Identity() + + assert np.array_equal((I + I + I)(f), 3 * f) + + +def test_add_scalar_to_identity(): + + f = np.array([1, 2, 3]) + I = Identity() + + Ip1 = I + 2 + + assert isinstance(Ip1, Add) + assert isinstance(Ip1.left, Identity) + assert isinstance(Ip1.right, FieldOperator) + assert np.array_equal(Ip1(f), 3 * f) + + Ip1 = 2 + I + + assert isinstance(Ip1, Add) + assert isinstance(Ip1.left, Identity) + assert isinstance(Ip1.right, FieldOperator) + assert np.array_equal(Ip1(f), 3 * f) + + +def test_add_scalars(): + + f = np.array([1, 2, 3]) + S = FieldOperator(2) + + assert np.array_equal((S + S)(f), 4 * f) + + +def test_sub_scalars(): + + f = np.array([1, 2, 3]) + S = FieldOperator(2) + + actual = (S - 2 * S)(f) + expected = -2 * f + np.testing.assert_array_almost_equal(actual, expected) + + +def test_multiply_scalars(): + + f = np.array([1, 2, 3]) + S = FieldOperator(3) + + assert np.array_equal((S * S)(f), 9 * f) + + +def test_multiply_scalar_with_number(): + + f = np.array([1, 2, 3]) + S = FieldOperator(3) + + assert np.array_equal((S * 2)(f), 6 * f) + assert np.array_equal((2 * S)(f), 6 * f) + + +def test_multiply_scalar_with_field(): + + f = np.array([1, 2, 3]) + g = np.array([-1, 1, -1]) + S = FieldOperator(3) + + assert np.array_equal((S * g)(f), [-3, 6, -9]) + assert np.array_equal((g * S)(f), [-3, 6, -9]) + + +def test_mix_addition_and_multiplication(): + f = np.array([1, 2, 3]) + S2 = FieldOperator(2) + S3 = FieldOperator(3) + + assert np.array_equal((S2 + S3 * 4)(f), 2 * f + 3 * 4 * f) + + +def test_simple_derivatives(): + + x = np.linspace(0, 1, 100) + dx = x[1] - x[0] + f = x**3 + + d_dx = Diff(dx) + + actual = d_dx(f) + np.testing.assert_array_almost_equal(actual, 3 * x**2, decimal=3) + + actual = d_dx(f, acc=4) + np.testing.assert_array_almost_equal(actual, 3 * x**2) + + d_dx = Diff(dx, 1, 0) + + actual = d_dx(f, acc=4) + np.testing.assert_array_almost_equal(actual, 3 * x**2) + + d2_dx2 = Diff(dx, 2, 0) + + actual = d2_dx2(f, acc=4) + np.testing.assert_array_almost_equal(actual, 6 * x) + + +def test_chained_derivatives(): + + x = np.linspace(0, 1, 100) + dx = x[1] - x[0] + f = x**3 + + d_dx = Diff(dx) + + d2_dx2 = d_dx * d_dx + actual = d2_dx2(f, acc=4) + np.testing.assert_array_almost_equal(actual, 6 * x) + + +def test_harmonic(): + + x = np.linspace(0, 1, 100) + dx = x[1] - x[0] + f = x**3 + + T = -0.5 * Diff(dx, 2, 0) + V = 0.5 * x**2 + H = T + V + + actual = H(f, acc=4) + Tf = T(f, acc=4) + Vf = FieldOperator(V)(f) + np.testing.assert_array_almost_equal(Tf, -3 * x) + np.testing.assert_array_almost_equal(Vf, 0.5 * x**2 * f) + np.testing.assert_array_almost_equal(actual, -3 * x + 0.5 * x**2 * f) + + +def test_partial_diff(): + nx = 100 + x = np.linspace(0, np.pi, nx) + u = np.sin(x) + ux_ex = np.cos(x) + dx = x[1] - x[0] + fd = Diff(dx, 1, 0) + + ux = fd(u, spacing=dx, acc=4) + + assert_array_almost_equal(ux, ux_ex, decimal=5) + + ny = 80 + y = np.linspace(0, np.pi, ny) + dy = y[1] - y[0] + X, Y = np.meshgrid(x, y, indexing="ij") + + u = np.sin(X) * np.sin(Y) + uxy_ex = np.cos(X) * np.cos(Y) + + fd = Diff(dx, 1, 0) * Diff(dy, 1, 1) + + uxy = fd(u, acc=4) + + assert_array_almost_equal(uxy, uxy_ex, decimal=5) From 6ba937635f08ab8c76c0836b8488f6157351cf3e Mon Sep 17 00:00:00 2001 From: maroba Date: Wed, 11 Dec 2024 11:49:08 +0100 Subject: [PATCH 3/6] Almost done --- README.md | 91 +++++----- examples/examples-basic.ipynb | 17 +- findiff/__init__.py | 7 +- findiff/{legacy => }/coefs.py | 0 findiff/compatible.py | 34 ++++ findiff/current/__init__.py | 3 - findiff/current/operators.py | 122 -------------- findiff/grids.py | 29 ++++ findiff/legacy/__init__.py | 23 --- findiff/legacy/diff.py | 6 +- findiff/legacy/grids.py | 31 ---- findiff/legacy/operators.py | 102 ++++++------ findiff/operators.py | 204 +++++++++++++++++++++++ findiff/{legacy => }/pde.py | 0 findiff/{legacy => }/stencils.py | 0 findiff/{legacy => }/symbolic.py | 0 findiff/{legacy => }/utils.py | 0 findiff/{legacy => }/vector.py | 164 ++++++++++-------- tests/legacy/test_bugs.py | 9 +- tests/legacy/test_coefs.py | 11 +- tests/legacy/test_findiff.py | 8 +- tests/legacy/test_grids.py | 21 --- tests/legacy/test_pde.py | 2 +- tests/legacy/test_scaling.py | 2 +- tests/legacy/test_stencils.py | 4 +- tests/legacy/test_symbolic.py | 30 ++-- tests/legacy/test_utils.py | 7 +- tests/legacy/test_vector.py | 6 +- tests/test_diff.py | 69 +++++--- tests/test_operators.py | 107 +++++++++++- tests/test_pde.py | 277 +++++++++++++++++++++++++++++++ 31 files changed, 928 insertions(+), 458 deletions(-) rename findiff/{legacy => }/coefs.py (100%) create mode 100644 findiff/compatible.py delete mode 100644 findiff/current/__init__.py delete mode 100644 findiff/current/operators.py create mode 100644 findiff/grids.py delete mode 100644 findiff/legacy/grids.py create mode 100644 findiff/operators.py rename findiff/{legacy => }/pde.py (100%) rename findiff/{legacy => }/stencils.py (100%) rename findiff/{legacy => }/symbolic.py (100%) rename findiff/{legacy => }/utils.py (100%) rename findiff/{legacy => }/vector.py (70%) delete mode 100644 tests/legacy/test_grids.py create mode 100644 tests/test_pde.py diff --git a/README.md b/README.md index 402d776..7d165f3 100644 --- a/README.md +++ b/README.md @@ -32,76 +32,68 @@ You can find the documentation of the code including examples of application at ## Taking Derivatives -*findiff* allows to easily define derivative operators that you can apply to *numpy* arrays of any dimension. -The syntax for a simple derivative operator is +*findiff* allows to easily define derivative operators that you can apply to *numpy* arrays of +any dimension. -```python -FinDiff(axis, spacing, degree) -``` - -where `spacing` is the separation of grid points between neighboring grid points. Consider the 1D case +Consider the simple 1D case of a equidistant grid with a first derivative $\displaystyle \frac{\partial}{\partial x}$ along the only axis (0): ```python import numpy as np -from findiff import FinDiff +from findiff import Diff +# define the grid: x = np.linspace(0, 1, 100) -f = np.sin(x) # as an example +dx = x[1] - x[0] # the grid spacing -# time step dx -dx = x[1] - x[0] +# the array to differentiate: +f = np.sin(x) # as an example # Define the derivative: -d_dx = FinDiff(0, dx, 1) +d_dx = Diff(0, dx) # Apply it: df_dx = d_dx(f) ``` -Similary, you can define partial derivative operators along different axes or of higher degree, for example: - -| Math | *findiff* | | -| :---------------------------------------------------------: | ----------------------------------- | --------------------------------------- | -| $\displaystyle \frac{\partial}{\partial y}$ | ``FinDiff(1, dy, 1)`` | same as `` FinDiff(1, dy)`` | -| $\displaystyle \frac{\partial^4}{\partial y^4}$ | ``FinDiff(1, dy, 4)`` | any degree is possible | -| $\displaystyle \frac{\partial^3}{\partial x^2\partial z}$ | ``FinDiff((0, dx, 2), (2, dz, 1))`` | mixed also possible, one tuple per axis | -| $\displaystyle \frac{\partial}{\partial x_{10}}$ | ``FinDiff(10, dx10, 1)`` | number of axes not limited | +Similarly, you can define partial derivatives along other axes, for example, if $z$ is the 2-axis, we can write +$\frac{partial}{\partial z}$ as: -We can also take linear combinations of derivative operators, for example: - -$$ -2x \frac{\partial^3}{\partial x^2 \partial z} + 3 \sin(y)z^2 \frac{\partial^3}{\partial x \partial y^2} -$$ +```python +Diff(2, dz) +``` -is +`Diff` always creates a first derivative. For higher derivatives, you simply exponentiate them, for example for $\frac{\partial^2}{\partial_x^2}$ -```python -Coef(2*X) * FinDiff((0, dz, 2), (2, dz, 1)) + Coef(3*sin(Y)*Z**2) * FinDiff((0, dx, 1), (1, dy, 2)) +``` +d2_dx2 = Diff(0, dx)**2 ``` -where `X, Y, Z` are *numpy* arrays with meshed grid points. +and apply it as before. -Chaining differential operators is also possible, e.g. +You can also define more general differential operators intuitively, like $$ -\left( \frac{\partial}{\partial x} - \frac{\partial}{\partial y} \right) -\cdot \left( \frac{\partial}{\partial x} + \frac{\partial}{\partial y} \right) -= \frac{\partial^2}{\partial x^2} - \frac{\partial^2}{\partial y^2} +2x \frac{\partial^3}{\partial x^2 \partial z} + 3 \sin(y)z^2 \frac{\partial^3}{\partial x \partial y^2} $$ -can be written as + +which can be written as ```python -(FinDiff(0, dx) - FinDiff(1, dy)) * (FinDiff(0, dx) + FinDiff(1, dy)) -``` +# define the operator +diff_op = 2 * X * Diff(0)**2 * Diff(2) + 3 * sin(Y) * Z**2 * Diff(0) * Diff(1) -and +# set the grid you use (equidistant here) +diff_op.set_grid({0: dx, 1: dy, 2: dz}) -```python -FinDiff(0, dx, 2) - FinDiff(1, dy, 2) +# apply the operator +result = diff_op(f) ``` +where `X, Y, Z` are *numpy* arrays with meshed grid points. Here you see that you can also define your grid +lazily. + Of course, standard operators from vector calculus like gradient, divergence and curl are also available as shortcuts. @@ -113,10 +105,21 @@ When constructing an instance of `FinDiff`, you can request the desired accuracy order by setting the keyword argument `acc`. For example: ```python -d2_dx2 = findiff.FinDiff(0, dy, 2, acc=4) -d2f_dx2 = d2_dx2(f) +d_dx = Diff(0, dy, acc=4) +df_dx = d2_dx2(f) +``` + +Alternatively, you can also split operator definition and configuration: + +```python +d_dx = Diff(0, dx) +d_dx.set_accuracy(2) +df_dx = d2_dx2(f) ``` +which comes in handy if you have a complicated expression of differential operators, because then you +can specify it on the whole expression and it will be passed down to all basic operators. + If not specified, second order accuracy will be taken by default. ## Finite Difference Coefficients @@ -164,7 +167,7 @@ For a given _FinDiff_ differential operator, you can get the matrix representati using the `matrix(shape)` method, e.g. for a small 1D grid of 10 points: ```python -d2_dx2 = FinDiff(0, dx, 2) +d2_dx2 = Diff(0, dx)**2 mat = d2_dx2.matrix((10,)) # this method returns a scipy sparse matrix print(mat.toarray()) ``` @@ -276,13 +279,13 @@ u(0) = 0, \hspace{1em} u(10) = 1 $$ ```python -from findiff import FinDiff, Id, PDE +from findiff import Diff, Id, PDE shape = (300, ) t = numpy.linspace(0, 10, shape[0]) dt = t[1]-t[0] -L = FinDiff(0, dt, 2) - FinDiff(0, dt, 1) + Coef(5) * Id() +L = Diff(0, dt)**2 - Diff(0, dt) + 5 * Id() f = numpy.cos(2*t) bc = BoundaryConditions(shape) diff --git a/examples/examples-basic.ipynb b/examples/examples-basic.ipynb index 8974be5..82fa188 100644 --- a/examples/examples-basic.ipynb +++ b/examples/examples-basic.ipynb @@ -20,15 +20,18 @@ }, { "cell_type": "code", - "execution_count": 1, "metadata": { - "collapsed": false + "collapsed": false, + "jupyter": { + "is_executing": true + } }, - "outputs": [], "source": [ "import numpy as np\n", - "from findiff import FinDiff, coefficients, Coefficient" - ] + "from findiff import Diff" + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -79,9 +82,7 @@ "collapsed": true }, "outputs": [], - "source": [ - "d2_dx2 = FinDiff(0, dx, 2)" - ] + "source": "d2_dx2 = Diff(0, dx=dx, acc=2)" }, { "cell_type": "markdown", diff --git a/findiff/__init__.py b/findiff/__init__.py index 6d8763e..6eedddc 100644 --- a/findiff/__init__.py +++ b/findiff/__init__.py @@ -21,4 +21,9 @@ __version__ = "1.0.0.dev" -from .current import Diff, Identity +from .legacy import * +from .operators import Diff, Identity +from .pde import PDE, BoundaryConditions +from .compatible import Coef, Coefficient, FinDiff, Id +from .coefs import coefficients +from .vector import Gradient, Divergence, Curl, Laplacian diff --git a/findiff/legacy/coefs.py b/findiff/coefs.py similarity index 100% rename from findiff/legacy/coefs.py rename to findiff/coefs.py diff --git a/findiff/compatible.py b/findiff/compatible.py new file mode 100644 index 0000000..6026fed --- /dev/null +++ b/findiff/compatible.py @@ -0,0 +1,34 @@ +"""Provides an interface to obsolete classes for backward compatibility.""" + +from findiff import Diff +from findiff.operators import FieldOperator, Identity + + +def FinDiff(*args, **kwargs): + + if len(args) > 3: + raise ValueError("FinDiff accepts not more than 3 positional arguments.") + + def diff_from_tuple(tpl): + if len(tpl) == 3: + axis, h, order = tpl + return Diff(axis, h, **kwargs) ** order + elif len(tpl) == 2: + axis, h = tpl + return Diff(axis, h, **kwargs) + + if isinstance(args[0], (list, tuple)): + diffs = [] + for tpl in args: + diffs.append(diff_from_tuple(tpl)) + fd = diffs[0] + for diff in diffs[1:]: + fd = fd * diff + return fd + + return diff_from_tuple(args) + + +# Define aliasses for backward compatibility +Coefficient = Coef = FieldOperator +Id = Identity diff --git a/findiff/current/__init__.py b/findiff/current/__init__.py deleted file mode 100644 index 6f8596a..0000000 --- a/findiff/current/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -API_VERSION = 1 - -from .operators import * diff --git a/findiff/current/operators.py b/findiff/current/operators.py deleted file mode 100644 index f0ab4e7..0000000 --- a/findiff/current/operators.py +++ /dev/null @@ -1,122 +0,0 @@ -import numbers -from abc import ABC, abstractmethod - -import numpy as np -from scipy import sparse - - -class Node(ABC): - __array_priority__ = 100 # Make sure custom multiplication is called over numpy's - - @abstractmethod - def __call__(self, f, *args, **kwargs): - pass - - def __add__(self, other): - return Add(self, other) - - def __radd__(self, other): - return Add(self, other) - - def __sub__(self, other): - return Add(FieldOperator(-1) * other, self) - - def __rsub__(self, other): - return Add(FieldOperator(-1) * other, self) - - def __mul__(self, other): - return Mul(self, other) - - def __rmul__(self, other): - return Mul(other, self) - - -class FieldOperator(Node): - """An operator that multiplies an array pointwise.""" - - def __init__(self, value): - self.value = value - - def __call__(self, f, *args, **kwargs): - if isinstance(f, (numbers.Number, np.ndarray)): - return self.value * f - return self.value * super().__call__(f, *args, **kwargs) - - -class ScalarOperator(FieldOperator): - def __init__(self, value): - if not isinstance(value, numbers.Number): - raise ValueError("Expected number, got " + str(type(value))) - super().__init__(value) - - -class Identity(ScalarOperator): - def __init__(self): - super().__init__(1) - - -class Add(Node): - def __init__(self, left, right): - if isinstance(left, (numbers.Number, np.ndarray)): - left = FieldOperator(left) - if isinstance(right, (numbers.Number, np.ndarray)): - right = FieldOperator(right) - self.left = left - self.right = right - - def __call__(self, f, *args, **kwargs): - return self.left(f, *args, **kwargs) + self.right(f, *args, **kwargs) - - -class Mul(Node): - def __init__(self, left, right): - if isinstance(left, (numbers.Number, np.ndarray)): - left = FieldOperator(left) - if isinstance(right, (numbers.Number, np.ndarray)): - right = FieldOperator(right) - self.left = left - self.right = right - - def __call__(self, f, *args, **kwargs): - return self.left(self.right(f, *args, **kwargs), *args, **kwargs) - - -class Diff(Node): - - DEFAULT_ACC = 2 - - def __init__(self, spacing, order=1, axis=0, acc=DEFAULT_ACC): - if spacing <= 0: - raise ValueError("spacing must be greater than zero") - if order <= 0: - raise ValueError("order must be greater than zero") - - self.spacing = spacing - self.axis = axis - self.order = order - self.acc = acc - - self._fd = None - - def __call__(self, f, *args, **kwargs): - - if "acc" in kwargs: - # allow to pass down new accuracy deep in expression tree - new_acc = kwargs["acc"] - - if new_acc != self.acc: - self._fd = None - - self.acc = new_acc - - if self._fd is None: - self._build_differentiator() - - if isinstance(f, Node): - f = f(*args, **kwargs) - return self._fd(f) - - def _build_differentiator(self): - from findiff.legacy import FinDiff - - self._fd = FinDiff(self.axis, self.spacing, self.order, acc=self.acc) diff --git a/findiff/grids.py b/findiff/grids.py new file mode 100644 index 0000000..f47fabd --- /dev/null +++ b/findiff/grids.py @@ -0,0 +1,29 @@ +from abc import ABC + + +class Grid(ABC): + pass + + +class EquidistantGrid(Grid): + + def __init__(self, spacing: dict): + for axis, h in spacing.items(): + if h <= 0: + raise ValueError("spacing must be greater than zero") + + # key -> value == axis -> spacing for axis + self.spacing = spacing + + +class TensorProductGrid(Grid): + + def __init__(self, axes_coords: dict): + if len(axes_coords) == 0: + raise ValueError("axes_coords cannot be empty") + for axis, c in axes_coords.items(): + if len(c.shape) != 1: + raise ValueError( + "each item in axes_coords must have only one dimension" + ) + self.coords = axes_coords diff --git a/findiff/legacy/__init__.py b/findiff/legacy/__init__.py index 3647c29..e69de29 100644 --- a/findiff/legacy/__init__.py +++ b/findiff/legacy/__init__.py @@ -1,23 +0,0 @@ -from .coefs import coefficients -from .operators import FinDiff, Coef, Coefficient, Identity -from .pde import PDE, BoundaryConditions -from .symbolic import SymbolicMesh, SymbolicDiff -from .vector import Gradient, Divergence, Curl, Laplacian - -API_VERSION = 0 - -__all__ = [ - "coefficients", - "FinDiff", - "Identity", - "Coef", - "Coefficient", - "PDE", - "BoundaryConditions", - "SymbolicMesh", - "SymbolicDiff", - "Gradient", - "Divergence", - "Curl", - "Laplacian", -] diff --git a/findiff/legacy/diff.py b/findiff/legacy/diff.py index 79d5415..281a9a6 100644 --- a/findiff/legacy/diff.py +++ b/findiff/legacy/diff.py @@ -4,9 +4,9 @@ import scipy.sparse as sparse -from .coefs import coefficients, coefficients_non_uni -from .stencils import StencilSet -from .utils import * +from ..coefs import coefficients, coefficients_non_uni +from ..stencils import StencilSet +from ..utils import * DEFAULT_ACC = 2 diff --git a/findiff/legacy/grids.py b/findiff/legacy/grids.py deleted file mode 100644 index a43baab..0000000 --- a/findiff/legacy/grids.py +++ /dev/null @@ -1,31 +0,0 @@ -import numpy as np - - -class Grid(object): - pass - - -class UniformGrid(Grid): - - def __init__(self, shape, spac, center=None): - - if not hasattr(shape, '__len__'): - self.shape = shape, - self.ndims = 1 - else: - self.shape = shape - self.ndims = len(shape) - - if not hasattr(spac, '__len__'): - self.spac = spac, - else: - self.spac = spac - - if center is None: - self.center = np.zeros(self.ndims) - else: - assert len(center) == self.ndims - self.center = np.array(center) - - def spacing(self, axis): - return self.spac[axis] diff --git a/findiff/legacy/operators.py b/findiff/legacy/operators.py index 372ea6c..00e87b9 100644 --- a/findiff/legacy/operators.py +++ b/findiff/legacy/operators.py @@ -1,75 +1,74 @@ -from .utils import * +from ..utils import * from .diff import Coef, Id, Diff, Plus, Minus, Mul, DEFAULT_ACC, LinearMap -from .stencils import StencilSet +from ..stencils import StencilSet class FinDiff(LinearMap): - r""" A representation of a general linear differential operator expressed in finite differences. + r"""A representation of a general linear differential operator expressed in finite differences. - FinDiff objects can be added with other FinDiff objects. They can be multiplied by - objects of type Coefficient. + FinDiff objects can be added with other FinDiff objects. They can be multiplied by + objects of type Coefficient. - FinDiff is callable, i.e. to apply the derivative, just call the object on the array to - differentiate. + FinDiff is callable, i.e. to apply the derivative, just call the object on the array to + differentiate. - :param args: variable number of tuples. Defines what derivative to take. - If only one tuple is given, you can leave away the tuple parentheses. + :param args: variable number of tuples. Defines what derivative to take. + If only one tuple is given, you can leave away the tuple parentheses. - Each tuple has the form + Each tuple has the form - `(axis, spacing, count)` for uniform grids + `(axis, spacing, count)` for uniform grids - `(axis, count)` for non-uniform grids. + `(axis, count)` for non-uniform grids. - `axis` is the dimension along which to take derivative. + `axis` is the dimension along which to take derivative. - `spacing` is the grid spacing of the uniform grid along that axis. + `spacing` is the grid spacing of the uniform grid along that axis. - `count` is the order of the derivative, which is optional an defaults to 1. + `count` is the order of the derivative, which is optional an defaults to 1. - :param kwargs: variable number of keyword arguments + :param kwargs: variable number of keyword arguments - Allowed keywords: + Allowed keywords: - `acc`: even integer - The desired accuracy order. Default is acc=2. + `acc`: even integer + The desired accuracy order. Default is acc=2. - This class is actually deprecated and will be replaced by the Diff class in the future. + This class is actually deprecated and will be replaced by the Diff class in the future. - **Example**: + **Example**: - For this example, we want to operate on some 3D array f: + For this example, we want to operate on some 3D array f: - >>> import numpy as np - >>> x, y, z = [np.linspace(-1, 1, 100) for _ in range(3)] - >>> X, Y, Z = np.meshgrid(x, y, z, indexing='ij') - >>> f = X**2 + Y**2 + Z**2 + >>> import numpy as np + >>> x, y, z = [np.linspace(-1, 1, 100) for _ in range(3)] + >>> X, Y, Z = np.meshgrid(x, y, z, indexing='ij') + >>> f = X**2 + Y**2 + Z**2 - To create :math:`\\frac{\\partial f}{\\partial x}` on a uniform grid with spacing dx, dy - along the 0th axis or 1st axis, respectively, instantiate a FinDiff object and call it: + To create :math:`\\frac{\\partial f}{\\partial x}` on a uniform grid with spacing dx, dy + along the 0th axis or 1st axis, respectively, instantiate a FinDiff object and call it: - >>> d_dx = FinDiff(0, dx) - >>> d_dy = FinDiff(1, dx) - >>> result = d_dx(f) + >>> d_dx = FinDiff(0, dx) + >>> d_dy = FinDiff(1, dx) + >>> result = d_dx(f) - For :math:`\\frac{\\partial^2 f}{\\partial x^2}` or :math:`\\frac{\\partial^2 f}{\\partial y^2}`: + For :math:`\\frac{\\partial^2 f}{\\partial x^2}` or :math:`\\frac{\\partial^2 f}{\\partial y^2}`: - >>> d2_dx2 = FinDiff(0, dx, 2) - >>> d2_dy2 = FinDiff(1, dy, 2) - >>> result_2 = d2_dx2(f) - >>> result_3 = d2_dy2(f) + >>> d2_dx2 = FinDiff(0, dx, 2) + >>> d2_dy2 = FinDiff(1, dy, 2) + >>> result_2 = d2_dx2(f) + >>> result_3 = d2_dy2(f) - For :math:`\\frac{\\partial^4 f}{\partial x \\partial^2 y \\partial z}`, do: + For :math:`\\frac{\\partial^4 f}{\partial x \\partial^2 y \\partial z}`, do: - >>> op = FinDiff((0, dx), (1, dy, 2), (2, dz)) - >>> result_4 = op(f) + >>> op = FinDiff((0, dx), (1, dy, 2), (2, dz)) + >>> result_4 = op(f) """ - def __init__(self, *args, **kwargs): self.acc = None self.spac = None @@ -79,18 +78,18 @@ def __call__(self, rhs, *args, **kwargs): return self.apply(rhs, *args, **kwargs) def apply(self, rhs, *args, **kwargs): - if 'acc' not in kwargs: + if "acc" not in kwargs: if self.acc is None: acc = DEFAULT_ACC else: acc = self.acc - kwargs['acc'] = acc + kwargs["acc"] = acc - if len(args) == 0 and 'h' not in kwargs: + if len(args) == 0 and "h" not in kwargs: if self.uniform: - args = self.spac, + args = (self.spac,) else: - args = self.coords, + args = (self.coords,) return self.pds(rhs, *args, **kwargs) def stencil(self, shape): @@ -125,10 +124,10 @@ def __mul__(self, other): def _eval_args(self, args, kwargs): spac = {} - if 'acc' in kwargs: - self.acc = kwargs['acc'] + if "acc" in kwargs: + self.acc = kwargs["acc"] - if isinstance(args[0], tuple): # mixed partial derivative + if isinstance(args[0], tuple): # mixed partial derivative pds = None for arg in args: @@ -138,7 +137,7 @@ def _eval_args(self, args, kwargs): axis, h = arg order = 1 else: - raise ValueError('Format: (axis, spacing, order=1)') + raise ValueError("Format: (axis, spacing, order=1)") spac[axis] = h if pds is None: pds = Diff(axis, order) @@ -152,14 +151,14 @@ def _eval_args(self, args, kwargs): axis, h = args order = 1 else: - raise ValueError('Format: (axis, spacing, order=1)') + raise ValueError("Format: (axis, spacing, order=1)") pds = Diff(axis, order) spac[axis] = h # Check if spac is really the spacing and not the coordinates (nonuniform case) for a, s in spac.items(): - if hasattr(s, '__len__'): + if hasattr(s, "__len__"): self.coords = spac self.uniform = False break @@ -174,4 +173,3 @@ def _eval_args(self, args, kwargs): # Alias for backward compatibility Coefficient = Coef Identity = Id - diff --git a/findiff/operators.py b/findiff/operators.py new file mode 100644 index 0000000..5f10ed6 --- /dev/null +++ b/findiff/operators.py @@ -0,0 +1,204 @@ +import numbers +from abc import ABC, abstractmethod + +import numpy as np +from scipy import sparse + +from findiff.grids import EquidistantGrid, Grid, TensorProductGrid +from findiff.stencils import StencilSet + + +class Node(ABC): + __array_priority__ = 100 # Makes sure custom multiplication is called over numpy's + + def __init__(self, *args, **kwargs): + self.children = [] + + @abstractmethod + def __call__(self, f, *args, **kwargs): + pass + + @abstractmethod + def matrix(self, shape): + pass + + def stencil(self, shape): + return StencilSet(self, shape) + + def __add__(self, other): + return Add(self, other) + + def __radd__(self, other): + return Add(self, other) + + def __sub__(self, other): + return Add(ScalarOperator(-1) * other, self) + + def __rsub__(self, other): + return Add(ScalarOperator(-1) * other, self) + + def __mul__(self, other): + return Mul(self, other) + + def __rmul__(self, other): + return Mul(other, self) + + @property + def grid(self): + return getattr(self, "_grid", None) + + def set_grid(self, grid): + if isinstance(grid, dict): + self._grid = EquidistantGrid(grid) + for child in self.children: + child.set_grid(self._grid) + else: + # Grid may be set lazily + self._grid = grid + for child in self.children: + child.set_grid(self._grid) + + def set_accuracy(self, acc): + self.acc = acc + for child in self.children: + child.set_accuracy(acc) + + +class FieldOperator(Node): + """An operator that multiplies an array pointwise.""" + + def __init__(self, value): + super().__init__() + self.value = value + + def __call__(self, f, *args, **kwargs): + if isinstance(f, (numbers.Number, np.ndarray)): + return self.value * f + return self.value * super().__call__(f, *args, **kwargs) + + def matrix(self, shape): + if isinstance(self.value, np.ndarray): + diag_values = self.value.reshape(-1) + return sparse.diags(diag_values) + elif isinstance(self.value, numbers.Number): + siz = np.prod(shape) + return sparse.diags(self.value * np.ones(siz)) + + +class ScalarOperator(FieldOperator): + def __init__(self, value): + if not isinstance(value, numbers.Number): + raise ValueError("Expected number, got " + str(type(value))) + super().__init__(value) + + def matrix(self, shape): + siz = np.prod(shape) + mat = sparse.lil_matrix((siz, siz)) + diag = list(range(siz)) + mat[diag, diag] = self.value + return sparse.csr_matrix(mat) + + +class Identity(ScalarOperator): + def __init__(self): + super().__init__(1) + + +class BinaryOperation(Node): + + @property + def left(self): + return self.children[0] + + @property + def right(self): + return self.children[1] + + +class Add(BinaryOperation): + def __init__(self, left, right): + if isinstance(left, (numbers.Number, np.ndarray)): + left = FieldOperator(left) + if isinstance(right, (numbers.Number, np.ndarray)): + right = FieldOperator(right) + super().__init__() + self.children = [left, right] + + def __call__(self, f, *args, **kwargs): + return self.left(f, *args, **kwargs) + self.right(f, *args, **kwargs) + + def matrix(self, shape): + return self.left.matrix(shape) + self.right.matrix(shape) + + +class Mul(BinaryOperation): + def __init__(self, left, right): + if isinstance(left, (numbers.Number, np.ndarray)): + left = FieldOperator(left) + if isinstance(right, (numbers.Number, np.ndarray)): + right = FieldOperator(right) + super().__init__() + self.children = [left, right] + + def __call__(self, f, *args, **kwargs): + return self.left(self.right(f, *args, **kwargs), *args, **kwargs) + + def matrix(self, shape): + return self.left.matrix(shape) * self.right.matrix(shape) + + +class Diff(Node): + + DEFAULT_ACC = 2 + + def __init__(self, axis=0, grid=None, acc=DEFAULT_ACC): + super().__init__() + + if isinstance(grid, numbers.Number): + grid = {axis: grid} + elif hasattr(grid, "shape") and hasattr(grid, "__len__"): + grid = TensorProductGrid({axis: grid}) + + self.set_grid(grid) + + self.axis = axis + self.order = 1 + self.acc = acc + + self._fd = None + + def __call__(self, f, *args, **kwargs): + + if "acc" in kwargs: + # allow to pass down new accuracy deep in expression tree + new_acc = kwargs["acc"] + if new_acc != self.acc: + self._fd = None + self.set_accuracy(new_acc) + + if self._fd is None: + self._build_differentiator() + + if isinstance(f, Node): + f = f(*args, **kwargs) + return self._fd(f) + + def _build_differentiator(self): + from findiff.legacy.operators import FinDiff + + if isinstance(self.grid, EquidistantGrid): + spacing = self.grid.spacing[self.axis] + self._fd = FinDiff(self.axis, spacing, self.order, acc=self.acc) + elif isinstance(self.grid, TensorProductGrid): + coords = self.grid.coords[self.axis] + self._fd = FinDiff(self.axis, coords, self.order, acc=self.acc) + + def matrix(self, shape): + if not self._fd: + self._build_differentiator() + return self._fd.matrix(shape) + + def __pow__(self, power): + new_diff = Diff(self.axis, acc=self.acc, grid=self.grid) + new_diff.order *= power + return new_diff diff --git a/findiff/legacy/pde.py b/findiff/pde.py similarity index 100% rename from findiff/legacy/pde.py rename to findiff/pde.py diff --git a/findiff/legacy/stencils.py b/findiff/stencils.py similarity index 100% rename from findiff/legacy/stencils.py rename to findiff/stencils.py diff --git a/findiff/legacy/symbolic.py b/findiff/symbolic.py similarity index 100% rename from findiff/legacy/symbolic.py rename to findiff/symbolic.py diff --git a/findiff/legacy/utils.py b/findiff/utils.py similarity index 100% rename from findiff/legacy/utils.py rename to findiff/utils.py diff --git a/findiff/legacy/vector.py b/findiff/vector.py similarity index 70% rename from findiff/legacy/vector.py rename to findiff/vector.py index b81669b..6b54523 100644 --- a/findiff/legacy/vector.py +++ b/findiff/vector.py @@ -1,28 +1,28 @@ """A module for the common differential operators of vector calculus""" import numpy as np -from .operators import FinDiff -from.diff import Diff + +from .legacy.operators import FinDiff class VectorOperator(object): """Base class for all vector differential operators. - Shall not be instantiated directly, but through the child classes. + Shall not be instantiated directly, but through the child classes. """ def __init__(self, **kwargs): """Constructor for the VectorOperator base class. - - kwargs: - ------- - - h list with the grid spacings of an N-dimensional uniform grid - - coords list of 1D arrays with the coordinate values along the N axes. - This is used for non-uniform grids. - - Either specify "h" or "coords", not both. - + + kwargs: + ------- + + h list with the grid spacings of an N-dimensional uniform grid + + coords list of 1D arrays with the coordinate values along the N axes. + This is used for non-uniform grids. + + Either specify "h" or "coords", not both. + """ if "acc" in kwargs: @@ -30,7 +30,9 @@ def __init__(self, **kwargs): else: self.acc = 2 - if "spac" in kwargs or "h" in kwargs: # necessary for backward compatibility 0.5.2 => 0.6 + if ( + "spac" in kwargs or "h" in kwargs + ): # necessary for backward compatibility 0.5.2 => 0.6 if "spac" in kwargs: kw = "spac" else: @@ -42,7 +44,9 @@ def __init__(self, **kwargs): if "coords" in kwargs: coords = kwargs.pop("coords") self.ndims = self.__get_dimension(coords) - self.components = [FinDiff((k, coords[k], 1), **kwargs) for k in range(self.ndims)] + self.components = [ + FinDiff((k, coords[k], 1), **kwargs) for k in range(self.ndims) + ] def __get_dimension(self, coords): return len(coords) @@ -51,18 +55,18 @@ def __get_dimension(self, coords): class Gradient(VectorOperator): r""" The N-dimensional gradient. - + .. math:: \nabla = \left(\frac{\partial}{\partial x_0}, \frac{\partial}{\partial x_1}, ... , \frac{\partial}{\partial x_{N-1}}\right) :param kwargs: exactly one of *h* and *coords* must be specified - - *h* - list with the grid spacings of an N-dimensional uniform grid + + *h* + list with the grid spacings of an N-dimensional uniform grid *coords* list of 1D arrays with the coordinate values along the N axes. This is used for non-uniform grids. - + *acc* accuracy order, must be positive integer, default is 2 """ @@ -73,17 +77,17 @@ def __init__(self, **kwargs): def __call__(self, f): """ Applies the N-dimensional gradient to the array f. - + :param f: ``numpy.ndarray`` - + Array to apply the gradient to. It represents a scalar function, - so it must have N axes for the N independent variables. - + so it must have N axes for the N independent variables. + :returns: ``numpy.ndarray`` - - The gradient of f, which has N+1 axes, i.e. it is + + The gradient of f, which has N+1 axes, i.e. it is an array of N arrays of N axes each. - + """ if not isinstance(f, np.ndarray): @@ -103,22 +107,22 @@ def __call__(self, f): class Divergence(VectorOperator): r""" The N-dimensional divergence. - + .. math:: - + {\rm \bf div} = \nabla \cdot - + :param kwargs: exactly one of *h* and *coords* must be specified - *h* - list with the grid spacings of an N-dimensional uniform grid + *h* + list with the grid spacings of an N-dimensional uniform grid *coords* list of 1D arrays with the coordinate values along the N axes. This is used for non-uniform grids. - + *acc* accuracy order, must be positive integer, default is 2 - + """ def __init__(self, **kwargs): @@ -129,19 +133,23 @@ def __call__(self, f): Applies the divergence to the array f. :param f: ``numpy.ndarray`` - + a vector function of N variables, so its array has N+1 axes. - + :returns: ``numpy.ndarray`` - + the divergence, which is a scalar function of N variables, so it's array dimension has N axes - + """ if not isinstance(f, np.ndarray) and not isinstance(f, list): - raise TypeError("Function to differentiate must be numpy.ndarray or list of numpy.ndarrays") + raise TypeError( + "Function to differentiate must be numpy.ndarray or list of numpy.ndarrays" + ) if len(f.shape) != self.ndims + 1 and f.shape[0] != self.ndims: - raise ValueError("Divergence can only be applied to vector functions of the same dimension") + raise ValueError( + "Divergence can only be applied to vector functions of the same dimension" + ) result = np.zeros(f.shape[1:]) @@ -153,33 +161,37 @@ def __call__(self, f): class Curl(VectorOperator): r""" - The curl operator. - + The curl operator. + .. math:: - + {\rm \bf rot} = \nabla \times - + Is only defined for 3D. - + :param kwargs: exactly one of *h* and *coords* must be specified - *h* - list with the grid spacings of a 3-dimensional uniform grid + *h* + list with the grid spacings of a 3-dimensional uniform grid *coords* list of 1D arrays with the coordinate values along the 3 axes. This is used for non-uniform grids. - + *acc* accuracy order, must be positive integer, default is 2 - + """ def __init__(self, **kwargs): super().__init__(**kwargs) if self.ndims != 3: - raise ValueError("Curl operation is only defined in 3 dimensions. {} were given.".format(self.ndims)) + raise ValueError( + "Curl operation is only defined in 3 dimensions. {} were given.".format( + self.ndims + ) + ) def __call__(self, f): """ @@ -196,44 +208,54 @@ def __call__(self, f): """ if not isinstance(f, np.ndarray) and not isinstance(f, list): - raise TypeError("Function to differentiate must be numpy.ndarray or list of numpy.ndarrays") + raise TypeError( + "Function to differentiate must be numpy.ndarray or list of numpy.ndarrays" + ) if len(f.shape) != self.ndims + 1 and f.shape[0] != self.ndims: - raise ValueError("Curl can only be applied to vector functions of the three dimensions") + raise ValueError( + "Curl can only be applied to vector functions of the three dimensions" + ) result = np.zeros(f.shape) - result[0] += self.components[1](f[2], acc=self.acc) - self.components[2](f[1], acc=self.acc) - result[1] += self.components[2](f[0], acc=self.acc) - self.components[0](f[2], acc=self.acc) - result[2] += self.components[0](f[1], acc=self.acc) - self.components[1](f[0], acc=self.acc) + result[0] += self.components[1](f[2], acc=self.acc) - self.components[2]( + f[1], acc=self.acc + ) + result[1] += self.components[2](f[0], acc=self.acc) - self.components[0]( + f[2], acc=self.acc + ) + result[2] += self.components[0](f[1], acc=self.acc) - self.components[1]( + f[0], acc=self.acc + ) return result class Laplacian(object): r""" - The N-dimensional Laplace operator. + The N-dimensional Laplace operator. - .. math:: + .. math:: - {\rm \bf \nabla^2} = \sum_{k=0}^{N-1} \frac{\partial^2}{\partial x_k^2} + {\rm \bf \nabla^2} = \sum_{k=0}^{N-1} \frac{\partial^2}{\partial x_k^2} - :param kwargs: exactly one of *h* and *coords* must be specified + :param kwargs: exactly one of *h* and *coords* must be specified - *h* - list with the grid spacings of an N-dimensional uniform grid - *coords* - list of 1D arrays with the coordinate values along the N axes. - This is used for non-uniform grids. + *h* + list with the grid spacings of an N-dimensional uniform grid + *coords* + list of 1D arrays with the coordinate values along the N axes. + This is used for non-uniform grids. - *acc* - accuracy order, must be positive integer, default is 2 + *acc* + accuracy order, must be positive integer, default is 2 - """ + """ """A representation of the Laplace operator in arbitrary dimensions using finite difference schemes""" - def __init__(self, h=[1.], acc=2): + def __init__(self, h=[1.0], acc=2): h = wrap_in_ndarray(h) self._parts = [FinDiff((k, h[k], 2), acc=acc) for k in range(len(h))] @@ -262,8 +284,8 @@ def __call__(self, f): def wrap_in_ndarray(value): """Wraps the argument in a numpy.ndarray. - If value is a scalar, it is converted in a list first. - If value is array-like, the shape is conserved. + If value is a scalar, it is converted in a list first. + If value is array-like, the shape is conserved. """ diff --git a/tests/legacy/test_bugs.py b/tests/legacy/test_bugs.py index d4a72c3..7dbd272 100644 --- a/tests/legacy/test_bugs.py +++ b/tests/legacy/test_bugs.py @@ -1,12 +1,9 @@ -import sys - -sys.path.insert(1, "..") - import unittest + import numpy as np -from findiff.legacy import FinDiff -import findiff.legacy as findiff +import findiff +from findiff import FinDiff class TestOldBugs(unittest.TestCase): diff --git a/tests/legacy/test_coefs.py b/tests/legacy/test_coefs.py index 12af237..61e497d 100644 --- a/tests/legacy/test_coefs.py +++ b/tests/legacy/test_coefs.py @@ -1,15 +1,12 @@ -import sys - -sys.path.insert(1, "..") - import unittest -from findiff.legacy import coefficients -from findiff.legacy.coefs import coefficients_non_uni -from findiff.legacy.coefs import calc_coefs import numpy as np from sympy import Rational +from findiff import coefficients +from findiff.coefs import calc_coefs +from findiff.coefs import coefficients_non_uni + class TestCoefs(unittest.TestCase): diff --git a/tests/legacy/test_findiff.py b/tests/legacy/test_findiff.py index 5bf51c3..0b6858a 100644 --- a/tests/legacy/test_findiff.py +++ b/tests/legacy/test_findiff.py @@ -1,11 +1,9 @@ -import sys - -sys.path.insert(1, "..") - import unittest + import numpy as np from numpy.testing import assert_array_almost_equal, assert_array_equal -from findiff.legacy.operators import FinDiff, Coef, Identity + +from findiff import FinDiff, Coef, Identity class FinDiffTest(unittest.TestCase): diff --git a/tests/legacy/test_grids.py b/tests/legacy/test_grids.py deleted file mode 100644 index de07530..0000000 --- a/tests/legacy/test_grids.py +++ /dev/null @@ -1,21 +0,0 @@ -import unittest - -import numpy as np - -from findiff.legacy.grids import UniformGrid - - -class TestUniformGrid(unittest.TestCase): - - def test_init_2d(self): - shape = 30, 30 - spac = 0.1, 0.2 - center = 2, 3 - grid = UniformGrid(shape, spac, center) - self.assertEqual(0.1, grid.spacing(0)) - np.testing.assert_array_equal(center, grid.center) - - def test_init_1d(self): - grid = UniformGrid(30, 0.1) - self.assertEqual(0.1, grid.spacing(0)) - np.testing.assert_array_equal((0,), grid.center) diff --git a/tests/legacy/test_pde.py b/tests/legacy/test_pde.py index 0f8392d..99da46e 100644 --- a/tests/legacy/test_pde.py +++ b/tests/legacy/test_pde.py @@ -4,7 +4,7 @@ import unittest from findiff.legacy.operators import FinDiff, Identity, Coef -from findiff.legacy.pde import * +from findiff.pde import * # import matplotlib.pyplot as plt # from mpl_toolkits import mplot3d diff --git a/tests/legacy/test_scaling.py b/tests/legacy/test_scaling.py index 6605e4e..8fcead0 100644 --- a/tests/legacy/test_scaling.py +++ b/tests/legacy/test_scaling.py @@ -5,7 +5,7 @@ from math import log import unittest import numpy as np -from findiff.legacy.operators import FinDiff +from findiff import FinDiff class TestScaling(unittest.TestCase): diff --git a/tests/legacy/test_stencils.py b/tests/legacy/test_stencils.py index cf89727..ca40f62 100644 --- a/tests/legacy/test_stencils.py +++ b/tests/legacy/test_stencils.py @@ -2,9 +2,9 @@ import numpy as np -from findiff.legacy import Identity, FinDiff +from findiff import Identity, FinDiff -from findiff.legacy.stencils import Stencil +from findiff.stencils import Stencil class TestStencilOperations(unittest.TestCase): diff --git a/tests/legacy/test_symbolic.py b/tests/legacy/test_symbolic.py index ea4b30e..95d15f7 100644 --- a/tests/legacy/test_symbolic.py +++ b/tests/legacy/test_symbolic.py @@ -2,7 +2,7 @@ from sympy import IndexedBase, Symbol, symbols, latex -from findiff.legacy.symbolic import SymbolicMesh, SymbolicDiff +from findiff.symbolic import SymbolicMesh, SymbolicDiff class TestSymbolicMesh(unittest.TestCase): @@ -10,8 +10,8 @@ class TestSymbolicMesh(unittest.TestCase): def test_parse_symbolic_mesh(self): # 1D mesh = SymbolicMesh(coord="x", equidistant=True) - x, = mesh.coord - dx, = mesh.spacing + (x,) = mesh.coord + (dx,) = mesh.spacing self.assertEqual(IndexedBase, type(x)) self.assertEqual(Symbol, type(dx)) @@ -44,17 +44,17 @@ def test_create_symbol(self): mesh = SymbolicMesh(coord="x,y", equidistant=True) n, m = symbols("n, m") actual = latex(mesh.create_symbol("u")[n, m]) - #expected = "u{}_{n}{}_{m}" + # expected = "u{}_{n}{}_{m}" expected = "{u}_{n,m}" self.assertEqual(latex(actual), latex(expected)) # both indices up - #mesh = SymbolicMesh(coord="x,y", equidistant=True) - #n, m = symbols("n, m") - #u = mesh.create_symbol("u", pos="u,u") - #actual = latex(u[n, m]) - #expected = "u{}^{n}{}^{m}" - #self.assertEqual(actual, expected) + # mesh = SymbolicMesh(coord="x,y", equidistant=True) + # n, m = symbols("n, m") + # u = mesh.create_symbol("u", pos="u,u") + # actual = latex(u[n, m]) + # expected = "u{}^{n}{}^{m}" + # self.assertEqual(actual, expected) class TestDiff(unittest.TestCase): @@ -78,9 +78,7 @@ def test_call(self): expected = (u[n + 1] - u[n - 1]) / (2 * mesh.spacing[0]) - self.assertEqual( - 0, (expected - actual).simplify() - ) + self.assertEqual(0, (expected - actual).simplify()) # 2D mesh = SymbolicMesh("x, y") @@ -92,10 +90,8 @@ def test_call(self): expected = (u[n, m + 1] - u[n, m - 1]) / (2 * mesh.spacing[1]) - self.assertEqual( - 0, (expected - actual).simplify() - ) + self.assertEqual(0, (expected - actual).simplify()) -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/tests/legacy/test_utils.py b/tests/legacy/test_utils.py index 018f613..42ad9b9 100644 --- a/tests/legacy/test_utils.py +++ b/tests/legacy/test_utils.py @@ -1,9 +1,6 @@ -import sys - -sys.path.insert(1, "..") - import unittest -from findiff.legacy.utils import * + +from findiff.utils import * class TestUtils(unittest.TestCase): diff --git a/tests/legacy/test_vector.py b/tests/legacy/test_vector.py index b135dfd..045d49b 100644 --- a/tests/legacy/test_vector.py +++ b/tests/legacy/test_vector.py @@ -1,11 +1,7 @@ -import sys - -sys.path.insert(1, "..") - import unittest from numpy.testing import assert_array_almost_equal import numpy as np -from findiff.legacy.vector import Gradient, Divergence, Curl, Laplacian +from findiff import Gradient, Divergence, Curl, Laplacian class TestGradient(unittest.TestCase): diff --git a/tests/test_diff.py b/tests/test_diff.py index 782e00e..1142707 100644 --- a/tests/test_diff.py +++ b/tests/test_diff.py @@ -14,7 +14,7 @@ def test_partial_d_dx(): u = x**2 expected = 2 * x - fd = Diff(dx, 1, 0) + fd = Diff(0, dx) actual = fd(u) assert_array_almost_equal(expected, actual) @@ -27,11 +27,8 @@ def test_partial_d2_dx2(): u = x**2 expected = 2 - fd = Diff( - dx, - 2, - 0, - ) + fd = Diff(0, dx) ** 2 + actual = fd(u, dx) assert_array_almost_equal(expected, actual) @@ -44,8 +41,8 @@ def test_partial_d_dx_acc(): u = x**3 expected = 3 * x**2 - fd = Diff(dx, 1, 0) - actual = fd(u, dx) + fd = Diff(0, dx) + actual = fd(u) np.testing.assert_raises( AssertionError, assert_array_almost_equal, expected, actual ) @@ -60,16 +57,13 @@ def test_partial_d2_dxdy(): (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) u = np.sin(X) * np.sin(Y) - fd = Diff(dx, 1, 0) * Diff(dy, 1, 1) - actual = fd(u, h={0: dx, 1: dy}) + fd = Diff(0) * Diff(1) + fd.set_grid({0: dx, 1: dy}) + actual = fd(u) expected = np.cos(X) * np.cos(Y) assert_array_almost_equal(expected, actual, decimal=3) - fd = Diff(dy, 1, 1) * Diff(dy, 1, 0) - actual = fd(u) - assert_array_almost_equal(expected, actual, decimal=3) - def test_partial_d_dx_on_2d_array(): @@ -77,7 +71,7 @@ def test_partial_d_dx_on_2d_array(): (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) u = np.sin(X) * np.sin(Y) - fd = Diff(dy, 1, 1) + fd = Diff(1, dy) actual = fd(u) expected = np.sin(X) * np.cos(Y) @@ -90,7 +84,7 @@ def test_add_two_diffs(): (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) u = np.sin(X) * np.sin(Y) - fd = Diff(dx, 2, 0) + Diff(dy, 2, 1) + fd = Diff(0, dx) ** 2 + Diff(1, dy) ** 2 actual = fd(u) expected = -2 * np.sin(X) * np.sin(Y) @@ -103,7 +97,7 @@ def test_sub_two_diffs(): (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) u = np.sin(X) * np.sin(Y) - fd = Diff(dx, 2, 0) - Diff(dy, 2, 1) + fd = Diff(0, dx) ** 2 - Diff(1, dy) ** 2 actual = fd(u) expected = np.zeros(shape) @@ -116,7 +110,7 @@ def test_multiply_with_scalar(): (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) u = np.sin(X) * np.sin(Y) - fd = 3 * Diff(dx, 2, 0) + fd = 3 * Diff(0, dx) ** 2 actual = fd(u, acc=4) expected = -3 * np.sin(X) * np.sin(Y) @@ -129,7 +123,7 @@ def test_multiply_with_variable_from_left(): (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) u = np.sin(X) * np.sin(Y) - fd = 3 * X * Y * Diff(dx, 2, 0) + fd = 3 * X * Y * Diff(0, dx) ** 2 actual = fd(u, acc=4) expected = -3 * X * Y * np.sin(X) * np.sin(Y) @@ -142,7 +136,7 @@ def test_identity(): (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) u = np.sin(X) * np.sin(Y) - fd = Diff(dx, 2, 0) + Identity() + fd = Diff(0, dx) ** 2 + Identity() actual = fd(u, acc=4) expected = np.zeros_like(u) @@ -154,7 +148,7 @@ def test_linear_comb(): (x, y), (dx, dy), (X, Y) = make_grid(shape, ((0, 1), (0, 1))) u = np.sin(X) * np.sin(Y) - fd = 3 * X * Diff(dx, 2, 0) - Y * (Diff(dy, 2, 1) + Diff(dx, 1, 0)) + fd = 3 * X * Diff(0, dx) ** 2 - Y * (Diff(1, dy) ** 2 + Diff(0, dx)) actual = fd(u, acc=4) expected = -3 * X * u - Y * (-sin(X) * sin(Y) + cos(X) * sin(Y)) @@ -174,3 +168,36 @@ def make_grid(shape, edges): coords = np.meshgrid(*axes, indexing="ij") spacings = [axes[k][1] - axes[k][0] for k in range(len(shape))] return axes, spacings, coords + + +def test_nonuniform_1d_different_accs(): + x = np.r_[np.arange(0, 4, 0.05), np.arange(4, 10, 1)] + f = np.exp(-(x**2)) + + d_dx = Diff(0, x, acc=4) + f_x = d_dx(f) + assert_array_almost_equal(-2 * x * np.exp(-(x**2)), f_x, decimal=4) + + # same, but this time with default acc + x = np.linspace(0, 1, 100) + f = np.exp(-(x**2)) + d_dx = Diff(0, x) + f_x = d_dx(f) + assert_array_almost_equal(-2 * x * np.exp(-(x**2)), f_x, decimal=4) + + +def test_3d_different_accs(): + x = np.r_[np.arange(0, 4, 0.05), np.arange(4, 10, 1)] + y = np.r_[np.arange(0, 4, 0.05), np.arange(4, 10, 1)] + z = np.r_[np.arange(0, 4.5, 0.05), np.arange(4.5, 10, 1)] + X, Y, Z = np.meshgrid(x, y, z, indexing="ij") + f = np.exp(-(X**2) - Y**2 - Z**2) + + d_dy = Diff(1, y, acc=4) + fy = d_dy(f) + fye = -2 * Y * np.exp(-(X**2) - Y**2 - Z**2) + assert_array_almost_equal(fy, fye, decimal=4) + + d_dy = Diff(1, y, acc=6) + fy = d_dy(f) + assert_array_almost_equal(fy, fye, decimal=4) diff --git a/tests/test_operators.py b/tests/test_operators.py index 080b91d..c18b80c 100644 --- a/tests/test_operators.py +++ b/tests/test_operators.py @@ -1,7 +1,7 @@ import numpy as np from numpy.testing import assert_array_almost_equal -from findiff.current import Add, FieldOperator +from findiff.operators import Add, FieldOperator from findiff import Diff, Identity @@ -100,7 +100,7 @@ def test_simple_derivatives(): dx = x[1] - x[0] f = x**3 - d_dx = Diff(dx) + d_dx = Diff(0, dx) actual = d_dx(f) np.testing.assert_array_almost_equal(actual, 3 * x**2, decimal=3) @@ -108,37 +108,66 @@ def test_simple_derivatives(): actual = d_dx(f, acc=4) np.testing.assert_array_almost_equal(actual, 3 * x**2) - d_dx = Diff(dx, 1, 0) + d_dx = Diff(0, dx) actual = d_dx(f, acc=4) np.testing.assert_array_almost_equal(actual, 3 * x**2) - d2_dx2 = Diff(dx, 2, 0) + d2_dx2 = Diff(0, dx) ** 2 actual = d2_dx2(f, acc=4) np.testing.assert_array_almost_equal(actual, 6 * x) +def test_pow_diff(): + + d = Diff(0) + assert d.order == 1 + + d2 = d**2 + + assert d2.order == 2 + assert d.order == 1 + + def test_chained_derivatives(): x = np.linspace(0, 1, 100) dx = x[1] - x[0] f = x**3 - d_dx = Diff(dx) + d_dx = Diff(0, dx) d2_dx2 = d_dx * d_dx actual = d2_dx2(f, acc=4) np.testing.assert_array_almost_equal(actual, 6 * x) +def test_set_grid_lazily(): + + x = np.linspace(0, 1, 100) + dx = x[1] - x[0] + f = x**3 + + # no grid defined here: + d_dx = Diff(0) + + diff_op = d_dx * d_dx + d_dx + + # now set grid for complete differential operator + diff_op.set_grid({0: dx}) + + actual = diff_op(f, acc=4) + np.testing.assert_array_almost_equal(actual, 6 * x + 3 * x**2) + + def test_harmonic(): x = np.linspace(0, 1, 100) dx = x[1] - x[0] f = x**3 - T = -0.5 * Diff(dx, 2, 0) + T = -0.5 * Diff(0, dx) ** 2 V = 0.5 * x**2 H = T + V @@ -156,9 +185,9 @@ def test_partial_diff(): u = np.sin(x) ux_ex = np.cos(x) dx = x[1] - x[0] - fd = Diff(dx, 1, 0) + fd = Diff(0, dx, acc=4) - ux = fd(u, spacing=dx, acc=4) + ux = fd(u) assert_array_almost_equal(ux, ux_ex, decimal=5) @@ -170,8 +199,68 @@ def test_partial_diff(): u = np.sin(X) * np.sin(Y) uxy_ex = np.cos(X) * np.cos(Y) - fd = Diff(dx, 1, 0) * Diff(dy, 1, 1) + # mixed partial derivative d2_dxdy + fd = Diff(0, dx) * Diff(1, dy) uxy = fd(u, acc=4) assert_array_almost_equal(uxy, uxy_ex, decimal=5) + + +def test_matrix_1d(): + + x = np.linspace(0, 6, 7) + d2_dx2 = Diff(0, x[1] - x[0]) ** 2 + u = x**2 + + mat = d2_dx2.matrix(u.shape) + + np.testing.assert_array_almost_equal(2 * np.ones_like(x), mat.dot(u.reshape(-1))) + + +def test_matrix_2d(): + thr = np.get_printoptions()["threshold"] + lw = np.get_printoptions()["linewidth"] + np.set_printoptions(threshold=np.inf) + np.set_printoptions(linewidth=500) + x, y = [np.linspace(0, 4, 5)] * 2 + X, Y = np.meshgrid(x, y, indexing="ij") + laplace = Diff(0, x[1] - x[0]) ** 2 + Diff(0, y[1] - y[0]) ** 2 + # d = FinDiff(1, y[1]-y[0], 2) + u = X**2 + Y**2 + + mat = laplace.matrix(u.shape) + + np.testing.assert_array_almost_equal( + 4 * np.ones_like(X).reshape(-1), mat.dot(u.reshape(-1)) + ) + + np.set_printoptions(threshold=thr) + np.set_printoptions(linewidth=lw) + + +def test_matrix_2d_mixed(): + x, y = [np.linspace(0, 5, 6), np.linspace(0, 6, 7)] + X, Y = np.meshgrid(x, y, indexing="ij") + d2_dxdy = Diff(0, x[1] - x[0]) * Diff(1, y[1] - y[0]) + u = X**2 * Y**2 + + mat = d2_dxdy.matrix(u.shape) + expected = d2_dxdy(u).reshape(-1) + + actual = mat.dot(u.reshape(-1)) + np.testing.assert_array_almost_equal(expected, actual) + + +def test_matrix_1d_coeffs(): + shape = (11,) + x = np.linspace(0, 10, 11) + dx = x[1] - x[0] + + L = x * Diff(0, dx) ** 2 + + u = np.random.rand(*shape).reshape(-1) + + actual = L.matrix(shape).dot(u) + expected = L(u).reshape(-1) + np.testing.assert_array_almost_equal(expected, actual) diff --git a/tests/test_pde.py b/tests/test_pde.py new file mode 100644 index 0000000..34c216d --- /dev/null +++ b/tests/test_pde.py @@ -0,0 +1,277 @@ +import numpy as np + +from findiff import Diff, BoundaryConditions, PDE, Identity +from tests.test_diff import make_grid + + +def test_1d_dirichlet_hom(): + + shape = (11,) + + x = np.linspace(0, 1, 11) + dx = x[1] - x[0] + L = Diff(0, dx) ** 2 + + bc = BoundaryConditions(shape) + + bc[0] = 1 + bc[-1] = 2 + + pde = PDE(L, np.zeros_like(x), bc) + u = pde.solve() + expected = x + 1 + np.testing.assert_array_almost_equal(expected, u) + + +def test_1d_dirichlet_inhom(): + + nx = 21 + shape = (nx,) + + x = np.linspace(0, 1, nx) + dx = x[1] - x[0] + L = Diff(0, dx) ** 2 + + bc = BoundaryConditions(shape) + + bc[0] = 1 + bc[-1] = 2 + + pde = PDE(L, 6 * x, bc) + + u = pde.solve() + expected = x**3 + 1 + np.testing.assert_array_almost_equal(expected, u) + + +def test_1d_neumann_hom(): + nx = 11 + shape = (nx,) + + x = np.linspace(0, 1, nx) + dx = x[1] - x[0] + L = Diff(0, dx) ** 2 + + bc = BoundaryConditions(shape) + + bc[0] = 1 + bc[-1] = Diff(0, dx), 2 + + pde = PDE(L, np.zeros_like(x), bc) + u = pde.solve() + expected = 2 * x + 1 + np.testing.assert_array_almost_equal(expected, u) + + +def test_2d_dirichlet_hom(): + shape = (11, 11) + + x, y = np.linspace(0, 1, shape[0]), np.linspace(0, 1, shape[1]) + dx, dy = x[1] - x[0], y[1] - y[0] + X, Y = np.meshgrid(x, y, indexing="ij") + L = Diff(0, dx) ** 2 + Diff(1, dy) ** 2 + + expected = X + 1 + + bc = BoundaryConditions(shape) + + bc[0, :] = 1 + bc[-1, :] = 2 + bc[:, 0] = X + 1 + bc[:, -1] = X + 1 + + pde = PDE(L, np.zeros_like(X), bc) + u = pde.solve() + + np.testing.assert_array_almost_equal(expected, u) + + +def test_2d_dirichlet_inhom(): + shape = (11, 11) + + x, y = np.linspace(0, 1, shape[0]), np.linspace(0, 1, shape[1]) + dx, dy = x[1] - x[0], y[1] - y[0] + X, Y = np.meshgrid(x, y, indexing="ij") + L = Diff(0, dx) ** 2 + Diff(1, dy) ** 2 + + expected = X**3 + Y**3 + 1 + f = 6 * X + 6 * Y + + bc = BoundaryConditions(shape) + + bc[0, :] = expected + bc[-1, :] = expected + bc[:, 0] = expected + bc[:, -1] = expected + + pde = PDE(L, f, bc) + u = pde.solve() + np.testing.assert_array_almost_equal(expected, u) + + +def test_2d_neumann_hom(): + shape = (31, 31) + + x, y = np.linspace(0, 1, shape[0]), np.linspace(0, 1, shape[1]) + dx, dy = x[1] - x[0], y[1] - y[0] + X, Y = np.meshgrid(x, y, indexing="ij") + L = Diff(0, dx) ** 2 + Diff(1, dy) ** 2 + + expected = X**2 + Y**2 + 1 + f = 4 * np.ones_like(X) + + bc = BoundaryConditions(shape) + + d_dy = Diff(1, dy) + + bc[0, :] = expected + bc[-1, :] = expected + bc[:, 0] = d_dy, 2 * Y + bc[:, -1] = d_dy, 2 * Y + + pde = PDE(L, f, bc) + u = pde.solve() + np.testing.assert_array_almost_equal(expected, u) + + +def test_1d_oscillator_free_dirichlet(): + n = 300 + shape = (n,) + t = np.linspace(0, 5, n) + dt = t[1] - t[0] + L = Diff(0, dt) ** 2 + Identity() + + bc = BoundaryConditions(shape) + + bc[0] = 1 + bc[-1] = 2 + + eq = PDE(L, np.zeros_like(t), bc) + u = eq.solve() + expected = np.cos(t) - (np.cos(5) - 2) * np.sin(t) / np.sin(5) + np.testing.assert_array_almost_equal(expected, u, decimal=4) + + +def test_1d_damped_osc_driv_dirichlet(): + n = 100 + shape = (n,) + t = np.linspace(0, 1, n) + dt = t[1] - t[0] + L = Diff(0, dt) ** 2 - Diff(0, dt) + Identity() + f = -3 * np.exp(-t) * np.cos(t) + 2 * np.exp(-t) * np.sin(t) + + expected = np.exp(-t) * np.sin(t) + + bc = BoundaryConditions(shape) + + bc[0] = expected[0] + bc[-1] = expected[-1] + + eq = PDE(L, f, bc) + u = eq.solve() + + np.testing.assert_array_almost_equal(expected, u, decimal=4) + + +def test_1d_oscillator_driv_neumann(): + n = 200 + shape = (n,) + t = np.linspace(0, 1, n) + dt = t[1] - t[0] + L = Diff(0, dt) ** 2 - Diff(0, dt) + Identity() + f = -3 * np.exp(-t) * np.cos(t) + 2 * np.exp(-t) * np.sin(t) + + expected = np.exp(-t) * np.sin(t) + + bc = BoundaryConditions(shape) + + bc[0] = Diff(0, dt), 1 + bc[-1] = expected[-1] + + eq = PDE(L, f, bc) + u = eq.solve() + + np.testing.assert_array_almost_equal(expected, u, decimal=4) + + +def test_1d_with_coeffs(): + n = 200 + shape = (n,) + t = np.linspace(0, 1, n) + dt = t[1] - t[0] + L = t * Diff(0, dt) ** 2 + f = 6 * t**2 + + bc = BoundaryConditions(shape) + + bc[0] = 0 + bc[-1] = 1 + + eq = PDE(L, f, bc) + u = eq.solve() + expected = t**3 + + np.testing.assert_array_almost_equal(expected, u, decimal=4) + + +def test_mixed_equation__with_coeffs_2d(): + shape = (41, 51) + + x, y = np.linspace(0, 1, shape[0]), np.linspace(0, 1, shape[1]) + dx, dy = x[1] - x[0], y[1] - y[0] + X, Y = np.meshgrid(x, y, indexing="ij") + L = Diff(0, dx) ** 2 + X * Y * Diff(0, dx) * Diff(1, dy) + Diff(1, dy) ** 2 + + expected = X**3 + Y**3 + 1 + f = 6 * (X + Y) + + bc = BoundaryConditions(shape) + + bc[0, :] = expected + bc[-1, :] = expected + bc[:, 0] = expected + bc[:, -1] = expected + + pde = PDE(L, f, bc) + u = pde.solve() + np.testing.assert_array_almost_equal(expected, u, decimal=4) + + +def test_2d_inhom_const_coefs_dirichlet_all(): + shape = (41, 50) + (x, y), (dx, dy), (X, Y) = make_grid(shape, edges=[(-1, 1), (-1, 1)]) + + expected = X**3 + Y**3 + X * Y + 1 + + L = 3 * Diff(0, dx) ** 2 + 2 * Diff(0, dx) * Diff(1, dy) + Diff(1, dy) ** 2 + f = 2 + 18 * X + 6 * Y + + bc = BoundaryConditions(shape) + bc[0, :] = expected + bc[-1, :] = expected + bc[:, 0] = expected + bc[:, -1] = expected + + pde = PDE(L, f, bc) + actual = pde.solve() + np.testing.assert_array_almost_equal(expected, actual, decimal=4) + + +def test_2d_inhom_var_coefs_dirichlet_all(): + shape = (41, 50) + (x, y), (dx, dy), (X, Y) = make_grid(shape, edges=[(-1, 1), (-1, 1)]) + + expected = X**3 + Y**3 + X * Y + 1 + + L = 3 * X * Diff(0, dx) ** 2 + 2 * Y * Diff(0, dx) * Diff(1, dy) + Diff(1, dy) ** 2 + f = 18 * X**2 + 8 * Y + + bc = BoundaryConditions(shape) + bc[0, :] = expected + bc[-1, :] = expected + bc[:, 0] = expected + bc[:, -1] = expected + + pde = PDE(L, f, bc) + actual = pde.solve() + np.testing.assert_array_almost_equal(expected, actual, decimal=4) From 42dd1962c4847238fd63280f7d96c589a78e8316 Mon Sep 17 00:00:00 2001 From: maroba Date: Wed, 11 Dec 2024 12:11:11 +0100 Subject: [PATCH 4/6] Prepare version 0.11.0 --- .github/workflows/check.yml | 1 - CHANGES.md | 4 + findiff/__init__.py | 2 +- findiff/legacy/diff.py | 265 +----------------- findiff/legacy/operators.py | 23 +- tests/legacy/__init__.py | 0 tests/{legacy => }/test_bugs.py | 0 tests/{legacy => }/test_coefs.py | 0 tests/{legacy => }/test_findiff.py | 0 tests/{legacy/test_pde.py => test_pde_old.py} | 5 +- tests/{legacy => }/test_scaling.py | 0 tests/{legacy => }/test_stencils.py | 0 tests/{legacy => }/test_symbolic.py | 0 tests/{legacy => }/test_utils.py | 0 tests/{legacy => }/test_vector.py | 0 15 files changed, 17 insertions(+), 283 deletions(-) delete mode 100644 tests/legacy/__init__.py rename tests/{legacy => }/test_bugs.py (100%) rename tests/{legacy => }/test_coefs.py (100%) rename tests/{legacy => }/test_findiff.py (100%) rename tests/{legacy/test_pde.py => test_pde_old.py} (98%) rename tests/{legacy => }/test_scaling.py (100%) rename tests/{legacy => }/test_stencils.py (100%) rename tests/{legacy => }/test_symbolic.py (100%) rename tests/{legacy => }/test_utils.py (100%) rename tests/{legacy => }/test_vector.py (100%) diff --git a/.github/workflows/check.yml b/.github/workflows/check.yml index df68973..73b37a5 100644 --- a/.github/workflows/check.yml +++ b/.github/workflows/check.yml @@ -2,7 +2,6 @@ name: Checks on: push: - branches: [ master ] jobs: build: diff --git a/CHANGES.md b/CHANGES.md index ae00975..cf8d1d8 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,5 +1,9 @@ # Change Log +## Version 0.11.x + +- More comfortable API. (Old API still available) + ## Version 0.10.x - Create symbolic representations of finite difference schemes diff --git a/findiff/__init__.py b/findiff/__init__.py index 6eedddc..33f1d32 100644 --- a/findiff/__init__.py +++ b/findiff/__init__.py @@ -18,7 +18,7 @@ - Version 1.0: Completely remodeled API (backward compatibility is maintained, though) """ -__version__ = "1.0.0.dev" +__version__ = "0.11.0" from .legacy import * diff --git a/findiff/legacy/diff.py b/findiff/legacy/diff.py index 281a9a6..68c1341 100644 --- a/findiff/legacy/diff.py +++ b/findiff/legacy/diff.py @@ -1,226 +1,15 @@ import itertools import numbers -import operator import scipy.sparse as sparse from ..coefs import coefficients, coefficients_non_uni -from ..stencils import StencilSet from ..utils import * DEFAULT_ACC = 2 -class Operator(object): - """Base class for all operator classes""" - - pass - - -class BinaryOperator(Operator): - """Base class for all binary operators (like addition) that allow to combine or chain differential operators""" - - def __init__(self, left, right): - self.left = left - self.right = right - - def apply(self, rhs, *args, **kwargs): - raise NotImplementedError - - def __call__(self, rhs, *args, **kwargs): - return self.apply(rhs, *args, **kwargs) - - def stencil(self, shape): - return StencilSet(self, shape) - - -class Plus(BinaryOperator): - """A class to add two differential operators""" - - def __init__(self, left, right): - super().__init__(left, right) - - def __add__(self, other): - return Plus(self, other) - - def __radd__(self, other): - return Plus(self, other) - - def __sub__(self, other): - return Minus(self, other) - - def __rsub__(self, other): - return Minus(self, other) - - def __mul__(self, other): - return Mul(self, other) - - def __rmul__(self, other): - return Mul(self, other) - - def apply(self, rhs, *args, **kwargs): - - if isinstance(self.right, LinearMap) or isinstance(self.right, BinaryOperator): - right = self.right.apply(rhs, *args, **kwargs) - else: - right = self.right - - if isinstance(self.left, LinearMap) or isinstance(self.left, BinaryOperator): - left = self.left.apply(rhs, *args, **kwargs) - else: - left = self.left * rhs - - return left + right - - def matrix(self, shape, *args, **kwargs): - left, right = self.left, self.right - if isinstance(self.left, Operator): - left = self.left.matrix(shape, *args, **kwargs) - if isinstance(self.right, Operator): - right = self.right.matrix(shape, *args, **kwargs) - return left + right - - -class Minus(BinaryOperator): - """A class to subtract two differential operators from each other""" - - def __init__(self, left, right): - super().__init__(left, right) - - def __add__(self, other): - return Plus(self, other) - - def __radd__(self, other): - return Plus(self, other) - - def __sub__(self, other): - return Minus(self, other) - - def __rsub__(self, other): - return Minus(self, other) - - def __mul__(self, other): - return Mul(self, other) - - def __rmul__(self, other): - return Mul(self, other) - - def apply(self, rhs, *args, **kwargs): - - if isinstance(self.right, LinearMap) or isinstance(self.right, BinaryOperator): - right = self.right.apply(rhs, *args, **kwargs) - else: - right = self.right - - if isinstance(self.left, LinearMap) or isinstance(self.left, BinaryOperator): - left = self.left.apply(rhs, *args, **kwargs) - else: - left = self.left * rhs - - return left - right - - def matrix(self, shape, *args, **kwargs): - left, right = self.left, self.right - if isinstance(self.left, Operator): - left = self.left.matrix(shape, *args, **kwargs) - if isinstance(self.right, Operator): - right = self.right.matrix(shape, *args, **kwargs) - return left - right - - -class Mul(BinaryOperator): - """A class to multiply (chain) two differential operators""" - - def __init__(self, left, right): - super().__init__(left, right) - self.oper = operator.mul - - def __add__(self, other): - return Plus(self, other) - - def __radd__(self, other): - return Plus(self, other) - - def __sub__(self, other): - return Minus(self, other) - - def __rsub__(self, other): - return Minus(self, other) - - def __mul__(self, other): - return Mul(self, other) - - def __rmul__(self, other): - return Mul(self, other) - - def apply(self, rhs, *args, **kwargs): - - if isinstance(self.right, LinearMap) or isinstance(self.right, BinaryOperator): - result = self.right.apply(rhs, *args, **kwargs) - else: - result = self.right * rhs - - if isinstance(self.left, LinearMap) or isinstance(self.left, BinaryOperator): - result = self.left.apply(result, *args, **kwargs) - else: - result = self.left * result - - return result - - def matrix(self, shape, *args, **kwargs): - """Matrix representation of given operator product on an equidistant grid of given shape. - - :param shape: tuple with the shape of the grid - :return: scipy sparse matrix representing the operator product - """ - - if isinstance(self.left, np.ndarray): - left = sparse.diags(self.left.reshape(-1), 0) - elif isinstance(self.left, LinearMap) or isinstance(self.left, BinaryOperator): - left = self.left.matrix(shape, *args, **kwargs) - else: - left = self.left * sparse.diags(np.ones(shape).reshape(-1), 0) - - if isinstance(self.right, np.ndarray): - right = sparse.diags(self.right.reshape(-1), 0) - elif isinstance(self.right, LinearMap) or isinstance( - self.right, BinaryOperator - ): - right = self.right.matrix(shape, *args, **kwargs) - else: - right = self.right * sparse.diags(np.ones(shape).reshape(-1), 0) - - return left.dot(right) - - -class LinearMap(Operator): - - def __init__(self, value): - self.value = value - - def __add__(self, other): - return Plus(self, other) - - def __radd__(self, other): - return Plus(self, other) - - def __sub__(self, other): - return Minus(self, other) - - def __rsub__(self, other): - return Minus(self, other) - - def __mul__(self, other): - return Mul(self, other) - - def __rmul__(self, other): - return Mul(self, other) - - def __call__(self, rhs, *args, **kwargs): - return self.apply(rhs, *args, **kwargs) - - -class Diff(LinearMap): +class Diff: """Representation of a single partial derivative based on finite differences. This class is usually not used directly by the user, but is wrapped in @@ -249,6 +38,9 @@ def __init__(self, axis, order, **kwargs): if "acc" in kwargs: self.acc = kwargs["acc"] + def __call__(self, rhs, *args, **kwargs): + return self.apply(rhs, *args, **kwargs) + def apply(self, u, *args, **kwargs): """Applies the partial derivative to a numpy array.""" @@ -512,52 +304,3 @@ def _shift_slice(self, sl, off, max_index): raise IndexError("Shift slice out of bounds") return slice(sl.start + off, sl.stop + off, sl.step) - - -class Id(LinearMap): - """The identity operator. When applied to an array, returns the same array (not a copy)""" - - def __init__(self): - self.value = 1 - - def apply(self, rhs, *args, **kwargs): - return rhs - - def matrix(self, shape): - """Matrix representation of the identity operator, i.e. identity matrix of given shape. - - :param shape: Shape of the arrays to which Id shall be applied - :type shape: tuple of ints - :return: Sparse identity matrix. - :rtype: scipy.sparse.csr_matrix - """ - - siz = np.prod(shape) - mat = sparse.lil_matrix((siz, siz)) - diag = list(range(siz)) - mat[diag, diag] = 1 - return sparse.csr_matrix(mat) - - -class Coef(object): - """ - Encapsulates a constant (number) or variable (N-dimensional coordinate array) value to multiply with a linear operator - - :param value: a number or an numpy.ndarray with meshed coordinates - - **Example**: - - The following example defines the differential operator - - .. math:: 2x\\frac{\\partial^3}{\\partial x^2\\partial z} - - >>> X, Y, Z, U = numpy.meshgrid(x, y, z, u, indexing="ij") - >>> diff_op = Coef(2*X) * FinDiff((0, dx, 2), (2, dz, 1)) - - """ - - def __init__(self, value): - self.value = value - - def __mul__(self, other): - return Mul(self.value, other) diff --git a/findiff/legacy/operators.py b/findiff/legacy/operators.py index 00e87b9..529c08f 100644 --- a/findiff/legacy/operators.py +++ b/findiff/legacy/operators.py @@ -1,9 +1,10 @@ -from ..utils import * -from .diff import Coef, Id, Diff, Plus, Minus, Mul, DEFAULT_ACC, LinearMap -from ..stencils import StencilSet +from findiff.legacy.diff import Diff +from findiff.stencils import StencilSet +DEFAULT_ACC = 2 -class FinDiff(LinearMap): + +class FinDiff: r"""A representation of a general linear differential operator expressed in finite differences. FinDiff objects can be added with other FinDiff objects. They can be multiplied by @@ -112,15 +113,6 @@ def matrix(self, shape, h=None, acc=None): def set_accuracy(self, acc): self.pds.set_accuracy(acc) - def __add__(self, other): - return Plus(self, other) - - def __sub__(self, other): - return Minus(self, other) - - def __mul__(self, other): - return Mul(self, other) - def _eval_args(self, args, kwargs): spac = {} @@ -168,8 +160,3 @@ def _eval_args(self, args, kwargs): break return pds - - -# Alias for backward compatibility -Coefficient = Coef -Identity = Id diff --git a/tests/legacy/__init__.py b/tests/legacy/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/tests/legacy/test_bugs.py b/tests/test_bugs.py similarity index 100% rename from tests/legacy/test_bugs.py rename to tests/test_bugs.py diff --git a/tests/legacy/test_coefs.py b/tests/test_coefs.py similarity index 100% rename from tests/legacy/test_coefs.py rename to tests/test_coefs.py diff --git a/tests/legacy/test_findiff.py b/tests/test_findiff.py similarity index 100% rename from tests/legacy/test_findiff.py rename to tests/test_findiff.py diff --git a/tests/legacy/test_pde.py b/tests/test_pde_old.py similarity index 98% rename from tests/legacy/test_pde.py rename to tests/test_pde_old.py index 99da46e..5aa9316 100644 --- a/tests/legacy/test_pde.py +++ b/tests/test_pde_old.py @@ -1,10 +1,11 @@ import sys +import numpy as np + sys.path.insert(1, "..") import unittest -from findiff.legacy.operators import FinDiff, Identity, Coef -from findiff.pde import * +from findiff import * # import matplotlib.pyplot as plt # from mpl_toolkits import mplot3d diff --git a/tests/legacy/test_scaling.py b/tests/test_scaling.py similarity index 100% rename from tests/legacy/test_scaling.py rename to tests/test_scaling.py diff --git a/tests/legacy/test_stencils.py b/tests/test_stencils.py similarity index 100% rename from tests/legacy/test_stencils.py rename to tests/test_stencils.py diff --git a/tests/legacy/test_symbolic.py b/tests/test_symbolic.py similarity index 100% rename from tests/legacy/test_symbolic.py rename to tests/test_symbolic.py diff --git a/tests/legacy/test_utils.py b/tests/test_utils.py similarity index 100% rename from tests/legacy/test_utils.py rename to tests/test_utils.py diff --git a/tests/legacy/test_vector.py b/tests/test_vector.py similarity index 100% rename from tests/legacy/test_vector.py rename to tests/test_vector.py From 40b976dcf10e2d47e750a2ea8a52be7b844f4bfc Mon Sep 17 00:00:00 2001 From: maroba Date: Wed, 11 Dec 2024 12:27:57 +0100 Subject: [PATCH 5/6] use pyproject.toml --- .gitignore | 1 + README.md | 2 +- pyproject.toml | 23 ++++++++++++++++------- requirements.txt | 6 ------ setup.py | 25 ------------------------- 5 files changed, 18 insertions(+), 39 deletions(-) delete mode 100644 requirements.txt delete mode 100644 setup.py diff --git a/.gitignore b/.gitignore index 385af19..e69bd33 100644 --- a/.gitignore +++ b/.gitignore @@ -106,3 +106,4 @@ ENV/ # Visual Studio Code .vscode +venv2 diff --git a/README.md b/README.md index 7d165f3..174f1d7 100644 --- a/README.md +++ b/README.md @@ -57,7 +57,7 @@ df_dx = d_dx(f) ``` Similarly, you can define partial derivatives along other axes, for example, if $z$ is the 2-axis, we can write -$\frac{partial}{\partial z}$ as: +$\frac{\partial}{\partial z}$ as: ```python Diff(2, dz) diff --git a/pyproject.toml b/pyproject.toml index 1b66482..14add72 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,17 +1,25 @@ -# pyproject.toml +[build-system] +requires = ["setuptools>=64", "setuptools_scm[toml]>=6.2"] +build-backend = "setuptools.build_meta" + [project] name = "findiff" authors = [ - { name = "Matthias Baer", email = "matthias.r.baer@googlemail.com" }, + {name = "Matthias Baer"}, +] +maintainers = [ + {name = "Matthias Baer", email="matthias.r.baer@googlemail.com"}, ] description = "A Python package for finite difference derivatives in any number of dimensions." -dynamic = ["version"] readme = "README.md" +license = {text = "MIT"} dependencies = [ "numpy", "scipy", "sympy" ] +dynamic = ["version"] + requires-python = ">=3.8" keywords = ["finite-differences", "numerical-derivatives", "scientific-computing"] classifiers = [ @@ -27,13 +35,14 @@ classifiers = [ "Operating System :: OS Independent", ] -[build-system] -requires = ["setuptools", "wheel"] -build-backend = "setuptools.build_meta" - [project.urls] Homepage = "https://github.com/maroba/findiff" +source = "https://github.com/maroba/findiff" Issues = "https://github.com/maroba/findiff/issues" +tracker = "https://github.com/maroba/findiff/issues" + +[tool.setuptools.packages.find] +include = ["findiff"] [tool.setuptools.dynamic] version = { attr = "findiff.__version__" } \ No newline at end of file diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index 67f91ef..0000000 --- a/requirements.txt +++ /dev/null @@ -1,6 +0,0 @@ -numpy>=1.26.4 -scipy>=1.12.0 -Sphinx>=7.2.6 -sympy - - diff --git a/setup.py b/setup.py deleted file mode 100644 index affa99b..0000000 --- a/setup.py +++ /dev/null @@ -1,25 +0,0 @@ -from setuptools import setup, find_packages - - -setup( - name="findiff", - packages=find_packages(where=".", include=["findiff*"]), - long_description="""A Python package for finite difference derivatives in any number of dimensions. - - Features: - - * Differentiate arrays of any number of dimensions along any axis - * Partial derivatives of any desired order - * Accuracy order can be specified - * Accurate treatment of grid boundary - * Includes standard operators from vector calculus like gradient, divergence and curl - * Can handle uniform and non-uniform grids - * Can handle arbitrary linear combinations of derivatives with constant and variable coefficients - * Fully vectorized for speed - * Calculate raw finite difference coefficients for any order and accuracy for uniform and non-uniform grids - * _New in version 0.7:_ Generate matrix representations of arbitrary linear differential operators - * _New in version 0.8:_ Solve partial differential equations with Dirichlet or Neumann boundary conditions - * _New in version 0.9:_ Generate differential operators for generic stencils - * _New in version 0.10:_ Create symbolic representations of finite difference schemes - """, -) From 71e40a6b894116a3cbc2d3382feefaffdbe99716 Mon Sep 17 00:00:00 2001 From: maroba Date: Wed, 11 Dec 2024 12:29:23 +0100 Subject: [PATCH 6/6] Update README --- README.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/README.md b/README.md index 174f1d7..7ead90f 100644 --- a/README.md +++ b/README.md @@ -44,13 +44,12 @@ from findiff import Diff # define the grid: x = np.linspace(0, 1, 100) -dx = x[1] - x[0] # the grid spacing # the array to differentiate: f = np.sin(x) # as an example # Define the derivative: -d_dx = Diff(0, dx) +d_dx = Diff(0, x[1] - x[0]) # Apply it: df_dx = d_dx(f)