Truncation error eradicated :). Check helmholtz next

Author: hirish99

sumpy/recurrence.py

@@ -473,19 +473,23 @@ def get_off_axis_expression(pde, taylor_order=4) -> [sp.Expr, int]:
     r"""
     A function that takes in as input a pde, and outputs
     the Taylor expression that gives the n th derivative
-    as a truncated taylor_order th order Taylor series with respect to x_1 and
+    as a truncated taylor_order th order Taylor series with respect to x_0 and
     s(i) where s(i) comes from the off-axis recurrence. See
     get_reindexed_and_center_origin_off_axis_recurrence.
 
     Also outputs the number of coefficients it needs from nth order.
     So if it outputs 3 as the second return value, then it needs
     s(deriv_order), s(deriv_order-1), ..., s(deriv_order-3).
+
+    YOU CANNOT SUB N < START_ORDER INTO THE OUTPUTTED EXPRESSION.
+    I CANNOT REARRANGE THE EXPRESSION IN THIS CASE TO HAVE INDICES
+    LOWER THAN THE SUBSITUTED N VALUE.
     """
     s = sp.Function("s")
     n = sp.symbols("n")
     deriv_order = n
 
-    t_recurrence = get_reindexed_and_center_origin_off_axis_recurrence(pde)[2]
+    start_order, _, t_recurrence = get_reindexed_and_center_origin_off_axis_recurrence(pde)
     var = _make_sympy_vec("x", 2)
     exp = 0
     for i in range(taylor_order+1):
@@ -507,4 +511,4 @@ def get_off_axis_expression(pde, taylor_order=4) -> [sp.Expr, int]:
 
     idx_l, _ = _extract_idx_terms_from_recurrence(exp)
 
-    return exp*(-1)**n, -min(idx_l)
+    return exp*(-1)**n, -min(idx_l), start_order

sumpy/recurrence_qbx.py

@@ -153,7 +153,7 @@ def generate_lamb_expr(i, n_initial):
 
     ### NEW CODE - COMPUTE OFF AXIS INTERACTIONS
     start_order, t_recur_order, t_recur = get_reindexed_and_center_origin_off_axis_recurrence(pde)
-    t_exp, t_exp_order = get_off_axis_expression(pde, 8)
+    t_exp, t_exp_order, _ = get_off_axis_expression(pde, 8)
     storage_taylor_order = max(t_recur_order, t_exp_order+1)
 
     storage_taylor = [np.zeros((n_p, n_p))] * storage_taylor_order
@@ -176,15 +176,15 @@ def gen_lamb_expr_t_recur(i, start_order):
         return sp.lambdify(arg_list, lamb_expr_symb)
 
 
-    def gen_lamb_expr_t_exp(i, t_exp_order):
+    def gen_lamb_expr_t_exp(i, t_exp_order, start_order):
         arg_list = []
         for j in range(t_exp_order, -1, -1):
             # pylint: disable-next=not-callable
             arg_list.append(s(i-j))
         for j in range(ndim):
             arg_list.append(var[j])
 
-        if i < t_exp_order:
+        if i < start_order:
             lamb_expr_symb_deriv = sp.diff(g_x_y, var_t[0], i)
             for j in range(ndim):
                 lamb_expr_symb_deriv = lamb_expr_symb_deriv.subs(var_t[j], 0)
@@ -203,7 +203,7 @@ def gen_lamb_expr_t_exp(i, t_exp_order):
         storage_taylor.pop(0)
         storage_taylor.append(lamb_expr_t_recur(*a1) + np.zeros((n_p, n_p)))
 
-        lamb_expr_t_exp = gen_lamb_expr_t_exp(i, t_exp_order)
+        lamb_expr_t_exp = gen_lamb_expr_t_exp(i, t_exp_order, start_order)
         a2 = [*storage_taylor[-(t_exp_order+1):], *coord]
 
         interactions_off_axis += lamb_expr_t_exp(*a2) * radius**i/math.factorial(i)
@@ -231,7 +231,6 @@ def generate_true(i):
         a4 = [*coord]
         s_new_true = lamb_expr_true(*a4)
         interactions_true += s_new_true * radius**i/math.factorial(i)
-
     ###############
 
     #slope of line y = mx
@@ -251,14 +250,13 @@ def generate_true(i):
     print("Y:", coord[1][mask_on_axis].reshape(-1)[np.argmax(relerr_on)])
 
     print("-------------------------")
-  
+
     if np.any(mask_off_axis):
         relerr_off = np.abs(interactions_off_axis[mask_off_axis]-interactions_true[mask_off_axis])/np.abs(interactions_off_axis[mask_off_axis])
         print("MAX OFF AXIS ERROR(", percent_off, "):", np.max(relerr_off))
         print(np.mean(relerr_off))
         print("X:", coord[0][mask_off_axis].reshape(-1)[np.argmax(relerr_off)])
-        print("Y:", coord[1][mask_off_axis].reshape(-1)[np.argmax(relerr_off)])  
-
+        print("Y:", coord[1][mask_off_axis].reshape(-1)[np.argmax(relerr_off)])
 
     interactions_total = np.zeros(coord[0].shape)
     interactions_total[mask_on_axis] = interactions_on_axis[mask_on_axis]

test/test_recurrence.py

@@ -257,7 +257,7 @@ def test_helmholtz_2d_off_axis(deriv_order, exp_order):
 
     # CHECK ACCURACY OF EXPRESSION FOR deriv_order
 
-    exp, exp_range = get_off_axis_expression(helmholtz2d, exp_order)
+    exp, exp_range, _ = get_off_axis_expression(helmholtz2d, exp_order)
     approx_deriv = exp.subs(n, deriv_order)
     for i in range(-exp_range+deriv_order, deriv_order+1):
         approx_deriv = approx_deriv.subs(s(i), ic[i])
@@ -317,7 +317,9 @@ def test_laplace_2d_off_axis(deriv_order, exp_order):
 
     # CHECK ACCURACY OF EXPRESSION FOR deriv_order
 
-    exp, exp_range = get_off_axis_expression(laplace2d, exp_order)
+    exp, exp_range, start_order = get_off_axis_expression(laplace2d, exp_order)
+    
+    assert deriv_order >= start_order
     approx_deriv = exp.subs(n, deriv_order)
     for i in range(-exp_range+deriv_order, deriv_order+1):
         approx_deriv = approx_deriv.subs(s(i), ic[i])

test/test_recurrence_qbx.py

@@ -294,15 +294,15 @@ def _construct_laplace_axis_2d(orders, resolutions):
     return err
 
 import matplotlib.pyplot as plt
-orders = [1]
+orders = [10, 16]
 #resolutions = range(200, 800, 200)
-resolutions = [400]
+resolutions = [400, 800, 1200]
 err_mat = _construct_laplace_axis_2d(orders, resolutions)
 
 for i in range(len(orders)):
-    plt.plot(resolutions, err_mat[i], label="order QBX="+str(orders[i]))
+    plt.scatter(resolutions, err_mat[i], label="order QBX="+str(orders[i]))
 plt.xlabel("Number of Nodes")
-plt.ylabel("Error")
-plt.title("2D Ellipse LP Eval Error")
+plt.ylabel("Relative Error")
+plt.title("2D Ellipse LP Eval Error (m=10)")
 plt.legend()
 plt.show()