From 9d620e2f74f44e790507a444ecadab7a204963f4 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sat, 7 Oct 2017 13:45:30 -0500
Subject: [PATCH] Fiddle with Helmholtz error model in
SimpleExpansionOrderFinder
---
sumpy/expansion/level_to_order.py | 19 +++++++++++++------
1 file changed, 13 insertions(+), 6 deletions(-)
diff --git a/sumpy/expansion/level_to_order.py b/sumpy/expansion/level_to_order.py
index 727f2d80..75a5a6c2 100644
--- a/sumpy/expansion/level_to_order.py
+++ b/sumpy/expansion/level_to_order.py
@@ -102,13 +102,15 @@ class SimpleExpansionOrderFinder(object):
.. math::
- C_{\text{helm}} \frac 1{p!} (hk)^{p+1},
+ C_{\text{helm}} \frac 1{p!}
+ \left(C_{\text{helm\_scale}} \cdot \frac{hk}{2\pi}\right)^{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, err_const_laplace=0.01, err_const_helmholtz=100,
+ scaling_const_helmholtz=4,
extra_order=1):
"""
:arg extra_order: order increase to accommodate, say, the taking of
@@ -118,6 +120,7 @@ class SimpleExpansionOrderFinder(object):
self.err_const_laplace = err_const_laplace
self.err_const_helmholtz = err_const_helmholtz
+ self.scaling_const_helmholtz = scaling_const_helmholtz
self.extra_order = extra_order
@@ -140,13 +143,17 @@ class SimpleExpansionOrderFinder(object):
tree.stick_out_factor
* tree.root_extent / (1 << level))
- factor = box_lengthscale * helmholtz_k / (2*float(np.pi))
+ factor = (
+ self.scaling_const_helmholtz
+ * box_lengthscale
+ * helmholtz_k
+ / (2*float(np.pi)))
from math import factorial
helm_order = 1
- helm_error = self.err_const_helmholtz * factor
+ helm_rec_error = self.err_const_helmholtz * factor
while True:
- helm_error = helm_error * factor / (helm_order+1)
+ helm_rec_error = helm_rec_error * factor / (helm_order+1)
if helm_order < 4:
# this may overflow for large orders
@@ -154,10 +161,10 @@ class SimpleExpansionOrderFinder(object):
1/factorial(helm_order+1)
* self.err_const_helmholtz
* factor**(helm_order+1))
- assert (abs(helm_error - helm_error_direct)
+ assert (abs(helm_rec_error - helm_error_direct)
< 1e-13 * abs(helm_error_direct))
- if helm_error < self.tol:
+ if helm_rec_error * helm_order**(tree.dimensions-1) < self.tol:
break
helm_order += 1
--
GitLab