inteq_tests.py

#!/usr/bin/env python3
"""Integral equation test infrastructure"""

from __future__ import division, absolute_import, print_function

__copyright__ = "Copyright (C) 2014 Andreas Kloeckner"

__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""


import numpy as np
import numpy.linalg as la
import pyopencl as cl
import pyopencl.clmath  # noqa
import gzip
import os
import pickle
import pytest
from pytools import Record

from functools import partial
from meshmode.mesh.generation import (  # noqa
        ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, WobblyCircle,
        make_curve_mesh)
from meshmode.discretization.visualization import make_visualizer
from pytential import bind, sym
from pytential.qbx import QBXTargetAssociationFailedException
from pyopencl.tools import pytest_generate_tests_for_pyopencl as pytest_generate_tests  # noqa


import logging
logger = logging.getLogger(__name__)

circle = partial(ellipse, 1)

try:
    import matplotlib.pyplot as pt
except ImportError:
    pass


def make_circular_point_group(
        ambient_dim, npoints, radius,
        center=np.array([0., 0.]), func=lambda x: x):
    t = func(np.linspace(0, 1, npoints, endpoint=False)) * (2 * np.pi)
    center = np.asarray(center)
    result = np.zeros((ambient_dim, npoints))
    result[:2, :] = (
            center[:, np.newaxis]
            + radius*np.vstack((np.cos(t), np.sin(t))))
    return result


# {{{ test cases

class IntEqTestCase:
    def __init__(self, helmholtz_k, bc_type, prob_side):
        """
        :arg prob_side: may be -1, +1, or ``'scat'`` for a scattering problem
        """

        if helmholtz_k is None:
            helmholtz_k = self.default_helmholtz_k

        self.helmholtz_k = helmholtz_k
        self.bc_type = bc_type
        self.prob_side = prob_side

    @property
    def k(self):
        return self.helmholtz_k

    def __str__(self):
        return (
                "name: %s, bc_type: %s, prob_side: %s, "
                "helmholtz_k: %s, qbx_order: %d, target_order: %d"
                % (
                    self.name, self.bc_type, self.prob_side, self.helmholtz_k,
                    self.qbx_order, self.target_order))

    fmm_backend = "sumpy"
    gmres_tol = 1e-14


class SphereIntEqTestCase(IntEqTestCase):
    resolutions = [1, 2]
    name = "sphere"

    def get_mesh(self, resolution, target_order):
        from meshmode.mesh.generation import generate_icosphere
        from meshmode.mesh.refinement import refine_uniformly
        mesh = refine_uniformly(
                generate_icosphere(1, target_order),
                resolution)

        return mesh

    fmm_backend = "fmmlib"
    use_refinement = True

    fmm_tol = 1e-4

    inner_radius = 0.4
    outer_radius = 5

    qbx_order = 3
    target_order = 8
    check_gradient = False
    check_tangential_deriv = False

    gmres_tol = 1e-3


