Issue with mask_off_axis

Author: hirish99

sumpy/recurrence.py

@@ -477,8 +477,8 @@ def get_off_axis_expression(pde, taylor_order=4) -> [sp.Expr, int]:
     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
+    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).
     """
     s = sp.Function("s")
@@ -506,6 +506,5 @@ def get_off_axis_expression(pde, taylor_order=4) -> [sp.Expr, int]:
         max_idx = max(idx_l) 
 
     idx_l, _ = _extract_idx_terms_from_recurrence(exp)
-    exp_range = (min(idx_l), max(idx_l))
 
-    return exp*(-1)**n, min(idx_l)
+    return exp*(-1)**n, -min(idx_l)

sumpy/recurrence_qbx.py

@@ -156,7 +156,7 @@ def generate_lamb_expr(i, n_initial):
     t_exp, t_exp_order = get_off_axis_expression(pde)
     storage_taylor_order = max(t_recur_order, t_exp_order+1)
 
-    storage = [np.zeros((n_p, n_p))] * storage_taylor_order
+    storage_taylor = [np.zeros((n_p, n_p))] * storage_taylor_order
 
     def gen_lamb_expr_t_recur(i, start_order):
         arg_list = []
@@ -170,11 +170,11 @@ def gen_lamb_expr_t_recur(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)
-            lamb_expr_symb = lamb_expr_symb_deriv
+            lamb_expr_symb = lamb_expr_symb_deriv.subs(var[0], 0)
         else:
             lamb_expr_symb = t_recur.subs(n, i)
 
-        return lamb_expr_symb
+        return sp.lambdify(arg_list, lamb_expr_symb)
 
 
     def gen_lamb_expr_t_exp(i, t_exp_order):
@@ -189,32 +189,71 @@ def gen_lamb_expr_t_exp(i, t_exp_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)
-            lamb_expr_symb = lamb_expr_symb_deriv.subs(var[0], 0)
+            lamb_expr_symb = lamb_expr_symb_deriv
         else:
             lamb_expr_symb = t_exp.subs(n, i)
 
-        return lamb_expr_symb
+        return sp.lambdify(arg_list, lamb_expr_symb)
 
     interactions_off_axis = 0
     for i in range(p+1):
         lamb_expr_t_recur = gen_lamb_expr_t_recur(i, start_order)
-        a1 = [*storage[:-t_recur_order], *coord]
+        a1 = [*storage_taylor[(-t_recur_order):], *coord]
 
         storage.pop(0)
-        storage.append(lamb_expr_t_recur(a1))
+        storage.append(lamb_expr_t_recur(*a1))
 
         lamb_expr_t_exp = gen_lamb_expr_t_exp(i, t_exp_order)
-        a2 = [*storage[:-(t_exp_order+1)], *coord]
+        a2 = [*storage_taylor[-(t_exp_order+1):], *coord]
 
-        interactions_off_axis += lamb_expr_t_exp(a2) * radius**i/math.factorial(i)
+        interactions_off_axis += lamb_expr_t_exp(*a2) * radius**i/math.factorial(i)
 
     ################
+    # Compute True Interactions
+    def generate_true(i):
+        arg_list = []
+        for j in range(ndim):
+            arg_list.append(var[j])
+
+        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)
+        lamb_expr_symb = lamb_expr_symb_deriv
+
+        #print("=============== ORDER = " + str(i))
+        #print(lamb_expr_symb)
+        return sp.lambdify(arg_list, lamb_expr_symb)#, sp.lambdify(arg_list, lamb_expr_symb_deriv)
+
+    interactions_true = 0
+    for i in range(p+1):
+        lamb_expr_true = 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
-    m = 1e10
+    m = 1e5
     mask_on_axis = m*np.abs(coord[0]) >= np.abs(coord[1])
     mask_off_axis = m*np.abs(coord[0]) < np.abs(coord[1])
 
+    print("-------------------------")
+
+    relerr_on = np.abs(interactions_on_axis[mask_on_axis]-interactions_true[mask_on_axis])/np.abs(interactions_on_axis[mask_on_axis])
+    print("MAX ON AXIS ERROR:", np.max(relerr_on))
+    print(np.mean(relerr_on))
+    print("X:", coord[0].reshape(-1)[np.argmax(relerr_on)])
+    print("Y:", coord[1].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:", np.max(relerr_off))
+        print(np.mean(relerr_off))
+        print("X:", coord[0].reshape(-1)[np.argmax(relerr_off)])
+        print("Y:", coord[1].reshape(-1)[np.argmax(relerr_off)])
+
     interactions_total = np.zeros(coord[0].shape)
     interactions_total[mask_on_axis] = interactions_on_axis[mask_on_axis]
     interactions_total[mask_off_axis] = interactions_off_axis[mask_off_axis]

test/test_recurrence.py

@@ -259,7 +259,7 @@ def test_helmholtz_2d_off_axis(deriv_order, 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):
+    for i in range(-exp_range+deriv_order, deriv_order+1):
         approx_deriv = approx_deriv.subs(s(i), ic[i])
 
     rat = coord_dict[var[0]]/coord_dict[var[1]]
@@ -273,7 +273,7 @@ def test_helmholtz_2d_off_axis(deriv_order, exp_order):
     print("RATIO(x0/x1): ", rat)
     #assert relerr <= prederror
 
-test_helmholtz_2d_off_axis(6, 4)
+test_helmholtz_2d_off_axis(5, 4)
 
 
 def test_laplace_2d_off_axis(deriv_order, exp_order):
@@ -319,7 +319,7 @@ def test_laplace_2d_off_axis(deriv_order, exp_order):
 
     exp, exp_range = get_off_axis_expression(laplace2d, exp_order)
     approx_deriv = exp.subs(n, deriv_order)
-    for i in range(exp_range+deriv_order, deriv_order+1):
+    for i in range(-exp_range+deriv_order, deriv_order+1):
         approx_deriv = approx_deriv.subs(s(i), ic[i])
 
     rat = coord_dict[var[0]]/coord_dict[var[1]]

test/test_recurrence_qbx.py

@@ -295,12 +295,15 @@ def _construct_laplace_axis_2d(orders, resolutions):
 
 import matplotlib.pyplot as plt
 orders = [10]
-resolutions = range(200, 800, 200)
+#resolutions = range(200, 800, 200)
+resolutions = [800]
 err_mat = _construct_laplace_axis_2d(orders, resolutions)
+"""
 for i in range(len(orders)):
     plt.plot(resolutions, err_mat[i], label="order ="+str(orders[i]))
 plt.xlabel("Number of Nodes")
 plt.ylabel("Error")
 plt.title("2D Ellipse LP Eval Error")
 plt.legend()
-plt.show()
+plt.show()
+"""