From 4f2db07b8290b1195a6d52d5b29af722c0cbdc2e Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 17 Aug 2020 20:57:43 -0500 Subject: [PATCH 01/11] Refactor droste --- volumential/droste.py | 185 +++++++++++++++++++++++++++--------------- 1 file changed, 120 insertions(+), 65 deletions(-) diff --git a/volumential/droste.py b/volumential/droste.py index 45845df..77773b1 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,13 @@ 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 + if self.dim is not None: self.n_q_points = quad_order ** self.dim else: @@ -157,6 +171,27 @@ 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 + pass + 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 +292,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 +311,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 +345,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 +623,7 @@ class DrosteBase(KernelCacheWrapper): * knl_scaling ) {id=result,dep=jac:mpoint:density:finish_kval} end - """ - ) + """) ] def get_target_points(self, queue=None): @@ -630,9 +686,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 +701,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[0], 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 +732,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 +776,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: @@ -762,36 +830,21 @@ class DrosteFull(DrosteBase): 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 + brick_quadrature_kwargs = self.make_brick_quadrature_kwargs() 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, + **brick_quadrature_kwargs, + **extra_kernel_kwargs) cheb_table = res["result"] @@ -804,6 +857,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): @@ -843,10 +897,13 @@ class DrosteReduced(DrosteBase): quad_order=None, case_vecs=None, n_brick_quad_points=50, + special_radial_quadrature=False, + nradial_quad_points=None, knl_symmetry_tags=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 @@ -1443,6 +1500,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 +1569,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, + special_radial_quadrature=False, nradial_quad_points=None, + knl_symmetry_tags=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 +1589,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 -- GitLab From 4b4568ae16a4c58493d22dc271700336db190ea3 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 17 Aug 2020 21:14:25 -0500 Subject: [PATCH 02/11] Fixes --- volumential/droste.py | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/volumential/droste.py b/volumential/droste.py index 77773b1..93b98de 100644 --- a/volumential/droste.py +++ b/volumential/droste.py @@ -187,7 +187,7 @@ class DrosteBase(KernelCacheWrapper): # 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 - pass + raise NotImplementedError else: return dict(quadrature_nodes=legendre_nodes, quadrature_weights=legendre_weights) @@ -897,9 +897,9 @@ class DrosteReduced(DrosteBase): quad_order=None, case_vecs=None, n_brick_quad_points=50, + knl_symmetry_tags=None, special_radial_quadrature=False, nradial_quad_points=None, - knl_symmetry_tags=None, ): super(DrosteReduced, self).__init__( integral_knl, quad_order, case_vecs, n_brick_quad_points, @@ -1333,22 +1333,17 @@ class DrosteReduced(DrosteBase): except Exception: pass knl = self.get_cached_optimized_kernel() + brick_quadrature_kwargs = self.make_brick_quadrature_kwargs() 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, **brick_quadrature_kwargs) + raw_cheb_table_case = res["result"] self.get_kernel_id = 1 @@ -1570,9 +1565,9 @@ class InverseDrosteReduced(DrosteReduced): """ def __init__(self, integral_knl, quad_order, case_vecs, - n_brick_quad_points=50, + n_brick_quad_points=50, knl_symmetry_tags=None, special_radial_quadrature=False, nradial_quad_points=None, - knl_symmetry_tags=None, auto_windowing=True): + auto_windowing=True): """ :param auto_windowing: auto-detect window radius. """ -- GitLab From 85f77cc3ba2eda4524b79a3326aab0ed085b2dad Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 18 Aug 2020 23:30:53 -0500 Subject: [PATCH 03/11] Set pycodestyle line length --- setup.cfg | 3 +++ 1 file changed, 3 insertions(+) diff --git a/setup.cfg b/setup.cfg index 77f81d3..65a7a4a 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 -- GitLab From 34526526a31fe3e57f640891f19113ec727fc5db Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 20 Aug 2020 10:24:55 -0500 Subject: [PATCH 04/11] Refactor brick quadrature --- setup.py | 1 + volumential/droste.py | 29 ++++------------------------- 2 files changed, 5 insertions(+), 25 deletions(-) diff --git a/setup.py b/setup.py index 5004b8d..8b037de 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 93b98de..a0f39b4 100644 --- a/volumential/droste.py +++ b/volumential/droste.py @@ -830,8 +830,8 @@ class DrosteFull(DrosteBase): t = np.array([pt[-1] for pt in q_points[: self.ntgt_points]]) # quad formula for each brick (normalized to [0,1]) - brick_quadrature_kwargs = self.make_brick_quadrature_kwargs() + knl = self.get_cached_optimized_kernel() evt, res = knl( queue, alpha=alpha, result=result_array, @@ -1317,16 +1317,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") @@ -2104,14 +2094,7 @@ 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 + brick_quadrature_kwargs = self.make_brick_quadrature_kwargs() # --------- call kernel 0 ---------- self.get_kernel_id = 0 @@ -2127,15 +2110,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, **brick_quadrature_kwargs ) # --------- call kernel 1 ---------- @@ -2152,15 +2133,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, **brick_quadrature_kwargs ) # --------- call kernel 2 ---------- -- GitLab From dfa48850dda39fb34c15c89bcebf921da2415a3b Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 12 Oct 2020 10:58:13 -0500 Subject: [PATCH 05/11] Use correct scaling for tanh-sinh --- volumential/droste.py | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/volumential/droste.py b/volumential/droste.py index a0f39b4..1018090 100644 --- a/volumential/droste.py +++ b/volumential/droste.py @@ -187,7 +187,33 @@ class DrosteBase(KernelCacheWrapper): # 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 - raise NotImplementedError + 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(100): + nodes = ts.calc_nodes(degree=deg, prec=prec) + if len(nodes) >= self.nradial_quad_points: + 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) + + # 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) -- GitLab From 3d7ca0376d131d8c10acf218acc0f93968f57c95 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 12 Oct 2020 11:14:20 -0500 Subject: [PATCH 06/11] Fix scaling for tanh-sinh when deg=1 --- volumential/droste.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/volumential/droste.py b/volumential/droste.py index 1018090..8061fc9 100644 --- a/volumential/droste.py +++ b/volumential/droste.py @@ -194,7 +194,7 @@ class DrosteBase(KernelCacheWrapper): 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(100): + for deg in range(1, 100): nodes = ts.calc_nodes(degree=deg, prec=prec) if len(nodes) >= self.nradial_quad_points: break @@ -204,6 +204,8 @@ class DrosteBase(KernelCacheWrapper): 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 -- GitLab From 92d400c4a3c0154fd8d59a873ef54a3b68b9fb00 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 13 Oct 2020 20:16:38 -0500 Subject: [PATCH 07/11] Allow tanh-sinh to change quadrature size --- volumential/droste.py | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/volumential/droste.py b/volumential/droste.py index 8061fc9..fece03e 100644 --- a/volumential/droste.py +++ b/volumential/droste.py @@ -114,6 +114,7 @@ class DrosteBase(KernelCacheWrapper): 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 @@ -197,6 +198,7 @@ class DrosteBase(KernelCacheWrapper): 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 @@ -732,7 +734,7 @@ class DrosteFull(DrosteBase): if self.special_radial_quadrature: quad_vars_subdomain = self.make_brick_domain( quad_vars[1:], self.nquad_points) & self.make_brick_domain( - quad_vars[0], self.nradial_quad_points) + quad_vars[:1], self.nradial_quad_points) else: quad_vars_subdomain = self.make_brick_domain( quad_vars, self.nquad_points) @@ -857,9 +859,6 @@ 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]) - brick_quadrature_kwargs = self.make_brick_quadrature_kwargs() - knl = self.get_cached_optimized_kernel() evt, res = knl( queue, alpha=alpha, result=result_array, @@ -871,7 +870,7 @@ class DrosteFull(DrosteBase): np.float64, copy=True), n_cases=self.ncases, nfunctions=self.nfunctions, quad_order=self.ntgt_points, nlevels=nlevels, - **brick_quadrature_kwargs, + **self.brick_quadrature_kwargs, **extra_kernel_kwargs) cheb_table = res["result"] @@ -967,8 +966,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) @@ -1351,7 +1359,6 @@ class DrosteReduced(DrosteBase): except Exception: pass knl = self.get_cached_optimized_kernel() - brick_quadrature_kwargs = self.make_brick_quadrature_kwargs() evt, res = knl( queue, alpha=alpha, result=result_array, root_brick=root_brick, @@ -1360,7 +1367,7 @@ class DrosteReduced(DrosteBase): 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, **brick_quadrature_kwargs) + **extra_kernel_kwargs, **self.brick_quadrature_kwargs) raw_cheb_table_case = res["result"] @@ -2121,9 +2128,6 @@ class InverseDrosteReduced(DrosteReduced): ] ) - # quad formula for each brick (normalized to [0,1]) - brick_quadrature_kwargs = self.make_brick_quadrature_kwargs() - # --------- call kernel 0 ---------- self.get_kernel_id = 0 try: -- GitLab From aac867f718a378edb2e2fad285c2220c4ccbc4e8 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 13 Oct 2020 21:09:00 -0500 Subject: [PATCH 08/11] Tweak table adaptive build logic --- volumential/nearfield_potential_table.py | 171 +++++++++++------------ 1 file changed, 78 insertions(+), 93 deletions(-) diff --git a/volumential/nearfield_potential_table.py b/volumential/nearfield_potential_table.py index cb0bcbd..bcc863b 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,14 @@ 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,60 +877,24 @@ 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 @@ -926,15 +917,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 @@ -963,15 +949,19 @@ class NearFieldInteractionTable(object): 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") + 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 not be " + "adaptively refined.") max_n_quad_pts = 200 + resid = np.inf while True: n_brick_quad_points = n_brick_quad_points + 3 + 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,16 +971,15 @@ 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 = np.max(np.abs(data1 - data0)) / np.max(np.abs(data1)) data0 = data1 @@ -1219,13 +1208,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() -- GitLab From ab20cc9233cc20f249938e0eb48149c964389835 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 13 Oct 2020 21:24:26 -0500 Subject: [PATCH 09/11] More agressive adaptation of brick orders --- volumential/nearfield_potential_table.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/volumential/nearfield_potential_table.py b/volumential/nearfield_potential_table.py index bcc863b..a40bc46 100644 --- a/volumential/nearfield_potential_table.py +++ b/volumential/nearfield_potential_table.py @@ -940,7 +940,7 @@ 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("Residual at level %d equals to %f" % (nlev, resid)) # }}} End adaptively determine number of levels @@ -951,17 +951,24 @@ class NearFieldInteractionTable(object): table_tol = np.finfo(self.dtype).eps * 64 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 not be " + "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 - logger.warn(f"Trying n_brick_quad_points = {n_brick_quad_points}, " - f"resid = {resid}") + 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 " -- GitLab From 143d29eca0b25fbfa957db5874065f3daeab90d3 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 13 Oct 2020 21:40:58 -0500 Subject: [PATCH 10/11] Minor fixes --- volumential/droste.py | 7 +++++-- volumential/nearfield_potential_table.py | 15 +++++++++------ 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/volumential/droste.py b/volumential/droste.py index fece03e..382fa2f 100644 --- a/volumential/droste.py +++ b/volumential/droste.py @@ -831,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 @@ -2148,7 +2151,7 @@ class InverseDrosteReduced(DrosteReduced): nfunctions=self.nfunctions, quad_order=self.ntgt_points, nlevels=nlevels, - **extra_kernel_kwargs, **brick_quadrature_kwargs + **extra_kernel_kwargs, **self.brick_quadrature_kwargs ) # --------- call kernel 1 ---------- @@ -2171,7 +2174,7 @@ class InverseDrosteReduced(DrosteReduced): nfunctions=self.nfunctions, quad_order=self.ntgt_points, nlevels=nlevels, - **extra_kernel_kwargs, **brick_quadrature_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 a40bc46..2d5fda4 100644 --- a/volumential/nearfield_potential_table.py +++ b/volumential/nearfield_potential_table.py @@ -858,7 +858,8 @@ class NearFieldInteractionTable(object): nlev = 1 if "special_radial_brick_quadrature" in kwargs: - special_radial_brick_quadrature = kwargs.pop("special_radial_brick_quadrature") + 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 @@ -963,12 +964,14 @@ class NearFieldInteractionTable(object): 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}") + 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}") + 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 " -- GitLab From 22d9b24333237fd73e275b584d0598c98eb6918a Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 13 Oct 2020 21:49:01 -0500 Subject: [PATCH 11/11] Optimize convergence criteria --- volumential/nearfield_potential_table.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/volumential/nearfield_potential_table.py b/volumential/nearfield_potential_table.py index 2d5fda4..8cef4d2 100644 --- a/volumential/nearfield_potential_table.py +++ b/volumential/nearfield_potential_table.py @@ -902,7 +902,7 @@ class NearFieldInteractionTable(object): 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: @@ -941,15 +941,14 @@ class NearFieldInteractionTable(object): if resid >= table_tol: logger.warn("Adaptive level refinement failed to converge.") - logger.warn("Residual at level %d equals to %f" % (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 + 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 " @@ -991,16 +990,21 @@ class NearFieldInteractionTable(object): # 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 " @@ -1009,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.") -- GitLab