diff --git a/setup.cfg b/setup.cfg index 77f81d3bfbff0b437d7d65c0bcda8064887aa084..65a7a4a38d2c90be87ca3fe40e83f2e856869d45 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,3 +1,6 @@ [flake8] ignore = E126,E127,E128,E123,E226,E241,E242,E265,E402,W503 max-line-length=85 + +[pycodestyle] +max-line-length=85 diff --git a/setup.py b/setup.py index 5004b8d061f28813b3c3d673b86ea8e9c185eab9..8b037de997fa0153013ace3061e836840ce2b883 100644 --- a/setup.py +++ b/setup.py @@ -120,6 +120,7 @@ def main(): "pytential", "scipy", "sumpy", + "mpmath", ], extras_require={ "gmsh_support": ["gmsh_interop"], diff --git a/volumential/droste.py b/volumential/droste.py index 45845df53b4c281c023d31e05c1c8c8c4e34c0f8..382fa2fd647b53006ae4b43dfb33f66f4a7c980a 100644 --- a/volumential/droste.py +++ b/volumential/droste.py @@ -63,9 +63,16 @@ class DrosteBase(KernelCacheWrapper): .. attribute:: interaction_case_scls The relative sizes of the target box for each case. + + .. attribute:: special_radial_quadrature + + If True, the radial direction uses a different quadrature rule. + The radial direction is identified with the condition + ``ibrick_axis == iaxis`` and always uses iname ``q0``. """ - def __init__(self, integral_knl, quad_order, case_vecs, n_brick_quad_points=15): + def __init__(self, integral_knl, quad_order, case_vecs, n_brick_quad_points, + special_radial_quadrature, nradial_quad_points): from sumpy.kernel import Kernel if integral_knl is not None: @@ -101,6 +108,14 @@ class DrosteBase(KernelCacheWrapper): self.ntgt_points = quad_order self.nquad_points = n_brick_quad_points + if special_radial_quadrature and (self.dim == 1): + raise ValueError( + "please only use normal quadrature nodes for 1D quadrature") + + self.special_radial_quadrature = special_radial_quadrature + self.nradial_quad_points = nradial_quad_points + self.brick_quadrature_kwargs = self.make_brick_quadrature_kwargs() + if self.dim is not None: self.n_q_points = quad_order ** self.dim else: @@ -157,6 +172,56 @@ class DrosteBase(KernelCacheWrapper): + [pwaffs[var].lt_set(ubound_pwaff) for var in variables], ) + def make_brick_quadrature_kwargs(self): + """Produce 1D quadrature formulae used for each brick. + The rules returned are defined over [0, 1]. + """ + # sps.legendre blows up easily at high order + import scipy.special as sps + # legendre_nodes, _, legendre_weights = sps.legendre( + # nquad_points).weights.T + legendre_nodes, legendre_weights = sps.p_roots(self.nquad_points) + legendre_nodes = legendre_nodes * 0.5 + 0.5 + legendre_weights = legendre_weights * 0.5 + + if self.special_radial_quadrature: + # Since there is no exact degree -> nqpts map for tanh-sinh rule, + # we find the largest rule whose size dose not exceed the parameter + # nradial_quad_points + import mpmath + mp_ctx = mpmath.mp + mp_ctx.dps = 20 # decimal precision + mp_ctx.pretty = True + ts = mpmath.calculus.quadrature.TanhSinh(mp_ctx) + prec = int(np.log(10) / np.log(2) * mp_ctx.dps) # bits of precision + + for deg in range(1, 100): + nodes = ts.calc_nodes(degree=deg, prec=prec) + if len(nodes) >= self.nradial_quad_points: + self.nradial_quad_points = len(nodes) + break + + # extract quadrature formula over [-1, 1], note the 0.5**level scaling + ts_nodes = np.array([p[0] for p in nodes], dtype=np.float64) + ts_weights = np.array( + [p[1] * 2 * mp_ctx.power(0.5, deg) for p in nodes], + dtype=np.float64) + if deg == 1: + ts_weights *= 0.5 # peculiar scaling when deg=1 + + # map to [0, 1] + ts_nodes = ts_nodes * 0.5 + 0.5 + ts_weights *= 0.5 + + return dict(quadrature_nodes=legendre_nodes, + quadrature_weights=legendre_weights, + radial_quadrature_nodes=ts_nodes, + radial_quadrature_weights=ts_weights) + + else: + return dict(quadrature_nodes=legendre_nodes, + quadrature_weights=legendre_weights) + def get_sumpy_kernel_insns(self): # get sumpy kernel insns from sumpy.symbolic import make_sym_vector @@ -257,19 +322,11 @@ class DrosteBase(KernelCacheWrapper): "<> knl_val_post = knl_val {id=pp_kval}" ) - resknl = resknl.replace( - "PROD_QUAD_WEIGHT", - " * ".join( - [ - "quadrature_weights[QID]".replace("QID", qvar) - for qvar in self.quad_vars - ] - ), - ) - if self.dim == 1: + resknl = resknl.replace("TPLTGT_ASSIGNMENT", """target_nodes[t0]""") resknl = resknl.replace("QUAD_PT_ASSIGNMENT", """quadrature_nodes[q0]""") + resknl = resknl.replace("PROD_QUAD_WEIGHT", """quadrature_weights[q0]""") resknl = resknl.replace("DENSITY_VAL_ASSIGNMENT", """basis_eval0""") resknl = resknl.replace( "PREPARE_BASIS_VALS", @@ -284,10 +341,23 @@ class DrosteBase(KernelCacheWrapper): "TPLTGT_ASSIGNMENT", """if(iaxis == 0, target_nodes[t0], target_nodes[t1])""", ) - resknl = resknl.replace( - "QUAD_PT_ASSIGNMENT", - """if(iaxis == 0, quadrature_nodes[q0], quadrature_nodes[q1])""", - ) + if self.special_radial_quadrature: + resknl = resknl.replace( + "QUAD_PT_ASSIGNMENT", + """if(iaxis == ibrick_axis, + radial_quadrature_nodes[q0], quadrature_nodes[q1])""", + ) + resknl = resknl.replace( + "PROD_QUAD_WEIGHT", + """radial_quadrature_weights[q0] * quadrature_weights[q1]""") + else: + resknl = resknl.replace( + "QUAD_PT_ASSIGNMENT", + """if(iaxis == 0, quadrature_nodes[q0], quadrature_nodes[q1])""", + ) + resknl = resknl.replace( + "PROD_QUAD_WEIGHT", + """quadrature_weights[q0] * quadrature_weights[q1]""") resknl = resknl.replace( "DENSITY_VAL_ASSIGNMENT", """basis_eval0 * basis_eval1""" ) @@ -305,11 +375,28 @@ class DrosteBase(KernelCacheWrapper): """if(iaxis == 0, target_nodes[t0], if( iaxis == 1, target_nodes[t1], target_nodes[t2]))""", ) - resknl = resknl.replace( - "QUAD_PT_ASSIGNMENT", - """if(iaxis == 0, quadrature_nodes[q0], if( - iaxis == 1, quadrature_nodes[q1], quadrature_nodes[q2]))""", - ) + if self.special_radial_quadrature: + # depending on ibrick_axis, axes could be + # [q0, q1, q2], [q1, q0, q2], or [q1, q2, q0] + resknl = resknl.replace( + "QUAD_PT_ASSIGNMENT", + """if(iaxis == ibrick_axis, radial_quadrature_nodes[q0], if( + (iaxis == 0) or ((iaxis == 1) and (ibrick_axis == 0)), + quadrature_nodes[q1], + quadrature_nodes[q2]))""", + ) + resknl = resknl.replace( + "PROD_QUAD_WEIGHT", + """radial_w[q0] * w[q1] * w[q2]""".replace( + "w", "quadrature_weights")) + else: + resknl = resknl.replace( + "QUAD_PT_ASSIGNMENT", + """if(iaxis == 0, quadrature_nodes[q0], if( + iaxis == 1, quadrature_nodes[q1], quadrature_nodes[q2]))""") + resknl = resknl.replace( + "PROD_QUAD_WEIGHT", + """w[q0] * w[q1] * w[q2]""".replace("w", "quadrature_weights")) resknl = resknl.replace( "DENSITY_VAL_ASSIGNMENT", """basis_eval0 * basis_eval1 * basis_eval2""" @@ -566,8 +653,7 @@ class DrosteBase(KernelCacheWrapper): * knl_scaling ) {id=result,dep=jac:mpoint:density:finish_kval} end - """ - ) + """) ] def get_target_points(self, queue=None): @@ -630,9 +716,13 @@ class DrosteFull(DrosteBase): Build the full table directly. """ - def __init__(self, integral_knl, quad_order, case_vecs, n_brick_quad_points=50): + def __init__(self, integral_knl, quad_order, case_vecs, + n_brick_quad_points=50, + special_radial_quadrature=False, + nradial_quad_points=None): super(DrosteFull, self).__init__( - integral_knl, quad_order, case_vecs, n_brick_quad_points) + integral_knl, quad_order, case_vecs, n_brick_quad_points, + special_radial_quadrature, nradial_quad_points) self.name = "DrosteFull" def make_loop_domain(self): @@ -641,9 +731,16 @@ class DrosteFull(DrosteBase): basis_vars = self.make_basis_vars() basis_eval_vars = self.make_basis_eval_vars() pwaffs = self.make_pwaffs() # noqa: F841 + if self.special_radial_quadrature: + quad_vars_subdomain = self.make_brick_domain( + quad_vars[1:], self.nquad_points) & self.make_brick_domain( + quad_vars[:1], self.nradial_quad_points) + else: + quad_vars_subdomain = self.make_brick_domain( + quad_vars, self.nquad_points) self.loop_domain = ( self.make_brick_domain(tgt_vars, self.ntgt_points) - & self.make_brick_domain(quad_vars, self.nquad_points) + & quad_vars_subdomain & self.make_brick_domain(basis_vars, self.nfunctions) & self.make_brick_domain(basis_eval_vars, self.nfunctions) & self.make_brick_domain(["iside"], 2) @@ -665,7 +762,8 @@ class DrosteFull(DrosteBase): loopy_knl = lp.make_kernel( # NOQA [domain], - self.get_kernel_code() + self.get_sumpy_kernel_eval_insns(), + self.get_kernel_code() + + self.get_sumpy_kernel_eval_insns(), [ lp.ValueArg("alpha", np.float64), lp.ValueArg("n_cases, nfunctions, quad_order, dim", np.int32), @@ -708,10 +806,10 @@ class DrosteFull(DrosteBase): def call_loopy_kernel(self, queue, **kwargs): """ - :arg source_box_extent - :arg alpha - :arg nlevels - :arg extra_kernel_kwargs + :param source_box_extent: + :param alpha: + :param nlevels: + :param extra_kernel_kwargs: """ if "source_box_extent" in kwargs: @@ -733,6 +831,9 @@ class DrosteFull(DrosteBase): if "nlevels" in kwargs: nlevels = kwargs["nlevels"] assert nlevels > 0 + if nlevels > 1 and self.special_radial_quadrature: + logger.warn("When using tanh-sinh quadrature in the radial " + "direction, it is often best to use a single level.") else: # Single level is equivalent to Duffy transform nlevels = 1 @@ -761,37 +862,19 @@ class DrosteFull(DrosteBase): assert len(q_points) == self.ntgt_points ** self.dim t = np.array([pt[-1] for pt in q_points[: self.ntgt_points]]) - # quad formula for each brick (normalized to [0,1]) - # sps.legendre blows up easily at high order - import scipy.special as sps - - # legendre_nodes, _, legendre_weights = sps.legendre( - # nquad_points).weights.T - legendre_nodes, legendre_weights = sps.p_roots(self.nquad_points) - legendre_nodes = legendre_nodes * 0.5 + 0.5 - legendre_weights = legendre_weights * 0.5 - knl = self.get_cached_optimized_kernel() evt, res = knl( - queue, - alpha=alpha, - result=result_array, + queue, alpha=alpha, result=result_array, root_brick=root_brick, target_nodes=t.astype(np.float64, copy=True), - quadrature_nodes=legendre_nodes, - quadrature_weights=legendre_weights, interaction_case_vecs=self.interaction_case_vecs.astype( - np.float64, copy=True - ), + np.float64, copy=True), interaction_case_scls=self.interaction_case_scls.reshape(-1).astype( - np.float64, copy=True - ), - n_cases=self.ncases, - nfunctions=self.nfunctions, - quad_order=self.ntgt_points, - nlevels=nlevels, - **extra_kernel_kwargs - ) + np.float64, copy=True), + n_cases=self.ncases, nfunctions=self.nfunctions, + quad_order=self.ntgt_points, nlevels=nlevels, + **self.brick_quadrature_kwargs, + **extra_kernel_kwargs) cheb_table = res["result"] @@ -804,6 +887,7 @@ class DrosteFull(DrosteBase): self.integral_knl.__str__(), "quad_order-" + str(self.ntgt_points), "brick_order-" + str(self.nquad_points), + "brick_radial_order-" + str(self.nradial_quad_points), ) def __call__(self, queue, **kwargs): @@ -844,9 +928,12 @@ class DrosteReduced(DrosteBase): case_vecs=None, n_brick_quad_points=50, knl_symmetry_tags=None, + special_radial_quadrature=False, + nradial_quad_points=None, ): super(DrosteReduced, self).__init__( - integral_knl, quad_order, case_vecs, n_brick_quad_points) + integral_knl, quad_order, case_vecs, n_brick_quad_points, + special_radial_quadrature, nradial_quad_points) from volumential.list1_symmetry import CaseVecReduction @@ -882,8 +969,17 @@ class DrosteReduced(DrosteBase): basis_vars = self.make_basis_vars() basis_eval_vars = self.make_basis_eval_vars() pwaffs = self.make_pwaffs() + + if self.special_radial_quadrature: + quad_vars_subdomain = self.make_brick_domain( + quad_vars[1:], self.nquad_points) & self.make_brick_domain( + quad_vars[:1], self.nradial_quad_points) + else: + quad_vars_subdomain = self.make_brick_domain( + quad_vars, self.nquad_points) + loop_domain_common_parts = ( - self.make_brick_domain(quad_vars, self.nquad_points) + quad_vars_subdomain & self.make_brick_domain(basis_vars, self.nfunctions) & self.make_brick_domain(basis_eval_vars, self.nfunctions) & self.make_brick_domain(["iside"], 2) @@ -1260,16 +1356,6 @@ class DrosteReduced(DrosteBase): ] ) - # quad formula for each brick (normalized to [0,1]) - # sps.legendre blows up easily at high order - import scipy.special as sps - - # legendre_nodes, _, legendre_weights = sps.legendre( - # nquad_points).weights.T - legendre_nodes, legendre_weights = sps.p_roots(self.nquad_points) - legendre_nodes = legendre_nodes * 0.5 + 0.5 - legendre_weights = legendre_weights * 0.5 - self.get_kernel_id = 0 try: delattr(self, "_memoize_dic_get_cached_optimized_kernel") @@ -1277,21 +1363,15 @@ class DrosteReduced(DrosteBase): pass knl = self.get_cached_optimized_kernel() evt, res = knl( - queue, - alpha=alpha, - result=result_array, + queue, alpha=alpha, result=result_array, root_brick=root_brick, target_nodes=t.astype(np.float64, copy=True), - quadrature_nodes=legendre_nodes, - quadrature_weights=legendre_weights, interaction_case_vecs=base_case_vec.astype(np.float64, copy=True), interaction_case_scls=base_case_scl.astype(np.float64, copy=True), - n_cases=1, - nfunctions=self.nfunctions, - quad_order=self.ntgt_points, - nlevels=nlevels, - **extra_kernel_kwargs - ) + n_cases=1, nfunctions=self.nfunctions, + quad_order=self.ntgt_points, nlevels=nlevels, + **extra_kernel_kwargs, **self.brick_quadrature_kwargs) + raw_cheb_table_case = res["result"] self.get_kernel_id = 1 @@ -1443,6 +1523,7 @@ class DrosteReduced(DrosteBase): self.integral_knl.__str__(), "quad_order-" + str(self.ntgt_points), "brick_order-" + str(self.nquad_points), + "brick_radial_order-" + str(self.nradial_quad_points), "case-" + str(self.current_base_case), "kernel_id-" + str(self.get_kernel_id), "symmetry-" + symmetry_info, @@ -1511,20 +1592,17 @@ class InverseDrosteReduced(DrosteReduced): """ - def __init__( - self, - integral_knl, - quad_order, - case_vecs, - n_brick_quad_points=50, - knl_symmetry_tags=None, - auto_windowing=True): + def __init__(self, integral_knl, quad_order, case_vecs, + n_brick_quad_points=50, knl_symmetry_tags=None, + special_radial_quadrature=False, nradial_quad_points=None, + auto_windowing=True): """ :param auto_windowing: auto-detect window radius. """ super(InverseDrosteReduced, self).__init__( integral_knl, quad_order, case_vecs, - n_brick_quad_points, knl_symmetry_tags) + n_brick_quad_points, special_radial_quadrature, + nradial_quad_points, knl_symmetry_tags) self.auto_windowing = auto_windowing def get_cache_key(self): @@ -1534,9 +1612,9 @@ class InverseDrosteReduced(DrosteReduced): self.integral_knl.__str__(), "quad_order-" + str(self.ntgt_points), "brick_order-" + str(self.nquad_points), + "brick_radial_order-" + str(self.nradial_quad_points), "case-" + str(self.current_base_case), - "kernel_id-" + str(self.get_kernel_id), - ) + "kernel_id-" + str(self.get_kernel_id),) def codegen_basis_tgt_eval(self, iaxis): """Generate instructions to evaluate Chebyshev polynomial basis @@ -2053,16 +2131,6 @@ class InverseDrosteReduced(DrosteReduced): ] ) - # quad formula for each brick (normalized to [0,1]) - # sps.legendre blows up easily at high order - import scipy.special as sps - - # legendre_nodes, _, legendre_weights = sps.legendre( - # nquad_points).weights.T - legendre_nodes, legendre_weights = sps.p_roots(self.nquad_points) - legendre_nodes = legendre_nodes * 0.5 + 0.5 - legendre_weights = legendre_weights * 0.5 - # --------- call kernel 0 ---------- self.get_kernel_id = 0 try: @@ -2077,15 +2145,13 @@ class InverseDrosteReduced(DrosteReduced): result=result_array_0, root_brick=root_brick, target_nodes=t.astype(np.float64, copy=True), - quadrature_nodes=legendre_nodes, - quadrature_weights=legendre_weights, interaction_case_vecs=base_case_vec.astype(np.float64, copy=True), interaction_case_scls=base_case_scl.astype(np.float64, copy=True), n_cases=1, nfunctions=self.nfunctions, quad_order=self.ntgt_points, nlevels=nlevels, - **extra_kernel_kwargs + **extra_kernel_kwargs, **self.brick_quadrature_kwargs ) # --------- call kernel 1 ---------- @@ -2102,15 +2168,13 @@ class InverseDrosteReduced(DrosteReduced): result=result_array_1, root_brick=root_brick, target_nodes=t.astype(np.float64, copy=True), - quadrature_nodes=legendre_nodes, - quadrature_weights=legendre_weights, interaction_case_vecs=base_case_vec.astype(np.float64, copy=True), interaction_case_scls=base_case_scl.astype(np.float64, copy=True), n_cases=1, nfunctions=self.nfunctions, quad_order=self.ntgt_points, nlevels=nlevels, - **extra_kernel_kwargs + **extra_kernel_kwargs, **self.brick_quadrature_kwargs ) # --------- call kernel 2 ---------- diff --git a/volumential/nearfield_potential_table.py b/volumential/nearfield_potential_table.py index cb0bcbdaa74a257020fb2a6d60cc9bf4b04746bb..8cef4d2cf84de2e55e9d91c22124dbdfaaae9a27 100644 --- a/volumential/nearfield_potential_table.py +++ b/volumential/nearfield_potential_table.py @@ -800,6 +800,39 @@ class NearFieldInteractionTable(object): # {{{ build table via adding up a Droste of bricks + def get_droste_table_builder(self, n_brick_quad_points, + special_radial_brick_quadrature, + nradial_brick_quad_points, + use_symmetry=False, + knl_symmetry_tags=None): + if self.inverse_droste: + from volumential.droste import InverseDrosteReduced + + drf = InverseDrosteReduced( + self.integral_knl, self.quad_order, self.interaction_case_vecs, + n_brick_quad_points, knl_symmetry_tags, auto_windowing=False, + special_radial_quadrature=special_radial_brick_quadrature, + nradial_quad_points=nradial_brick_quad_points) + + else: + if not use_symmetry: + from volumential.droste import DrosteFull + + drf = DrosteFull( + self.integral_knl, self.quad_order, + self.interaction_case_vecs, n_brick_quad_points, + special_radial_quadrature=special_radial_brick_quadrature, + nradial_quad_points=nradial_brick_quad_points) + else: + from volumential.droste import DrosteReduced + + drf = DrosteReduced( + self.integral_knl, self.quad_order, self.interaction_case_vecs, + n_brick_quad_points, knl_symmetry_tags, + special_radial_quadrature=special_radial_brick_quadrature, + nradial_quad_points=nradial_brick_quad_points) + return drf + def build_table_via_droste_bricks( self, n_brick_quad_points=50, @@ -824,20 +857,15 @@ class NearFieldInteractionTable(object): else: nlev = 1 - # extra_kernel_kwargs = {} - # if "extra_kernel_kwargs" in kwargs: - # extra_kernel_kwargs = kwargs["extra_kernel_kwargs"] - - cheb_coefs = [ - self.get_mode_cheb_coeffs(mid, self.quad_order) - for mid in range(self.n_q_points) - ] - - # {{{ get the table builder - - if self.inverse_droste: - from volumential.droste import InverseDrosteReduced + if "special_radial_brick_quadrature" in kwargs: + special_radial_brick_quadrature = kwargs.pop( + "special_radial_brick_quadrature") + nradial_brick_quad_points = kwargs.pop("nradial_brick_quad_points") + else: + special_radial_brick_quadrature = False + nradial_brick_quad_points = None + if use_symmetry: if "knl_symmetry_tags" in kwargs: knl_symmetry_tags = kwargs["knl_symmetry_tags"] else: @@ -850,67 +878,31 @@ class NearFieldInteractionTable(object): ) knl_symmetry_tags = None - drf = InverseDrosteReduced( - self.integral_knl, - self.quad_order, - self.interaction_case_vecs, - n_brick_quad_points, - knl_symmetry_tags, - auto_windowing=False - ) - - else: - if not use_symmetry: - from volumential.droste import DrosteFull - - drf = DrosteFull( - self.integral_knl, - self.quad_order, - self.interaction_case_vecs, - n_brick_quad_points, - ) - else: - from volumential.droste import DrosteReduced - - if "knl_symmetry_tags" in kwargs: - knl_symmetry_tags = kwargs["knl_symmetry_tags"] - else: - # Maximum symmetry by default - logger.warning( - "use_symmetry is set to True, but knl_symmetry_tags is not " - "set. Using the default maximum symmetry. (Using maximum " - "symmetry for some kernels (e.g. derivatives of " - "LaplaceKernel will yield incorrect results)." - ) - knl_symmetry_tags = None - - drf = DrosteReduced( - self.integral_knl, - self.quad_order, - self.interaction_case_vecs, - n_brick_quad_points, - knl_symmetry_tags, - ) + # extra_kernel_kwargs = {} + # if "extra_kernel_kwargs" in kwargs: + # extra_kernel_kwargs = kwargs["extra_kernel_kwargs"] - # }}} End get the table builder + cheb_coefs = [ + self.get_mode_cheb_coeffs(mid, self.quad_order) + for mid in range(self.n_q_points) + ] # compute an initial table - data0 = drf( - queue, - source_box_extent=self.source_box_extent, - alpha=alpha, - nlevels=nlev, - # extra_kernel_kwargs=extra_kernel_kwargs, - cheb_coefs=cheb_coefs, - **kwargs - ) + drf = self.get_droste_table_builder( + n_brick_quad_points, + special_radial_brick_quadrature, nradial_brick_quad_points, + use_symmetry, knl_symmetry_tags) + data0 = drf(queue, source_box_extent=self.source_box_extent, + alpha=alpha, nlevels=nlev, + # extra_kernel_kwargs=extra_kernel_kwargs, + cheb_coefs=cheb_coefs, **kwargs) # {{{ adaptively determine number of levels resid = -1 missing_measure = 1 if adaptive_level: - table_tol = np.finfo(self.dtype).eps * 64 + table_tol = np.finfo(self.dtype).eps * 256 # 5e-14 for float64 logger.warn("Searching for nlevels since adaptive_level=True") while True: @@ -926,15 +918,10 @@ class NearFieldInteractionTable(object): break nlev = nlev + 1 - data1 = drf( - queue, - source_box_extent=self.source_box_extent, - alpha=alpha, - nlevels=nlev, - # extra_kernel_kwargs=extra_kernel_kwargs, - cheb_coefs=cheb_coefs, - **kwargs - ) + data1 = drf(queue, source_box_extent=self.source_box_extent, + alpha=alpha, nlevels=nlev, + # extra_kernel_kwargs=extra_kernel_kwargs, + cheb_coefs=cheb_coefs, **kwargs) resid = np.max(np.abs(data1 - data0)) / np.max(np.abs(data1)) data0 = data1 @@ -954,24 +941,36 @@ class NearFieldInteractionTable(object): if resid >= table_tol: logger.warn("Adaptive level refinement failed to converge.") - logger.warn("Residual at level %d equals to %d" % (nlev, resid)) + logger.warn(f"Residual at level {nlev} equals to {resid}") # }}} End adaptively determine number of levels # {{{ adaptively determine brick quad order - resid = -1 if adaptive_quadrature: - table_tol = np.finfo(self.dtype).eps * 64 - logger.warn( - "Searching for n_brick_quad_points since " - "adaptive_quadrature=True") + table_tol = np.finfo(self.dtype).eps * 256 # 5e-14 for float64 + logger.warn("Searching for n_brick_quad_points since " + "adaptive_quadrature=True. Note that if you are using " + "special radial quadrature, the radial order will also be " + "adaptively refined.") - max_n_quad_pts = 200 + max_n_quad_pts = 1000 + resid = np.inf while True: - n_brick_quad_points = n_brick_quad_points + 3 + n_brick_quad_points += max(int(n_brick_quad_points * 0.2), 3) + if special_radial_brick_quadrature: + nradial_brick_quad_points += max( + int(nradial_brick_quad_points * 0.2), 3) + logger.warn( + f"Trying n_brick_quad_points = {n_brick_quad_points}, " + f"nradial_brick_quad_points = {nradial_brick_quad_points}, " + f"resid = {resid}") + else: + logger.warn( + f"Trying n_brick_quad_points = {n_brick_quad_points}, " + f"resid = {resid}") if n_brick_quad_points > max_n_quad_pts: logger.warn( "Adaptive quadrature refinement terminated " @@ -981,27 +980,31 @@ class NearFieldInteractionTable(object): max_n_quad_pts - 1)) break - data1 = drf( - queue, - source_box_extent=self.source_box_extent, - alpha=alpha, - nlevels=nlev, - n_brick_quad_points=n_brick_quad_points, - # extra_kernel_kwargs=extra_kernel_kwargs, - cheb_coefs=cheb_coefs, - **kwargs - ) - + drf = self.get_droste_table_builder( + n_brick_quad_points, + special_radial_brick_quadrature, nradial_brick_quad_points, + use_symmetry, knl_symmetry_tags) + data1 = drf(queue, source_box_extent=self.source_box_extent, + alpha=alpha, nlevels=nlev, + n_brick_quad_points=n_brick_quad_points, + # extra_kernel_kwargs=extra_kernel_kwargs, + cheb_coefs=cheb_coefs, **kwargs) + + resid_prev = resid resid = np.max(np.abs(data1 - data0)) / np.max(np.abs(data1)) data0 = data1 - if abs(resid) < table_tol: + if resid < table_tol: logger.warn( "Adaptive quadrature " "converged at order %d with residual %e" % ( n_brick_quad_points - 1, resid)) break + if resid > resid_prev: + logger.warn("Non-monotonic residual, breaking..") + break + if np.isnan(resid): logger.warn( "Adaptive quadrature terminated " @@ -1010,8 +1013,8 @@ class NearFieldInteractionTable(object): if resid >= table_tol: logger.warn("Adaptive quadrature failed to converge.") - logger.warn("Residual at order %d equals to %d" % ( - n_brick_quad_points, resid)) + logger.warn(f"Residual at order {n_brick_quad_points} " + f"equals to {resid}") if resid < 0: logger.warn("Failed to perform quadrature order refinement.") @@ -1219,13 +1222,9 @@ class NearFieldInteractionTable(object): # only for getting kernel evaluation related stuff drf = InverseDrosteReduced( - self.integral_knl, - self.quad_order, - self.interaction_case_vecs, - n_brick_quad_points=0, - knl_symmetry_tags=[], - auto_windowing=False - ) + self.integral_knl, self.quad_order, + self.interaction_case_vecs, n_brick_quad_points=0, + knl_symmetry_tags=[], auto_windowing=False) # uses "dist[dim]", assigned to "knl_val" knl_insns = drf.get_sumpy_kernel_insns()