diff --git a/CHANGES.md b/CHANGES.md index 8894342..088411a 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -7,6 +7,7 @@ - Iterative solvers for PDEs (CG, GMRES, BiCGSTAB, LGMRES, MINRES) with preconditioner support - Eigenvalue problem support (eigs and eigsh methods) - Robin (mixed) boundary conditions for PDEs +- Fix order-of-accuracy loss for even-order derivatives on non-uniform grids (central stencil was undersized) ## Version 0.13.x diff --git a/findiff/coefs.py b/findiff/coefs.py index 4bfa4f1..eeddb98 100644 --- a/findiff/coefs.py +++ b/findiff/coefs.py @@ -159,13 +159,16 @@ def coefficients_non_uni(deriv, acc, coords, idx): _validate_acc(acc) num_central = 2 * math.floor((deriv + 1) / 2) - 1 + acc - num_side = num_central // 2 if deriv % 2 == 0: num_coef = num_central + 1 else: num_coef = num_central + # Non-uniform central stencils get no symmetry bonus, so size them like the + # one-sided schemes (num_coef points) to reach the requested order `acc`. + num_side = num_coef // 2 + if idx < num_side: matrix = _build_matrix_non_uniform(0, num_coef - 1, coords, idx) @@ -219,8 +222,10 @@ def calc_coefs_non_uni_batched(deriv, acc, coords): _validate_acc(acc) num_central = 2 * math.floor((deriv + 1) / 2) - 1 + acc - num_side = num_central // 2 num_coef = num_central + 1 if deriv % 2 == 0 else num_central + # Non-uniform central stencils get no symmetry bonus, so size them like the + # one-sided schemes (num_coef points) to reach the requested order `acc`. + num_side = num_coef // 2 N = len(coords) diff --git a/tests/test_coefs.py b/tests/test_coefs.py index 346bb27..4a380e5 100644 --- a/tests/test_coefs.py +++ b/tests/test_coefs.py @@ -170,12 +170,18 @@ def test_non_uniform(analytic_inv): x = np.linspace(0, 10, 100) dx = x[1] - x[0] - c_uni = coefficients(deriv=2, acc=2, analytic_inv=analytic_inv) - 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"] + # On a uniform grid the non-uniform scheme must reproduce the equidistant + # coefficients for the stencil it actually uses. + coefs_uni = ( + calc_coefs(2, list(c_non_uni["offsets"]), analytic_inv=analytic_inv)[ + "coefficients" + ] + / dx**2 + ) + np.testing.assert_array_almost_equal(coefs_non_uni, coefs_uni) diff --git a/tests/test_nonuniform_accuracy.py b/tests/test_nonuniform_accuracy.py new file mode 100644 index 0000000..910ee33 --- /dev/null +++ b/tests/test_nonuniform_accuracy.py @@ -0,0 +1,86 @@ +"""Accuracy-order regression tests for non-uniform grids. + +On a non-uniform grid the central stencil earns no symmetry bonus, so the +interior scheme for an even-order derivative must carry as many points as the +one-sided schemes to reach the requested order. Before that fix the even-order +interior schemes were one point short and only reached order ``acc - 1`` (the +non-uniform Laplacian was only first-order in the interior). Odd orders were +already correct. These tests pin the order down by polynomial exactness and by +a graded-grid convergence study. +""" + +import numpy as np +import pytest + +from findiff import Diff + + +def _graded_grid(n, ratio=2.0, length=3.0): + """Grid with alternating spacings h, ratio*h, ... (fixed O(1) asymmetry). + + A smoothly mapped grid hides the order loss (its asymmetry vanishes as the + grid is refined); a graded/boundary-layer mesh does not. + """ + steps = np.empty(n - 1) + steps[0::2] = 1.0 + steps[1::2] = ratio + x = np.concatenate([[0.0], np.cumsum(steps)]) + return x / x[-1] * length + + +@pytest.mark.parametrize("acc", [2, 4]) +@pytest.mark.parametrize("deriv", [1, 2, 3, 4]) +def test_nonuniform_polynomial_exactness(deriv, acc): + """An order-acc scheme is exact for polynomials of degree deriv+acc-1 at + every grid point. This failed for even derivatives before the fix.""" + x = _graded_grid(40) + degree = deriv + acc - 1 + p = np.poly1d(np.random.default_rng(0).standard_normal(degree + 1)) + f = p(x) + op = Diff(0, x, acc=acc) ** deriv if deriv > 1 else Diff(0, x, acc=acc) + got = op(f) + exact = p.deriv(deriv)(x) + scale = max(1.0, np.max(np.abs(exact))) + assert np.max(np.abs(got - exact)) / scale < 1e-8 + + +@pytest.mark.parametrize("deriv,acc", [(1, 2), (1, 4), (2, 2), (2, 4)]) +def test_nonuniform_convergence_order(deriv, acc): + """Interior convergence order on a graded grid reaches the requested acc. + + Even derivatives dropped to order acc-1 before the central stencil was + widened; odd derivatives were already correct. + """ + + def f_and_deriv(x): + f = np.exp(np.sin(2 * x)) + if deriv == 1: + return f, 2 * np.cos(2 * x) * f + return f, (2 * np.cos(2 * x)) ** 2 * f - 4 * np.sin(2 * x) * f + + errs = [] + for n in (41, 81, 161): + x = _graded_grid(n) + f, exact = f_and_deriv(x) + op = Diff(0, x, acc=acc) ** deriv if deriv > 1 else Diff(0, x, acc=acc) + got = op(f) + b = deriv + acc + 1 # skip a few boundary points + errs.append(np.max(np.abs(got[b:-b] - exact[b:-b]))) + orders = np.log2(np.array(errs[:-1]) / np.array(errs[1:])) + assert orders.min() > acc - 0.4 + + +def test_nonuniform_even_derivative_interior_matches_boundary_order(): + """Interior even-derivative error is no worse than the boundary error. + + The undersized central stencil made the interior *less* accurate than the + one-sided boundary stencils, which is backwards. + """ + x = _graded_grid(60) + f = np.exp(np.sin(2 * x)) + exact = (2 * np.cos(2 * x)) ** 2 * f - 4 * np.sin(2 * x) * f + got = (Diff(0, x, acc=2) ** 2)(f) + err = np.abs(got - exact) + interior = err[5:-5].max() + boundary = max(err[:5].max(), err[-5:].max()) + assert interior <= boundary diff --git a/tests/test_nonuniform_optimized.py b/tests/test_nonuniform_optimized.py index 682d231..d6c50d9 100644 --- a/tests/test_nonuniform_optimized.py +++ b/tests/test_nonuniform_optimized.py @@ -101,15 +101,17 @@ def test_batched_result_structure(self): assert result["center"]["coefficients"].shape[0] == len(x) - 2 * num_bndry def test_batched_even_derivative(self): - """Even derivatives have forward/backward stencils wider than central.""" + """Even-derivative central stencils get no symmetry bonus on a + non-uniform grid, so they carry deriv+acc points (one more than the + one-sided schemes) to reach order acc.""" x = np.linspace(0, 5, 20) result = calc_coefs_non_uni_batched(deriv=2, acc=2, coords=x) central_size = result["center"]["coefficients"].shape[1] forward_size = result["forward"]["coefficients"].shape[1] - # For even derivatives, forward/backward have one more point - assert forward_size == central_size + 1 + assert central_size == forward_size + 1 + assert central_size >= 2 + 2 # deriv + acc points needed for order acc def test_batched_invalid_inputs(self): """Invalid deriv/acc raise ValueError.""" diff --git a/tests/test_vector.py b/tests/test_vector.py index 767b81b..7bec262 100644 --- a/tests/test_vector.py +++ b/tests/test_vector.py @@ -196,13 +196,24 @@ def test_1d_non_uniform(self): laplace = Laplacian(coords=[x], acc=4) assert_array_almost_equal(laplace(f), expected, decimal=4) - def test_matches_uniform_result(self): - """Non-uniform Laplacian on a uniform grid should match the uniform version.""" + def test_uniform_grid_gives_correct_laplacian(self): + """On a uniform grid the coords= Laplacian is correct and at least as + accurate as the h= path. Its central stencils are wider (the non-uniform + scheme cannot rely on the symmetry bonus), so the two paths agree only to + the uniform scheme's order, not bit-for-bit.""" axes, h, [X, Y] = init_mesh(2, (50, 50)) f = np.sin(X) * np.cos(Y) + expected = -2 * np.sin(X) * np.cos(Y) result_uniform = Laplacian(h=h, acc=2)(f) result_nonuniform = Laplacian(coords=axes, acc=2)(f) - assert_array_almost_equal(result_nonuniform, result_uniform) + + interior = (slice(2, -2), slice(2, -2)) + assert_array_almost_equal( + result_nonuniform[interior], expected[interior], decimal=3 + ) + err_uniform = np.max(np.abs(result_uniform[interior] - expected[interior])) + err_nonuniform = np.max(np.abs(result_nonuniform[interior] - expected[interior])) + assert err_nonuniform <= err_uniform def test_h_and_coords_raises(self): axes, h, _ = init_mesh(2, (10, 10))