From cf0932d6807fee6d3ae1c0d3fae1a8ce6e9028b4 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Mon, 25 Sep 2017 17:22:50 -0500
Subject: [PATCH 1/3] FMM-level-to-order: also pass tree

---
 pytential/qbx/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py
index 60f439ef..a5a5d115 100644
--- a/pytential/qbx/__init__.py
+++ b/pytential/qbx/__init__.py
@@ -132,7 +132,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
             if fmm_order is False:
                 fmm_level_to_order = False
             else:
-                def fmm_level_to_order(level):
+                def fmm_level_to_order(tree, level):
                     return fmm_order
 
         # }}}
-- 
GitLab


From 8365c7d05720ebb5687962af7243a1ed0d823df0 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Mon, 25 Sep 2017 17:23:25 -0500
Subject: [PATCH 2/3] Add Helmholtz order finder to Maxwell test (for now)

---
 test/test_maxwell.py | 77 +++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 69 insertions(+), 8 deletions(-)

diff --git a/test/test_maxwell.py b/test/test_maxwell.py
index a27e066e..8eaabf0c 100644
--- a/test/test_maxwell.py
+++ b/test/test_maxwell.py
@@ -50,12 +50,12 @@ class MaxwellTestCase:
     fmm_backend = "fmmlib"
 
     def __init__(self, k, is_interior, resolutions, qbx_order,
-            fmm_order):
+            fmm_tolerance):
         self.k = k
         self.is_interior = is_interior
         self.resolutions = resolutions
         self.qbx_order = qbx_order
-        self.fmm_order = fmm_order
+        self.fmm_tolerance = fmm_tolerance
 
 
 class SphereTestCase(MaxwellTestCase):
@@ -186,15 +186,15 @@ class ElliptiPlaneTestCase(MaxwellTestCase):
 
 
 tc_int = SphereTestCase(k=1.2, is_interior=True, resolutions=[0, 1],
-        qbx_order=3, fmm_order=5)
+        qbx_order=3, fmm_tolerance=1e-4)
 tc_ext = SphereTestCase(k=1.2, is_interior=False, resolutions=[0, 1],
-        qbx_order=3, fmm_order=5)
+        qbx_order=3, fmm_tolerance=1e-4)
 
 tc_rc_ext = RoundedCubeTestCase(k=6.4, is_interior=False, resolutions=[0.1],
-        qbx_order=3, fmm_order=10)
+        qbx_order=3, fmm_tolerance=1e-4)
 
-tc_plane_ext = ElliptiPlaneTestCase(k=1.4, is_interior=False, resolutions=[0.2],
-        qbx_order=3, fmm_order=10)
+tc_plane_ext = ElliptiPlaneTestCase(k=0.4, is_interior=False, resolutions=[0.2],
+        qbx_order=3, fmm_tolerance=1e-3)
 
 
 class EHField(object):
@@ -211,6 +211,65 @@ class EHField(object):
         return self.field[3:]
 
 
+class HelmholtzFMMOrderFinder(object):
+    r"""
+    This models the error as:
+
+    .. math::
+
+        C_{\text{lap}} \left(\frac{\sqrt{d}}{3}\right)^{p+1}
+        +
+        C_{\text{helm}} \frac 1{p!} (hk)^{p+1},
+
+    where :math:`d` is the number of dimensions, :math:`p` is the expansion order,
+    :math:`h` is the box size, and :math:`k` is the wave number.
+    """
+
+    def __init__(self, tol, helmholtz_k, err_const_laplace=0.01,
+            err_const_helmholtz=100, extra_order=1):
+        """
+        :arg extra_order: order increase to accommodate, say, the taking of
+            oderivatives f the FMM expansions.
+        """
+        self.tol = tol
+        self.helmholtz_k = helmholtz_k
+
+        self.err_const_laplace = err_const_laplace
+        self.err_const_helmholtz = err_const_helmholtz
+
+        self.extra_order = extra_order
+
+    def __call__(self, tree, level):
+        laplace_order = int(np.ceil(
+                (np.log(self.tol) - np.log(self.err_const_laplace))
+                /
+                np.log(
+                    np.sqrt(tree.dimensions)/3
+                    ) - 1))
+
+        box_lengthscale = (
+            tree.stick_out_factor
+            * tree.root_extent / (1 << level))
+
+        from math import factorial
+        helm_order = 1
+        while True:
+            helm_error = (
+                    1/factorial(helm_order+1)
+                    * self.err_const_helmholtz
+                    * (box_lengthscale * self.helmholtz_k)**(helm_order+1))
+            if helm_error < self.tol:
+                break
+
+            helm_order += 1
+
+            if helm_order > 500:
+                raise ValueError("unable to find suitable order estimate "
+                        "for Helmholtz expansion")
+
+        return max(laplace_order, helm_order) + self.extra_order
+
+
 # {{{ driver
 
 @pytest.mark.parametrize("case", [
@@ -306,7 +365,9 @@ def test_pec_mfie_extinction(ctx_getter, case, visualize=False):
         qbx, _ = QBXLayerPotentialSource(
                 pre_scat_discr, fine_order=4*case.target_order,
                 qbx_order=case.qbx_order,
-                fmm_order=case.fmm_order, fmm_backend=case.fmm_backend
+                fmm_level_to_order=HelmholtzFMMOrderFinder(
+                    case.fmm_tolerance, case.k),
+                fmm_backend=case.fmm_backend
                 ).with_refinement(_expansion_disturbance_tolerance=0.05)
         h_max = qbx.h_max
 
-- 
GitLab


From b39240e2b34afb144ccc095c5794d2fb6ded351b Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Mon, 25 Sep 2017 23:59:55 -0500
Subject: [PATCH 3/3] Fix another occurrence of fmm_level_to_order for being
 passed a tree

---
 test/test_layer_pot.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py
index af642093..e58fba64 100644
--- a/test/test_layer_pot.py
+++ b/test/test_layer_pot.py
@@ -292,7 +292,7 @@ def test_unregularized_off_surface_fmm_vs_direct(ctx_getter):
             density_discr,
             fmm_order=False,
             )
-    fmm = direct.copy(fmm_level_to_order=lambda _: fmm_order)
+    fmm = direct.copy(fmm_level_to_order=lambda tree, level: fmm_order)
 
     sigma = density_discr.zeros(queue) + 1
 
-- 
GitLab