diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 0d058c820944180ed27bea27579051cb71fc9b5b..4ecf5253e6c0f62f3d0717d7952932f146fc899f 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -172,7 +172,9 @@ class QBXLayerPotentialSource(LayerPotentialSource): fmm_level_to_order=( fmm_level_to_order or self.fmm_level_to_order), target_stick_out_factor=( - target_stick_out_factor or self.target_stick_out_factor), + target_stick_out_factor + if target_stick_out_factor is not None + else self.target_stick_out_factor), base_resampler=base_resampler or self._base_resampler, debug=( diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index 7fd7f49d4e957d216d0fd83c680740fb11e40521..426de45d3789a02fd392fe5fadc92ab286e72b33 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -369,7 +369,7 @@ def drive_fmm(expansion_wrangler, src_weights): # Interface guidelines: Attributes of the tree are assumed to be known # to the expansion wrangler and should not be passed. - logger.debug("start qbx fmm") + logger.info("start qbx fmm") logger.debug("reorder source weights") src_weights = wrangler.reorder_sources(src_weights) @@ -522,7 +522,7 @@ def drive_fmm(expansion_wrangler, src_weights): # }}} - logger.debug("qbx fmm complete") + logger.info("qbx fmm complete") return result diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index af43bb94307ebce188706ea072a11ee419f69592..236a80415743d95e12d7818bda2bbf3b2ac8aae6 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -34,7 +34,7 @@ from pyopencl.tools import ( # noqa from functools import partial from meshmode.mesh.generation import ( # noqa - ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, + ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, WobblyCircle, make_curve_mesh) from sumpy.visualization import FieldPlotter from pytential import bind, sym, norm @@ -702,6 +702,13 @@ def get_starfish_mesh(refinement_increment, target_order): target_order) +def get_wobbly_circle_mesh(refinement_increment, target_order): + nelements = [3000, 5000, 7000][refinement_increment] + return make_curve_mesh(WobblyCircle.random(30, seed=30), + np.linspace(0, 1, nelements+1), + target_order) + + def get_sphere_mesh(refinement_increment, target_order): from meshmode.mesh.generation import generate_icosphere mesh = generate_icosphere(1, target_order) @@ -906,8 +913,12 @@ def test_off_surface_eval(ctx_getter, use_fmm, do_plot=False): pre_density_discr = Discretization( cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) - qbx, _ = QBXLayerPotentialSource(pre_density_discr, 4*target_order, qbx_order, - fmm_order=fmm_order).with_refinement() + qbx, _ = QBXLayerPotentialSource( + pre_density_discr, + 4*target_order, + qbx_order, + fmm_order=fmm_order, + ).with_refinement() density_discr = qbx.density_discr @@ -922,10 +933,11 @@ def test_off_surface_eval(ctx_getter, use_fmm, do_plot=False): (qbx, PointsTarget(fplot.points)), op)(queue, sigma=sigma) - print(fld_in_vol) - err = cl.clmath.fabs(fld_in_vol - (-1)) + linf_err = cl.array.max(err).get() + print("l_inf error:", linf_err) + if do_plot: fplot.show_scalar_in_matplotlib(fld_in_vol.get()) import matplotlib.pyplot as pt @@ -933,7 +945,86 @@ def test_off_surface_eval(ctx_getter, use_fmm, do_plot=False): pt.show() # FIXME: Why does the FMM only meet this sloppy tolerance? - assert (err < 1e-2).get().all() + assert linf_err < 1e-2 + +# }}} + + +# {{{ test off-surface eval vs direct + +def test_off_surface_eval_vs_direct(ctx_getter, do_plot=False): + logging.basicConfig(level=logging.INFO) + + cl_ctx = ctx_getter() + queue = cl.CommandQueue(cl_ctx) + + # prevent cache 'splosion + from sympy.core.cache import clear_cache + clear_cache() + + nelements = 300 + target_order = 8 + qbx_order = 3 + + mesh = make_curve_mesh(WobblyCircle.random(8, seed=30), + np.linspace(0, 1, nelements+1), + target_order) + + from pytential.qbx import QBXLayerPotentialSource + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + pre_density_discr = Discretization( + cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) + direct_qbx, _ = QBXLayerPotentialSource( + pre_density_discr, 4*target_order, qbx_order, + fmm_order=False, + target_stick_out_factor=0.05, + ).with_refinement() + fmm_qbx, _ = QBXLayerPotentialSource( + pre_density_discr, 4*target_order, qbx_order, + fmm_order=qbx_order + 3, + expansion_disks_in_tree_have_extent=True, + target_stick_out_factor=0.05, + ).with_refinement() + + fplot = FieldPlotter(np.zeros(2), extent=5, npoints=1000) + from pytential.target import PointsTarget + ptarget = PointsTarget(fplot.points) + from sumpy.kernel import LaplaceKernel + + op = sym.D(LaplaceKernel(2), sym.var("sigma"), qbx_forced_limit=None) + + from pytential.qbx import QBXTargetAssociationFailedException + try: + direct_density_discr = direct_qbx.density_discr + direct_sigma = direct_density_discr.zeros(queue) + 1 + direct_fld_in_vol = bind((direct_qbx, ptarget), op)( + queue, sigma=direct_sigma) + + except QBXTargetAssociationFailedException as e: + fplot.show_scalar_in_matplotlib(e.failed_target_flags.get(queue)) + import matplotlib.pyplot as pt + pt.show() + raise + + fmm_density_discr = fmm_qbx.density_discr + fmm_sigma = fmm_density_discr.zeros(queue) + 1 + fmm_fld_in_vol = bind((fmm_qbx, ptarget), op)(queue, sigma=fmm_sigma) + + err = cl.clmath.fabs(fmm_fld_in_vol - direct_fld_in_vol) + + linf_err = cl.array.max(err).get() + print("l_inf error:", linf_err) + + if do_plot: + #fplot.show_scalar_in_mayavi(0.1*cl.clmath.log10(1e-15 + err).get(queue)) + fplot.show_scalar_in_mayavi(fmm_fld_in_vol.get(queue)) + import mayavi.mlab as mlab + mlab.show() + + assert linf_err < 1e-3 # }}}