class BetterplaneIntEqTestCase(IntEqTestCase):
    name = "betterplane"

    default_helmholtz_k = 20
    resolutions = [0.25]

    fmm_backend = "fmmlib"
    use_refinement = True

    qbx_order = 4
    fmm_tol = 1e-5
    target_order = 7
    check_gradient = False
    check_tangential_deriv = False

    visualize_geometry = False

    scaled_max_curvature_threshold = 5
    expansion_disturbance_tolerance = 0.1
    refinement_maxiter = 20

    gmres_tol = 1e-6

    vis_grid_spacing = (0.025, 0.2, 0.025)
    vis_extend_factor = 0.2

    def __init__(self, helmholtz_k, bc_type, prob_side):
        IntEqTestCase.__init__(self, helmholtz_k, bc_type, prob_side)

    def get_mesh(self, resolution, target_order):
        # mesh order is fixed at 2
        del target_order

        assert resolution == 0.25

        try:
            with gzip.open(f"plane-mesh.pkl.gz", "rb") as inf:
                return pickle.load(inf)
        except FileNotFoundError:
            pass

        from pytools import download_from_web_if_not_present

        download_from_web_if_not_present(
                "https://raw.githubusercontent.com/inducer/geometries/a869fc3/"
                "surface-3d/betterplane.brep")

        from meshmode.mesh.io import generate_gmsh, ScriptWithFilesSource
        mesh = generate_gmsh(
                ScriptWithFilesSource("""
                    Merge "betterplane.brep";

                    Mesh.CharacteristicLengthMax = %(lcmax)f;
                    Mesh.ElementOrder = 2;
                    Mesh.CharacteristicLengthExtendFromBoundary = 0;

                    // 2D mesh optimization
                    // Mesh.Lloyd = 1;

                    l_superfine() = Unique(Abs(Boundary{ Surface{
                        27, 25, 17, 13, 18  }; }));
                    l_fine() = Unique(Abs(Boundary{ Surface{ 2, 6, 7}; }));
                    l_coarse() = Unique(Abs(Boundary{ Surface{ 14, 16  }; }));

                    // p() = Unique(Abs(Boundary{ Line{l_fine()}; }));
                    // Characteristic Length{p()} = 0.05;

                    Field[1] = Attractor;
                    Field[1].NNodesByEdge = 100;
                    Field[1].EdgesList = {l_superfine()};

                    Field[2] = Threshold;
                    Field[2].IField = 1;
                    Field[2].LcMin = 0.075;
                    Field[2].LcMax = %(lcmax)f;
                    Field[2].DistMin = 0.1;
                    Field[2].DistMax = 0.4;

                    Field[3] = Attractor;
                    Field[3].NNodesByEdge = 100;
                    Field[3].EdgesList = {l_fine()};

                    Field[4] = Threshold;
                    Field[4].IField = 3;
                    Field[4].LcMin = 0.1;
                    Field[4].LcMax = %(lcmax)f;
                    Field[4].DistMin = 0.15;
                    Field[4].DistMax = 0.4;

                    Field[5] = Attractor;
                    Field[5].NNodesByEdge = 100;
                    Field[5].EdgesList = {l_coarse()};

                    Field[6] = Threshold;
                    Field[6].IField = 5;
                    Field[6].LcMin = 0.2;
                    Field[6].LcMax = %(lcmax)f;
                    Field[6].DistMin = 0.2;
                    Field[6].DistMax = 0.4;

                    Field[7] = Min;
                    Field[7].FieldsList = {2, 4, 6};

                    Background Field = 7;
                    """ % {
                        "lcmax": resolution,
                        }, ["betterplane.brep"]), 2)

        # Flip elements--gmsh generates inside-out geometry.
        from meshmode.mesh.processing import perform_flips
        mesh = perform_flips(mesh, np.ones(mesh.nelements))

        with gzip.open("plane-mesh.pkl.gz", "wb") as outf:
            pickle.dump(mesh, outf)

        return mesh

    inner_radius = 0.2
    outer_radius = 15

# }}}


# {{{ lpot source from case

def lpot_source_from_case(cl_ctx, queue, case, resolution, base_lpot_kwargs):
    mesh = case.get_mesh(resolution, case.target_order)
    logger.info("%d elements" % mesh.nelements)

    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(case.target_order))

    source_order = 4*case.target_order

    refiner_extra_kwargs = {}

    qbx_lpot_kwargs = base_lpot_kwargs.copy()

    if case.fmm_backend is None:
        qbx_lpot_kwargs["fmm_order"] = False
    else:
        if hasattr(case, "fmm_tol"):
            from sumpy.expansion.level_to_order import\
                    SimpleExpansionOrderFinder
            qbx_lpot_kwargs["fmm_level_to_order"] = SimpleExpansionOrderFinder(
                    case.fmm_tol)

        elif hasattr(case, "fmm_order"):
            qbx_lpot_kwargs["fmm_order"] = case.fmm_order
        else:
            qbx_lpot_kwargs["fmm_order"] = case.qbx_order + 5

    qbx = QBXLayerPotentialSource(
            pre_density_discr,
            fine_order=source_order,
            qbx_order=case.qbx_order,
            fmm_backend="fmmlib",
            **qbx_lpot_kwargs)

    if case.use_refinement:
        if case.k != 0 and getattr(case, "refine_on_helmholtz_k", True):
            refiner_extra_kwargs["kernel_length_scale"] = 5/case.k

        if hasattr(case, "scaled_max_curvature_threshold"):
            refiner_extra_kwargs["_scaled_max_curvature_threshold"] = \
                    case.scaled_max_curvature_threshold

        if hasattr(case, "expansion_disturbance_tolerance"):
            refiner_extra_kwargs["_expansion_disturbance_tolerance"] = \
                    case.expansion_disturbance_tolerance

        if hasattr(case, "refinement_maxiter"):
            refiner_extra_kwargs["maxiter"] = case.refinement_maxiter

        logger.info(
                "%d elements before refinement",
                pre_density_discr.mesh.nelements)
        qbx, _ = qbx.with_refinement(**refiner_extra_kwargs)

        logger.info(
                "%d stage-1 elements after refinement",
                qbx.density_discr.mesh.nelements)
        logger.info(
                "%d stage-2 elements after refinement",
                qbx.stage2_density_discr.mesh.nelements)
        logger.info(
                "quad stage-2 elements have %d nodes",
                qbx.quad_stage2_density_discr.groups[0].nunit_nodes)

    return qbx

