Created plots for recurrenceqbx

Author: hirish99

sumpy/recurrence.py

@@ -374,7 +374,7 @@ def _get_initial_order_off_axis(recurrence):
 
     i = 0
     r_c = recurrence.subs(n, i)
-    while _check_neg_ind(r_c) or r_c == 0:
+    while (_check_neg_ind(r_c) or r_c == 0) or i % 2 != 0:
         i += 1
         r_c = recurrence.subs(n, i)
     return i
@@ -455,15 +455,6 @@ def get_reindexed_and_center_origin_off_axis_recurrence(pde: LinearPDESystemOper
     """
     var = _make_sympy_vec("x", 1)
     r_exp = _move_center_origin_source_arbitrary_expression(pde).subs(var[0], 0)
-
-    var = _make_sympy_vec("x", 2)
-    var_t = _make_sympy_vec("t", 2)
-    g_x_y = sp.log(sp.sqrt((var[0]-var_t[0])**2 + (var[1]-var_t[1])**2))
-    derivs = [sp.diff(g_x_y,
-                      var_t[0], i).subs(var_t[0], 0).subs(var_t[1], 0)
-                      for i in range(8)]
-    
-
     recur_order, recur = reindex_recurrence_relation(r_exp)
     start_order = _get_initial_order_off_axis(recur)
     return start_order, recur_order, recur

sumpy/recurrence_qbx.py

@@ -77,7 +77,7 @@ def _compute_rotated_shifted_coordinates(
 
 # ================ Recurrence LP Eval =================
 def recurrence_qbx_lp(sources, centers, normals, strengths, radius, pde, g_x_y,
-                      ndim, p) -> np.ndarray:
+                      ndim, p, off_axis_start=0) -> np.ndarray:
     r"""
     A function that computes a single-layer potential using a recurrence.
 
@@ -130,8 +130,9 @@ def generate_lamb_expr(i, n_initial):
         #print(lamb_expr_symb)
         return sp.lambdify(arg_list, lamb_expr_symb)#, sp.lambdify(arg_list, lamb_expr_symb_deriv)
 
-    interactions_on_axis = 0
+    
     coord = [cts_r_s[j] for j in range(ndim)]
+    interactions_on_axis = coord[0] * 0
     for i in range(p+1):
         lamb_expr = generate_lamb_expr(i, n_initial)
         a = [*storage, *coord]
@@ -156,6 +157,8 @@ def generate_lamb_expr(i, n_initial):
     t_exp, t_exp_order, _ = get_off_axis_expression(pde, 8)
     storage_taylor_order = max(t_recur_order, t_exp_order+1)
 
+    start_order = max(start_order, order)
+
     storage_taylor = [np.zeros((n_p, n_p))] * storage_taylor_order
     def gen_lamb_expr_t_recur(i, start_order):
         arg_list = []
@@ -210,6 +213,7 @@ def gen_lamb_expr_t_exp(i, t_exp_order, start_order):
 
     ################
     # Compute True Interactions
+    '''
     storage_taylor_true = [np.zeros((n_p, n_p))] * storage_taylor_order
     def generate_true(i):
         arg_list = []
@@ -231,32 +235,33 @@ 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
     m = 100
     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("-------------------------")
+    #print("-------------------------")
 
-    percent_on = np.sum(mask_on_axis)/(mask_on_axis.shape[0]*mask_on_axis.shape[1])
-    percent_off = 1-percent_on
+    #percent_on = np.sum(mask_on_axis)/(mask_on_axis.shape[0]*mask_on_axis.shape[1])
+    #percent_off = 1-percent_on
 
-    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(", percent_on, "):", np.max(relerr_on))
-    print(np.mean(relerr_on))
-    print("X:", coord[0][mask_on_axis].reshape(-1)[np.argmax(relerr_on)])
-    print("Y:", coord[1][mask_on_axis].reshape(-1)[np.argmax(relerr_on)])
+    #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(", percent_on, "):", np.max(relerr_on))
+    #print(np.mean(relerr_on))
+    #print("X:", coord[0][mask_on_axis].reshape(-1)[np.argmax(relerr_on)])
+    #print("Y:", coord[1][mask_on_axis].reshape(-1)[np.argmax(relerr_on)])
 
-    print("-------------------------")
+    #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)])
+    #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)])
 
     interactions_total = np.zeros(coord[0].shape)
     interactions_total[mask_on_axis] = interactions_on_axis[mask_on_axis]

test/investigate_normal_recurrence.ipynb

@@ -86,6 +86,8 @@ def generate_lamb_expr(i, n_initial):
     t_exp, t_exp_order, _ = get_off_axis_expression(pde, 8)
     storage_taylor_order = max(t_recur_order, t_exp_order+1)
 
+    start_order = max(start_order, order)
+
     storage_taylor = [np.zeros((1, n_p))] * storage_taylor_order
     def gen_lamb_expr_t_recur(i, start_order):
         arg_list = []
@@ -95,7 +97,7 @@ def gen_lamb_expr_t_recur(i, start_order):
         for j in range(ndim):
             arg_list.append(var[j])
 
-        if i < start_order+4:
+        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)
@@ -267,12 +269,12 @@ def create_suite_plot(relerr_on, relerr_off, relerr_comb, str_title):
 
 
 #======================== BIHARMONIC 2D ===================================
-interactions_on_axis, interactions_off_axis, interactions_true, interactions_total = produce_error_for_recurrences(mesh_points, biharmonic_pde, g_x_y_biharmonic, 7)
+interactions_on_axis, interactions_off_axis, interactions_true, interactions_total = produce_error_for_recurrences(mesh_points, biharmonic_pde, g_x_y_biharmonic, 12, m=1e2/2)
 
 relerr_on = np.abs((interactions_on_axis-interactions_true)/interactions_true)
 relerr_off = np.abs((interactions_off_axis-interactions_true)/interactions_true)
 relerr_comb = np.abs((interactions_total-interactions_true)/interactions_true)
 
-create_suite_plot(relerr_on, relerr_off+1e-20, relerr_comb+1e-20, "Biharmonic 2D: 7th Order Derivative Evaluation Error $(u_{recur}-u_{sympy})/u_{recur}$")
+create_suite_plot(relerr_on+1e-20, relerr_off+1e-20, relerr_comb+1e-20, "Biharmonic 2D: 8th Order Derivative Evaluation Error $(u_{recur}-u_{sympy})/u_{recur}$")
 
 plt.show()

test/plotting.py

@@ -139,15 +139,14 @@ def test_biharmonic2d():
     var = _make_sympy_vec("x", 2)
     var_t = _make_sympy_vec("t", 2)
     abs_dist = sp.sqrt((var[0]-var_t[0])**2 + (var[1]-var_t[1])**2)
-    w = make_identity_diff_op(2)
-
-    partial_1x = DerivativeIdentifier((4,0), 0)
-    partial_1y = DerivativeIdentifier((0,4), 0)
-    biharmonic_op = {partial_1x: 1, partial_1y: 1}
+    partial_4x = DerivativeIdentifier((4,0), 0)
+    partial_4y = DerivativeIdentifier((0,4), 0)
+    partial_2x2y = DerivativeIdentifier((2,2), 0)
+    biharmonic_op = {partial_4x: 1, partial_4y: 1, partial_2x2y:2}
     list_pde = immutabledict(biharmonic_op)
-
     biharmonic_pde = LinearPDESystemOperator(2, (list_pde,))
-    g_x_y = abs_dist**2 * (sp.log(abs_dist)-1)
+
+    g_x_y = abs_dist**2 * (sp.log(abs_dist))
 
     n_init, _, r = get_reindexed_and_center_origin_on_axis_recurrence(biharmonic_pde)
 
@@ -156,15 +155,15 @@ def test_biharmonic2d():
                                             for i in range(8)]
 
     x_coord = np.random.rand()  # noqa: NPY002
-    y_coord = np.random.rand()  # noqa: NPY002
+    y_coord = np.random.rand() * 1e-3  # noqa: NPY002
     coord_dict = {var[0]: x_coord, var[1]: y_coord}
     derivs = [d.subs(coord_dict) for d in derivs]
 
     n = sp.symbols("n")
     s = sp.Function("s")
 
     # pylint: disable-next=not-callable
-    subs_dict = {s(0): derivs[0], s(1): derivs[1], s(2): derivs[1], s(3): derivs[1]}
+    subs_dict = {s(0): derivs[0], s(1): derivs[1], s(2): derivs[2], s(3): derivs[3]}
     check = []
 
     assert n_init == 4
@@ -175,13 +174,13 @@ def test_biharmonic2d():
         subs_dict[s(i)] = derivs[i]
 
     f2 = sp.lambdify([var[0], var[1]], check[0])
-    assert abs(f2(x_coord, y_coord)) <= 1e-13
+    assert abs(f2(x_coord, y_coord)) <= 1e-11
     f3 = sp.lambdify([var[0], var[1]], check[1])
-    assert abs(f3(x_coord, y_coord)) <= 1e-13
+    assert abs(f3(x_coord, y_coord)) <= 1e-11
     f4 = sp.lambdify([var[0], var[1]], check[2])
-    assert abs(f4(x_coord, y_coord)) <= 1e-13
+    assert abs(f4(x_coord, y_coord)) <= 1e-11
     f5 = sp.lambdify([var[0], var[1]], check[3])
-    assert abs(f5(x_coord, y_coord)) <= 1e-12
+    assert abs(f5(x_coord, y_coord)) <= 1e-11
 
 test_biharmonic2d()
 
@@ -334,6 +333,69 @@ def test_helmholtz_2d_off_axis(deriv_order, exp_order):
     #assert relerr <= prederror
 
 
+def test_helmholtz_2d_off_axis(deriv_order, exp_order):
+    r"""
+    Tests off-axis recurrence code for orders up to 6 laplace2d.
+    """
+    w = make_identity_diff_op(2)
+    helmholtz2d = laplacian(w) + w
+
+    n = sp.symbols("n")
+    s = sp.Function("s")
+
+    var = _make_sympy_vec("x", 2)
+    var_t = _make_sympy_vec("t", 2)
+    abs_dist = sp.sqrt((var[0]-var_t[0])**2 +
+                       (var[1]-var_t[1])**2)
+    g_x_y = abs_dist**2*sp.log(abs_dist)
+    derivs = [sp.diff(g_x_y,
+                      var_t[0], i).subs(var_t[0], 0).subs(var_t[1], 0)
+                                               for i in range(6)]
+    
+    x_coord = 1e-2 * np.random.rand()  # noqa: NPY002
+    y_coord = np.random.rand()  # noqa: NPY002
+    coord_dict = {var[0]: x_coord, var[1]: y_coord}
+    start_order, recur_order, recur = get_reindexed_and_center_origin_off_axis_recurrence(helmholtz2d)
+
+    ic = []
+    #Generate ic
+
+    for i in range(start_order):
+        ic.append(derivs[i].subs(var[0], 0).subs(var[1], coord_dict[var[1]]))
+
+    n = sp.symbols("n")
+    for i in range(start_order, 15):
+        recur_eval = recur.subs(var[0], coord_dict[var[0]]).subs(var[1], coord_dict[var[1]]).subs(n, i)
+        for j in range(i-recur_order, i):
+            recur_eval = recur_eval.subs(s(j), ic[j])
+        ic.append(recur_eval)
+
+    ic = np.array(ic)
+
+    #true_ic = np.array([derivs[i].subs(var[0], 0).subs(var[1], coord_dict[var[1]]) for i in range(15)])
+    
+    #assert np.max(np.abs(ic[::2]-true_ic[::2])/np.abs(true_ic[::2])) < 1e-8
+    #print(np.max(np.abs(ic[::2]-true_ic[::2])/np.abs(true_ic[::2])))
+
+    # CHECK ACCURACY OF EXPRESSION FOR deriv_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])
+    
+    rat = coord_dict[var[0]]/coord_dict[var[1]]
+    if deriv_order + exp_order % 2 == 0:
+        prederror = abs((ic[deriv_order+exp_order+2] * coord_dict[var[0]]**(exp_order+2)/math.factorial(exp_order+2)).evalf())
+    else:
+        prederror = abs((ic[deriv_order+exp_order+1] * coord_dict[var[0]]**(exp_order+1)/math.factorial(exp_order+1)).evalf())
+    print("PREDICTED ERROR: ", prederror)
+    relerr = abs(((approx_deriv - derivs[deriv_order])/derivs[deriv_order]).subs(var[0], coord_dict[var[0]]).subs(var[1], coord_dict[var[1]]).evalf())
+    print("RELATIVE ERROR: ", relerr)
+    print("RATIO(x0/x1): ", rat)
+    #assert relerr <= prederror
+
+
 def test_laplace_2d_off_axis(deriv_order, exp_order):
     r"""
     Tests off-axis recurrence code for orders up to 6 laplace2d.

test/test_recurrence.py

@@ -98,7 +98,7 @@ def _qbx_lp_general(knl, sources, targets, centers, radius,
 
 def _create_ellipse(n_p):
     h = 9.688 / n_p
-    radius = 7*h * 1/40
+    radius = (h * 1/4) * 1.11
     t = np.linspace(0, 2 * np.pi, n_p, endpoint=False)
 
     unit_circle_param = np.exp(1j * t)
@@ -294,15 +294,21 @@ def _construct_laplace_axis_2d(orders, resolutions):
     return err
 
 import matplotlib.pyplot as plt
-orders = [10, 16]
+orders = [6, 7, 8, 9]
 #resolutions = range(200, 800, 200)
-resolutions = [400, 800, 1200]
+resolutions = [200, 300]
 err_mat = _construct_laplace_axis_2d(orders, resolutions)
 
+fig, ax = plt.subplots(figsize = (6, 6))
+ax.set_yscale("log")
+
+orders_fake = [12, 14, 16, 18]
+
 for i in range(len(orders)):
-    plt.scatter(resolutions, err_mat[i], label="order QBX="+str(orders[i]))
-plt.xlabel("Number of Nodes")
-plt.ylabel("Relative Error")
-plt.title("2D Ellipse LP Eval Error (m=10)")
-plt.legend()
+    ax.scatter(9.68845/np.array(resolutions), np.array(err_mat[i]), label="$p_{QBX}$="+str(orders_fake[i]))
+ax.set_xlabel("Mesh Resolution ($h$)")
+ax.set_ylabel("Relative Error ($L_{\infty}$)")
+ax.set_title("Laplace 2D: Ellipse SLP Boundary Evaluation Error $(u_{qbxrec}-u_{qbx})/u_{qbx}$")
+ax.legend()
 plt.show()
+#$m=100$, $r=0.28h$