3434"""
3535
3636import math
37- from typing import Sequence
37+ from collections . abc import Sequence
3838
3939import numpy as np
4040import sympy as sp
@@ -102,7 +102,9 @@ def recurrence_qbx_lp(sources, centers, normals, strengths, radius, pde, g_x_y,
102102 var_t = _make_sympy_vec ("t" , ndim )
103103
104104 # ------------ 5. Compute recurrence
105- n_initial , order , recurrence = get_reindexed_and_center_origin_on_axis_recurrence (pde )
105+ n_initial , order , recurrence = (
106+ get_reindexed_and_center_origin_on_axis_recurrence (pde )
107+ )
106108
107109 # ------------ 6. Set order p = 5
108110 n_p = sources .shape [1 ]
@@ -126,11 +128,10 @@ def generate_lamb_expr(i, n_initial):
126128 lamb_expr_symb = lamb_expr_symb_deriv
127129 else :
128130 lamb_expr_symb = recurrence .subs (n , i )
129- #print("=============== ORDER = " + str(i))
130- #print(lamb_expr_symb)
131- return sp .lambdify (arg_list , lamb_expr_symb )#, sp.lambdify(arg_list, lamb_expr_symb_deriv)
131+ # print("=============== ORDER = " + str(i))
132+ # print(lamb_expr_symb)
133+ return sp .lambdify (arg_list , lamb_expr_symb )
132134
133-
134135 coord = [cts_r_s [j ] for j in range (ndim )]
135136 interactions_on_axis = coord [0 ] * 0
136137 for i in range (p + 1 ):
@@ -142,24 +143,29 @@ def generate_lamb_expr(i, n_initial):
142143 s_new_true = true_lamb_expr(*a)
143144 arg_max = np.argmax(abs(s_new-s_new_true)/abs(s_new_true))
144145 print((s_new-s_new_true).reshape(-1)[arg_max]/s_new_true.reshape(-1)[arg_max])
145- print("x:", coord[0].reshape(-1)[arg_max], "y:", coord[1].reshape(-1)[arg_max],
146- "s_recur:", s_new.reshape(-1)[arg_max], "s_true:", s_new_true.reshape(-1)[arg_max], "order: ", i)
146+ print("x:", coord[0].reshape(-1)[arg_max],
147+ "y:", coord[1].reshape(-1)[arg_max],
148+ "s_recur:", s_new.reshape(-1)[arg_max],
149+ "s_true:", s_new_true.reshape(-1)[arg_max],
150+ "order: ", i)
147151 """
148152
149153 interactions_on_axis += s_new * radius ** i / math .factorial (i )
150154
151155 storage .pop (0 )
152156 storage .append (s_new )
153157
154-
155- ### NEW CODE - COMPUTE OFF AXIS INTERACTIONS
156- start_order , t_recur_order , t_recur = get_reindexed_and_center_origin_off_axis_recurrence (pde )
158+ # NEW CODE - COMPUTE OFF AXIS INTERACTIONS
159+ start_order , t_recur_order , t_recur = (
160+ get_reindexed_and_center_origin_off_axis_recurrence (pde )
161+ )
157162 t_exp , t_exp_order , _ = get_off_axis_expression (pde , 8 )
158163 storage_taylor_order = max (t_recur_order , t_exp_order + 1 )
159164
160165 start_order = max (start_order , order )
161166
162167 storage_taylor = [np .zeros ((n_p , n_p ))] * storage_taylor_order
168+
163169 def gen_lamb_expr_t_recur (i , start_order ):
164170 arg_list = []
165171 for j in range (t_recur_order , 0 , - 1 ):
@@ -178,7 +184,6 @@ def gen_lamb_expr_t_recur(i, start_order):
178184
179185 return sp .lambdify (arg_list , lamb_expr_symb )
180186
181-
182187 def gen_lamb_expr_t_exp (i , t_exp_order , start_order ):
183188 arg_list = []
184189 for j in range (t_exp_order , - 1 , - 1 ):
@@ -197,7 +202,6 @@ def gen_lamb_expr_t_exp(i, t_exp_order, start_order):
197202
198203 return sp .lambdify (arg_list , lamb_expr_symb )
199204
200-
201205 interactions_off_axis = 0
202206 for i in range (p + 1 ):
203207 lamb_expr_t_recur = gen_lamb_expr_t_recur (i , start_order )
@@ -213,7 +217,7 @@ def gen_lamb_expr_t_exp(i, t_exp_order, start_order):
213217
214218 ################
215219 # Compute True Interactions
216- '''
220+ """
217221 storage_taylor_true = [np.zeros((n_p, n_p))] * storage_taylor_order
218222 def generate_true(i):
219223 arg_list = []
@@ -225,52 +229,65 @@ def generate_true(i):
225229 lamb_expr_symb_deriv = lamb_expr_symb_deriv.subs(var_t[j], 0)
226230 lamb_expr_symb = lamb_expr_symb_deriv
227231
228- #print("=============== ORDER = " + str(i))
229- #print(lamb_expr_symb)
230- return sp.lambdify(arg_list, lamb_expr_symb)#, sp.lambdify(arg_list, lamb_expr_symb_deriv)
232+ # print("=============== ORDER = " + str(i))
233+ # print(lamb_expr_symb)
234+ return sp.lambdify(arg_list, lamb_expr_symb)
231235
232236 interactions_true = 0
233237 for i in range(p+1):
234238 lamb_expr_true = generate_true(i)
235239 a4 = [*coord]
236240 s_new_true = lamb_expr_true(*a4)
237241 interactions_true += s_new_true * radius**i/math.factorial(i)
238- '''
242+ """
239243 ###############
240244
241- #slope of line y = mx
245+ # slope of line y = mx
242246 m = 100
243247 mask_on_axis = m * np .abs (coord [0 ]) >= np .abs (coord [1 ])
244248 mask_off_axis = m * np .abs (coord [0 ]) < np .abs (coord [1 ])
245249
246- #print("-------------------------")
247-
248- #percent_on = np.sum(mask_on_axis)/(mask_on_axis.shape[0]*mask_on_axis.shape[1])
249- #percent_off = 1-percent_on
250-
251- #relerr_on = np.abs(interactions_on_axis[mask_on_axis]-interactions_true[mask_on_axis])/np.abs(interactions_on_axis[mask_on_axis])
252- #print("MAX ON AXIS ERROR(", percent_on, "):", np.max(relerr_on))
253- #print(np.mean(relerr_on))
254- #print("X:", coord[0][mask_on_axis].reshape(-1)[np.argmax(relerr_on)])
255- #print("Y:", coord[1][mask_on_axis].reshape(-1)[np.argmax(relerr_on)])
256-
257- #print("-------------------------")
258-
259- #if np.any(mask_off_axis):
260- # relerr_off = np.abs(interactions_off_axis[mask_off_axis]-interactions_true[mask_off_axis])/np.abs(interactions_off_axis[mask_off_axis])
261- # print("MAX OFF AXIS ERROR(", percent_off, "):", np.max(relerr_off))
262- # print(np.mean(relerr_off))
263- # print("X:", coord[0][mask_off_axis].reshape(-1)[np.argmax(relerr_off)])
264- # print("Y:", coord[1][mask_off_axis].reshape(-1)[np.argmax(relerr_off)])
250+ # print("-------------------------")
251+
252+ # percent_on = np.sum(mask_on_axis) / (
253+ # mask_on_axis.shape[0]*mask_on_axis.shape[1])
254+ # percent_off = 1-percent_on
255+
256+ # relerr_on = (
257+ # np.abs(
258+ # interactions_on_axis[mask_on_axis]
259+ # - interactions_true[mask_on_axis])
260+ # / np.abs(interactions_on_axis[mask_on_axis]))
261+ # print("MAX ON AXIS ERROR(", percent_on, "):", np.max(relerr_on))
262+ # print(np.mean(relerr_on))
263+ # print("X:", coord[0][mask_on_axis].reshape(-1)[np.argmax(relerr_on)])
264+ # print("Y:", coord[1][mask_on_axis].reshape(-1)[np.argmax(relerr_on)])
265+
266+ # print("-------------------------")
267+
268+ # if np.any(mask_off_axis):
269+ # relerr_off = (
270+ # np.abs(
271+ # interactions_off_axis[mask_off_axis]
272+ # - interactions_true[mask_off_axis])
273+ # / np.abs(interactions_off_axis[mask_off_axis]))
274+ # print("MAX OFF AXIS ERROR(", percent_off, "):",
275+ # np.max(relerr_off))
276+ # print(np.mean(relerr_off))
277+ # print("X:",
278+ # coord[0][mask_off_axis].reshape(-1)[np.argmax(relerr_off)])
279+ # print("Y:",
280+ # coord[1][mask_off_axis].reshape(-1)[np.argmax(relerr_off)])
265281
266282 interactions_total = np .zeros (coord [0 ].shape )
267283 interactions_total [mask_on_axis ] = interactions_on_axis [mask_on_axis ]
268284 interactions_total [mask_off_axis ] = interactions_off_axis [mask_off_axis ]
269285
270286 exp_res = (interactions_total * strengths [None , :]).sum (axis = 1 )
271- #exp_res_true = (interactions_true * strengths[None, :]).sum(axis=1)
287+ # exp_res_true = (interactions_true * strengths[None, :]).sum(axis=1)
272288
273- #relerr_total = np.max(np.abs(exp_res-exp_res_true)/np.abs(exp_res_true))
274- #print("OVERALL ERROR:", relerr_total)
289+ # relerr_total = np.max(
290+ # np.abs(exp_res-exp_res_true)/np.abs(exp_res_true))
291+ # print("OVERALL ERROR:", relerr_total)
275292
276- return exp_res
293+ return exp_res
0 commit comments