# }}}


# {{{ test backend

def run_int_eq_test(
        cl_ctx, queue, case, resolution, visualize, lpot_kwargs=None,
        output_dir=None):
    if lpot_kwargs is None:
        lpot_kwargs = {}
    qbx = lpot_source_from_case(cl_ctx, queue, case, resolution, lpot_kwargs)

    def output_file_path(fname):
        if output_dir is not None:
            os.makedirs(output_dir, exist_ok=True)
            return os.path.join(output_dir, fname)
        return fname

    density_discr = qbx.density_discr
    mesh = density_discr.mesh

    if (
            visualize
            and getattr(case, "visualize_geometry", False)):
        bdry_normals = bind(
                density_discr, sym.normal(mesh.ambient_dim)
                )(queue).as_vector(dtype=object)

        bdry_vis = make_visualizer(queue, density_discr, case.target_order)
        bdry_vis.write_vtk_file(output_file_path("geometry.vtu"), [
                ("normals", bdry_normals)
                ])

    # {{{ plot geometry

    if 0:
        if mesh.ambient_dim == 2:
            # show geometry, centers, normals
            nodes_h = density_discr.nodes().get(queue=queue)
            pt.plot(nodes_h[0], nodes_h[1], "x-")
            normal = (
                    bind(density_discr, sym.normal(2))(queue)
                    ).as_vector(np.object)
            pt.quiver(
                    nodes_h[0], nodes_h[1],
                    normal[0].get(queue), normal[1].get(queue))
            pt.gca().set_aspect("equal")
            pt.show()

        elif mesh.ambient_dim == 3:
            bdry_vis = make_visualizer(
                    queue, density_discr, case.target_order+3)

            bdry_normals = (
                    bind(density_discr, sym.normal(3))(queue)
                    ).as_vector(dtype=object)

            bdry_vis.write_vtk_file(
                    output_file_path("pre-solve-source-%s.vtu" % resolution),
                    [("bdry_normals", bdry_normals),])

        else:
            raise ValueError("invalid mesh dim")

    # }}}

    # {{{ set up operator

    from pytential.symbolic.pde.scalar import (
            DirichletOperator,
            NeumannOperator)

    from sumpy.kernel import LaplaceKernel, HelmholtzKernel
    if case.k:
        knl = HelmholtzKernel(mesh.ambient_dim)
        knl_kwargs = {"k": sym.var("k")}
        concrete_knl_kwargs = {"k": case.k}
    else:
        knl = LaplaceKernel(mesh.ambient_dim)
        knl_kwargs = {}
        concrete_knl_kwargs = {}

    if knl.is_complex_valued:
        dtype = np.complex128
    else:
        dtype = np.float64

    loc_sign = +1 if case.prob_side in [+1, "scat"] else -1

    if case.bc_type == "dirichlet":
        op = DirichletOperator(
                knl, loc_sign, use_l2_weighting=True,
                kernel_arguments=knl_kwargs)
    elif case.bc_type == "neumann":
        op = NeumannOperator(
                knl, loc_sign, use_l2_weighting=True,
                use_improved_operator=False, kernel_arguments=knl_kwargs)
    else:
        assert False

    op_u = op.operator(sym.var("u"))

    # }}}

    # {{{ set up test data

    if case.prob_side == -1:
        test_src_geo_radius = case.outer_radius
        test_tgt_geo_radius = case.inner_radius
    elif case.prob_side == +1:
        test_src_geo_radius = case.inner_radius
        test_tgt_geo_radius = case.outer_radius
    elif case.prob_side == "scat":
        test_src_geo_radius = case.outer_radius
        test_tgt_geo_radius = case.outer_radius
    else:
        raise ValueError("unknown problem_side")

    point_sources = make_circular_point_group(
            mesh.ambient_dim, 10, test_src_geo_radius,
            func=lambda x: x**1.5)
    test_targets = make_circular_point_group(
            mesh.ambient_dim, 20, test_tgt_geo_radius)

    np.random.seed(22)
    source_charges = np.random.randn(point_sources.shape[1])
    source_charges[-1] = -np.sum(source_charges[:-1])
    source_charges = source_charges.astype(dtype)
    assert np.sum(source_charges) < 1e-15

    source_charges_dev = cl.array.to_device(queue, source_charges)

    # }}}

    # {{{ establish BCs

    from pytential.source import PointPotentialSource
    from pytential.target import PointsTarget

    point_source = PointPotentialSource(cl_ctx, point_sources)

    pot_src = sym.IntG(
            # FIXME: qbx_forced_limit--really?
            knl, sym.var("charges"), qbx_forced_limit=None, **knl_kwargs)

    test_direct = bind((point_source, PointsTarget(test_targets)), pot_src)(
            queue, charges=source_charges_dev, **concrete_knl_kwargs)

    if case.bc_type == "dirichlet":
        bc = bind((point_source, density_discr), pot_src)(
                queue, charges=source_charges_dev, **concrete_knl_kwargs)

    elif case.bc_type == "neumann":
        bc = bind(
                (point_source, density_discr),
                sym.normal_derivative(
                    qbx.ambient_dim, pot_src, where=sym.DEFAULT_TARGET)
                )(queue, charges=source_charges_dev, **concrete_knl_kwargs)

    # }}}

    # {{{ solve

    bound_op = bind(qbx, op_u)

    rhs = bind(density_discr, op.prepare_rhs(sym.var("bc")))(queue, bc=bc)

    try:
        from pytential.solve import gmres
        gmres_result = gmres(
                bound_op.scipy_op(queue, "u", dtype, **concrete_knl_kwargs),
                rhs,
                tol=case.gmres_tol,
                progress=True,
                hard_failure=True,
                stall_iterations=50, no_progress_factor=1.05)
    except QBXTargetAssociationFailedException as e:
        bdry_vis = make_visualizer(queue, density_discr, case.target_order+3)
        bdry_vis.write_vtk_file(
                output_file_path("failed-targets-%s.vtu" % resolution),
                [("failed_targets", e.failed_target_flags),])
        raise

    logger.info("gmres state: %s", gmres_result.state)
    weighted_u = gmres_result.solution

    # }}}

    if case.prob_side != "scat":
        # {{{ error check

        points_target = PointsTarget(test_targets)
        bound_tgt_op = bind(
                (qbx, points_target),
                op.representation(sym.var("u")))

        test_via_bdry = bound_tgt_op(queue, u=weighted_u, k=case.k)

        err = test_via_bdry - test_direct

        err = err.get()
        test_direct = test_direct.get()
        test_via_bdry = test_via_bdry.get()

        # {{{ remove effect of net source charge

        if case.k == 0 and case.bc_type == "neumann" and loc_sign == -1:

            # remove constant offset in interior Laplace Neumann error
            tgt_ones = np.ones_like(test_direct)
            tgt_ones = tgt_ones/la.norm(tgt_ones)
            err = err - np.vdot(tgt_ones, err)*tgt_ones

        # }}}

        rel_err_2 = la.norm(err)/la.norm(test_direct)
        rel_err_inf = la.norm(err, np.inf)/la.norm(test_direct, np.inf)

        # }}}

        logger.info("rel_err_2: %g rel_err_inf: %g" % (rel_err_2, rel_err_inf))

    else:
        rel_err_2 = None
        rel_err_inf = None

    # {{{ test gradient

    if case.check_gradient and case.prob_side != "scat":
        bound_grad_op = bind(
                (qbx, points_target),
                op.representation(
                    sym.var("u"),
                    map_potentials=lambda pot: sym.grad(mesh.ambient_dim, pot),
                    qbx_forced_limit=None))

        # logger.info(bound_t_deriv_op.code)

        grad_from_src = bound_grad_op(
                queue, u=weighted_u, **concrete_knl_kwargs)

        grad_ref = (bind(
                (point_source, points_target),
                sym.grad(mesh.ambient_dim, pot_src)
                )(queue, charges=source_charges_dev, **concrete_knl_kwargs)
                )

        grad_err = (grad_from_src - grad_ref)

        rel_grad_err_inf = (
                la.norm(grad_err[0].get(), np.inf)
                /
                la.norm(grad_ref[0].get(), np.inf))

        logger.info("rel_grad_err_inf: %g" % rel_grad_err_inf)

    # }}}

    # {{{ test tangential derivative

    if case.check_tangential_deriv and case.prob_side != "scat":
        bound_t_deriv_op = bind(
                qbx,
                op.representation(
                    sym.var("u"),
                    map_potentials=partial(sym.tangential_derivative, 2),
                    qbx_forced_limit=loc_sign))

        # logger.info(bound_t_deriv_op.code)

        tang_deriv_from_src = bound_t_deriv_op(
                queue, u=weighted_u, **concrete_knl_kwargs).as_scalar().get()

        tang_deriv_ref = (bind(
                (point_source, density_discr),
                sym.tangential_derivative(2, pot_src)
                )(queue, charges=source_charges_dev, **concrete_knl_kwargs)
                .as_scalar().get())

        if 0:
            pt.plot(tang_deriv_ref.real)
            pt.plot(tang_deriv_from_src.real)
            pt.show()

        td_err = (tang_deriv_from_src - tang_deriv_ref)

        rel_td_err_inf = (
                la.norm(td_err, np.inf)
                / la.norm(tang_deriv_ref, np.inf))

        logger.info("rel_td_err_inf: %g" % rel_td_err_inf)

    else:
        rel_td_err_inf = None

    # }}}

    # {{{ any-D file plotting

    if visualize:
        bdry_vis = make_visualizer(queue, density_discr, case.target_order+3)

        bdry_normals = (
                bind(
                    density_discr, sym.normal(qbx.ambient_dim))(queue)
                ).as_vector(dtype=object)

        sym_sqrt_j = sym.sqrt_jac_q_weight(density_discr.ambient_dim)
        u = bind(density_discr, sym.var("u")/sym_sqrt_j)(queue, u=weighted_u)

        bdry_vis.write_vtk_file(
                output_file_path("source-%s.vtu" % resolution),
                [
                    ("u", u),
                    ("bc", bc),
                    # ("bdry_normals", bdry_normals),
                    ])

        from sumpy.visualization import make_field_plotter_from_bbox  # noqa
        from meshmode.mesh.processing import find_bounding_box

        vis_grid_spacing = (0.1, 0.1, 0.1)[:qbx.ambient_dim]
        if hasattr(case, "vis_grid_spacing"):
            vis_grid_spacing = case.vis_grid_spacing
        vis_extend_factor = 0.2
        if hasattr(case, "vis_extend_factor"):
            vis_grid_spacing = case.vis_grid_spacing

        fplot = make_field_plotter_from_bbox(
                find_bounding_box(mesh),
                h=vis_grid_spacing,
                extend_factor=vis_extend_factor)

        qbx_tgt_tol = qbx.copy(target_association_tolerance=0.15)
        from pytential.target import PointsTarget

        try:
            solved_pot = bind(
                    (qbx_tgt_tol, PointsTarget(fplot.points)),
                    op.representation(sym.var("u"))
                    )(queue, u=weighted_u, k=case.k)
        except QBXTargetAssociationFailedException as e:
            fplot.write_vtk_file(
                    output_file_path("failed-targets.vts"),
                    [
                        ("failed_targets", e.failed_target_flags.get(queue))
                        ])
            raise

        from sumpy.kernel import LaplaceKernel
        ones_density = density_discr.zeros(queue)
        ones_density.fill(1)
        indicator_func = bind(
                (
                    qbx_tgt_tol, PointsTarget(fplot.points)),
                -sym.D(
                    LaplaceKernel(density_discr.ambient_dim),
                    sym.var("sigma"),
                    qbx_forced_limit=None))
        indicator = indicator_func(queue, sigma=ones_density).get()

        solved_pot = solved_pot.get()

        true_pot = bind((point_source, PointsTarget(fplot.points)), pot_src)(
                queue, charges=source_charges_dev, **concrete_knl_kwargs).get()

        # fplot.show_scalar_in_mayavi(solved_pot.real, max_val=5)
        if case.prob_side == "scat":
            fplot.write_vtk_file(
                    output_file_path("potential-%s.vts" % resolution),
                    [
                        ("pot_scattered", solved_pot),
                        ("pot_incoming", -true_pot),
                        ("indicator", indicator),
                        ]
                    )
        else:
            fplot.write_vtk_file(
                    output_file_path("potential-%s.vts" % resolution),
                    [
                        ("solved_pot", solved_pot),
                        ("true_pot", true_pot),
                        ("indicator", indicator),
                        ]
                    )

    # }}}

    class Result(Record):
        pass

    return Result(
            h_max=qbx.h_max,
            rel_err_2=rel_err_2,
            rel_err_inf=rel_err_inf,
            rel_td_err_inf=rel_td_err_inf,
            gmres_result=gmres_result)

# }}}


# {{{ test frontend

SphereHelmholtzDirichletTestCase = SphereIntEqTestCase(5, "dirichlet", +1)


@pytest.mark.parametrize(
        "case",
        (SphereHelmholtzDirichletTestCase,))
def test_integral_equation(ctx_getter, case, visualize=False):
    cl_ctx = ctx_getter()
    queue = cl.CommandQueue(cl_ctx)

    from pytools.convergence import EOCRecorder
    print("qbx_order: %d, %s" % (case.qbx_order, case))

    lpot_kwargs = {"_use_target_specific_qbx": True}

    eoc_rec_target = EOCRecorder()
    eoc_rec_td = EOCRecorder()

    have_error_data = False
    for resolution in case.resolutions:
        result = run_int_eq_test(
                cl_ctx, queue, case, resolution,
                visualize=visualize, lpot_kwargs=lpot_kwargs)

        if result.rel_err_2 is not None:
            have_error_data = True
            eoc_rec_target.add_data_point(result.h_max, result.rel_err_2)

        if result.rel_td_err_inf is not None:
            eoc_rec_td.add_data_point(result.h_max, result.rel_td_err_inf)

    if case.bc_type == "dirichlet":
        tgt_order = case.qbx_order
    elif case.bc_type == "neumann":
        tgt_order = case.qbx_order-1
    else:
        assert False

    if have_error_data:
        print("TARGET ERROR:")
        print(eoc_rec_target)
        assert eoc_rec_target.order_estimate() > tgt_order - 1.3

        if case.check_tangential_deriv:
            print("TANGENTIAL DERIVATIVE ERROR:")
            print(eoc_rec_td)
            assert eoc_rec_td.order_estimate() > tgt_order - 2.3

# }}}


def plane_lpot_source(queue, lpot_kwargs):
    cl_ctx = queue.context
    lpot_kwargs = lpot_kwargs.copy()
    case = BetterplaneIntEqTestCase(20, "dirichlet", +1)
    return lpot_source_from_case(cl_ctx, queue, case, 0.25, lpot_kwargs)


def sphere_lpot_source(queue, lpot_kwargs):
    cl_ctx = queue.context
    lpot_kwargs = lpot_kwargs.copy()
    case = SphereIntEqTestCase(20, "dirichlet", +1)
    return lpot_source_from_case(cl_ctx, queue, case, 1, lpot_kwargs)


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    pytest.main([__file__])