diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 6a1c1be779afe96f06c9c59a747b9302267b79e1..ebf121320d6e6c910e9d05f849913191592c0531 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -64,7 +64,7 @@ Python 3 POCL Examples: - test -n "$SKIP_EXAMPLES" && exit - export PY_EXE=python3 - export PYOPENCL_TEST=portable:pthread - - export EXTRA_INSTALL="Cython pybind11 numpy mako pyvisfile matplotlib" + - export EXTRA_INSTALL="Cython pybind11 numpy mako git+git://github.com/inducer/pytools pyvisfile matplotlib" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-py-project-and-run-examples.sh - ". ./build-py-project-and-run-examples.sh" tags: diff --git a/.test-conda-env-py3-macos.yml b/.test-conda-env-py3-macos.yml index cbf0efad3ebc7688138e04becf0fbccdc79b4277..b4057a83810cc455a597561cbcd40e16133b3203 100644 --- a/.test-conda-env-py3-macos.yml +++ b/.test-conda-env-py3-macos.yml @@ -22,6 +22,7 @@ dependencies: - pip - pip: + - git+https://github.com/inducer/pytools - git+https://gitlab.tiker.net/inducer/boxtree - git+https://github.com/inducer/pymbolic - git+https://github.com/inducer/loopy diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index 750b00724d10913f44a8070d577ca040fe1a6251..81849ec0c582e2d89f5540e157897d4eea631577 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -18,6 +18,7 @@ dependencies: - pip - pip: + - git+https://github.com/inducer/pytools - git+https://gitlab.tiker.net/inducer/boxtree - git+https://github.com/inducer/pymbolic - git+https://github.com/inducer/loopy diff --git a/pytential/solve.py b/pytential/solve.py index 82d6f68ff18e25fcbbb1e975b1ce8dc4af12eeba..47d5b907cac3034377fa5256f641ead0bbe15671 100644 --- a/pytential/solve.py +++ b/pytential/solve.py @@ -301,7 +301,7 @@ def gmres(op, rhs, restart=None, tol=None, x0=None, inner_product=None, maxiter=None, hard_failure=None, no_progress_factor=None, stall_iterations=None, - callback=None, progress=False): + callback=None, progress=False, require_monotonicity=True): """Solve a linear system Ax=b by means of GMRES with restarts. @@ -338,7 +338,8 @@ def gmres(op, rhs, restart=None, tol=None, x0=None, dot=inner_product, maxiter=maxiter, hard_failure=hard_failure, no_progress_factor=no_progress_factor, - stall_iterations=stall_iterations, callback=callback) + stall_iterations=stall_iterations, callback=callback, + require_monotonicity=require_monotonicity) return result.copy(solution=chopper.chop(result.solution)) diff --git a/pytential/symbolic/pde/scalar.py b/pytential/symbolic/pde/scalar.py index 882b1b3c53cb53830bb94796b887a71494fdae32..763a403f6a497e07583a3657cc287d05838e4ad8 100644 --- a/pytential/symbolic/pde/scalar.py +++ b/pytential/symbolic/pde/scalar.py @@ -26,6 +26,7 @@ __doc__ = """ .. autoclass:: L2WeightedPDEOperator .. autoclass:: DirichletOperator .. autoclass:: NeumannOperator +.. autoclass:: BiharmonicClampedPlateOperator """ @@ -33,6 +34,7 @@ from pytential import sym from pytential.symbolic.primitives import ( cse, sqrt_jac_q_weight, QWeight, area_element) +from sumpy.kernel import DirectionalSourceDerivative import numpy as np # noqa @@ -64,6 +66,13 @@ class L2WeightedPDEOperator(object): def prepare_rhs(self, b): return self.get_sqrt_weight()*b + def get_density_var(self, name): + """ + Returns a symbolic variable/array corresponding to the density with the + given name. + """ + return sym.var(name) + # }}} @@ -305,6 +314,90 @@ class NeumannOperator(L2WeightedPDEOperator): + ones_contribution )) + +class BiharmonicClampedPlateOperator: + r"""IE operator and field representation for solving clamped plate Biharmonic + equation where, + + .. math:: + \begin{align*} + \Delta^2 u &= 0 \text{ on } D \\ + u &= g_1 \text{ on } \delta D \\ + \frac{\partial u}{\partial \nu} &= g_2 \text{ on } \delta D. + \end{align*} + + This operator assumes that the boundary data :math:`g_1, g_2` are + represented as column vectors and vertically stacked. + + .. note :: This operator supports only interior problem. + + Ref: Farkas, Peter. Mathematical foundations for fast algorithms for the + biharmonic equation. Technical Report 765, Department of Computer Science, + Yale University, 1990. + + .. automethod:: representation + .. automethod:: operator + """ + + def __init__(self, knl, loc_sign): + self.knl = knl + self.loc_sign = loc_sign + if loc_sign != -1: + raise ValueError("loc_sign!=-1 is not supported.") + + def prepare_rhs(self, b): + return b + + def get_density_var(self, name): + """ + Returns a symbolic variable of length 2 corresponding to the density with + the given name. + """ + return sym.make_sym_vector(name, 2) + + def representation(self, + sigma, map_potentials=None, qbx_forced_limit=None): + + if map_potentials is None: + def map_potentials(x): # pylint:disable=function-redefined + return x + + def dv(knl): + return DirectionalSourceDerivative(knl, "normal_dir") + + def dt(knl): + return DirectionalSourceDerivative(knl, "tangent_dir") + + normal_dir = sym.normal(2).as_vector() + tangent_dir = np.array([-normal_dir[1], normal_dir[0]]) + + k1 = map_potentials(sym.S(dv(dv(dv(self.knl))), sigma[0], + kernel_arguments={"normal_dir": normal_dir}, + qbx_forced_limit=qbx_forced_limit)) + \ + 3*map_potentials(sym.S(dv(dt(dt(self.knl))), sigma[0], + kernel_arguments={"normal_dir": normal_dir, + "tangent_dir": tangent_dir}, + qbx_forced_limit=qbx_forced_limit)) + + k2 = -map_potentials(sym.S(dv(dv(self.knl)), sigma[1], + kernel_arguments={"normal_dir": normal_dir}, + qbx_forced_limit=qbx_forced_limit)) + \ + map_potentials(sym.S(dt(dt(self.knl)), sigma[1], + kernel_arguments={"tangent_dir": tangent_dir}, + qbx_forced_limit=qbx_forced_limit)) + + return k1 + k2 + + def operator(self, sigma): + """ + Returns the two second kind integral equations. + """ + rep = self.representation(sigma, qbx_forced_limit='avg') + rep_diff = sym.normal_derivative(2, rep) + int_eq1 = sigma[0]/2 + rep + int_eq2 = -sym.mean_curvature(2)*sigma[0] + sigma[1]/2 + rep_diff + return np.array([int_eq1, int_eq2]) + # }}} diff --git a/requirements.txt b/requirements.txt index dd15a69ebf6dade1ca2fb2df0bdc3644b6b8c6d8..625deb28d7e5a04253bc3ffb5c12b2ba263362b5 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,5 @@ numpy +git+git://github.com/inducer/pytools git+git://github.com/inducer/pymbolic sympy==1.1.1 git+https://github.com/inducer/modepy diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py index b00efdd7c0b9e1111e3dd490418dd99e3fddfb0c..82faf4cc118e97e1c8dd3464a711cec484ea77a4 100644 --- a/test/test_scalar_int_eq.py +++ b/test/test_scalar_int_eq.py @@ -40,6 +40,7 @@ from meshmode.discretization.visualization import make_visualizer from sumpy.symbolic import USE_SYMENGINE from pytential import bind, sym from pytential.qbx import QBXTargetAssociationFailedException +from sumpy.kernel import LaplaceKernel, HelmholtzKernel, BiharmonicKernel import logging logger = logging.getLogger(__name__) @@ -65,10 +66,6 @@ def make_circular_point_group(ambient_dim, npoints, radius, class IntEqTestCase: - @property - def default_helmholtz_k(self): - raise NotImplementedError - @property def name(self): raise NotImplementedError @@ -81,27 +78,37 @@ class IntEqTestCase: def target_order(self): raise NotImplementedError - def __init__(self, helmholtz_k, bc_type, prob_side): + def __init__(self, knl_class_or_helmholtz_k, bc_type, + prob_side, knl_kwargs={}, **kwargs): """ :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 + if not isinstance(knl_class_or_helmholtz_k, type): + if knl_class_or_helmholtz_k == 0: + knl_class = LaplaceKernel + else: + knl_kwargs = {"k": knl_class_or_helmholtz_k} + knl_class = HelmholtzKernel + else: + knl_class = knl_class_or_helmholtz_k + self.knl_class = knl_class + self.knl_kwargs = knl_kwargs + self.knl_kwargs_syms = dict((k, sym.var(k)) for k in self.knl_kwargs.keys()) + for k, v in self.knl_kwargs.items(): + setattr(self, k, v) + for k, v in kwargs.items(): + setattr(self, k, v) def __str__(self): + kwargs = ", ".join("%s: %s" % (k, v) for k, v in self.knl_kwargs.items()) + if kwargs: + kwargs = " , " + kwargs 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)) + "knl_class: %s, qbx_order: %d, target_order: %d%s" + % (self.name, self.bc_type, self.prob_side, self.knl_class, + self.qbx_order, self.target_order, kwargs)) fmm_backend = "sumpy" gmres_tol = 1e-14 @@ -323,7 +330,6 @@ class ElliptiplaneIntEqTestCase(IntEqTestCase): class BetterplaneIntEqTestCase(IntEqTestCase): name = "betterplane" - default_helmholtz_k = 20 resolutions = [0.2] # refine_on_helmholtz_k = False @@ -469,7 +475,8 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): fmm_backend=case.fmm_backend, **qbx_lpot_kwargs) if case.use_refinement: - if case.k != 0 and getattr(case, "refine_on_helmholtz_k", True): + if case.knl_class == HelmholtzKernel and \ + getattr(case, "refine_on_helmholtz_k", True): refiner_extra_kwargs["kernel_length_scale"] = 5/case.k if hasattr(case, "scaled_max_curvature_threshold"): @@ -537,18 +544,14 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): # {{{ 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 = {} + DirichletOperator, + NeumannOperator, + BiharmonicClampedPlateOperator, + ) + + knl = case.knl_class(mesh.ambient_dim) + knl_kwargs_syms = case.knl_kwargs_syms + concrete_knl_kwargs = case.knl_kwargs if knl.is_complex_valued: dtype = np.complex128 @@ -559,14 +562,16 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): if case.bc_type == "dirichlet": op = DirichletOperator(knl, loc_sign, use_l2_weighting=True, - kernel_arguments=knl_kwargs) + kernel_arguments=knl_kwargs_syms) elif case.bc_type == "neumann": op = NeumannOperator(knl, loc_sign, use_l2_weighting=True, - use_improved_operator=False, kernel_arguments=knl_kwargs) + use_improved_operator=False, kernel_arguments=knl_kwargs_syms) + elif case.bc_type == "clamped_plate": + op = BiharmonicClampedPlateOperator(knl, loc_sign) else: assert False - op_u = op.operator(sym.var("u")) + op_u = op.operator(op.get_density_var("u")) # }}} @@ -609,7 +614,7 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): pot_src = sym.IntG( # FIXME: qbx_forced_limit--really? - knl, sym.var("charges"), qbx_forced_limit=None, **knl_kwargs) + knl, sym.var("charges"), qbx_forced_limit=None, **knl_kwargs_syms) test_direct = bind((point_source, PointsTarget(test_targets)), pot_src)( queue, charges=source_charges_dev, **concrete_knl_kwargs) @@ -625,12 +630,22 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): qbx.ambient_dim, pot_src, dofdesc=sym.DEFAULT_TARGET) )(queue, charges=source_charges_dev, **concrete_knl_kwargs) + elif case.bc_type == "clamped_plate": + bc_u = bind((point_source, density_discr), pot_src)( + queue, charges=source_charges_dev, **concrete_knl_kwargs) + bc_du = bind( + (point_source, density_discr), + sym.normal_derivative( + qbx.ambient_dim, pot_src, dofdesc=sym.DEFAULT_TARGET) + )(queue, charges=source_charges_dev, **concrete_knl_kwargs) + bc = [bc_u, bc_du] + # }}} # {{{ solve bound_op = bind(qbx, op_u) - rhs = bind(density_discr, op.prepare_rhs(sym.var("bc")))(queue, bc=bc) + rhs = bind(density_discr, op.prepare_rhs(op.get_density_var("bc")))(queue, bc=bc) try: from pytential.solve import gmres @@ -660,7 +675,7 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): from sumpy.tools import build_matrix mat = build_matrix( bound_op.scipy_op( - queue, arg_name="u", dtype=dtype, k=case.k)) + queue, arg_name="u", dtype=dtype, **concrete_knl_kwargs)) w, v = la.eig(mat) if 0: pt.imshow(np.log10(1e-20+np.abs(mat))) @@ -678,9 +693,9 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): points_target = PointsTarget(test_targets) bound_tgt_op = bind((qbx, points_target), - op.representation(sym.var("u"))) + op.representation(op.get_density_var("u"))) - test_via_bdry = bound_tgt_op(queue, u=weighted_u, k=case.k) + test_via_bdry = bound_tgt_op(queue, u=weighted_u, **concrete_knl_kwargs) err = test_via_bdry - test_direct @@ -690,7 +705,8 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): # {{{ remove effect of net source charge - if case.k == 0 and case.bc_type == "neumann" and loc_sign == -1: + if case.knl_class == LaplaceKernel and case.bc_type == "neumann" \ + and loc_sign == -1: # remove constant offset in interior Laplace Neumann error tgt_ones = np.ones_like(test_direct) @@ -715,7 +731,7 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): if case.check_gradient and case.prob_side != "scat": bound_grad_op = bind((qbx, points_target), op.representation( - sym.var("u"), + op.get_density_var("u"), map_potentials=lambda pot: sym.grad(mesh.ambient_dim, pot), qbx_forced_limit=None)) @@ -745,7 +761,7 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): if case.check_tangential_deriv and case.prob_side != "scat": bound_t_deriv_op = bind(qbx, op.representation( - sym.var("u"), + op.get_density_var("u"), map_potentials=lambda pot: sym.tangential_derivative(2, pot), qbx_forced_limit=loc_sign)) @@ -785,7 +801,8 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): .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) + u = bind(density_discr, op.get_density_var("u")/sym_sqrt_j)(queue, + u=weighted_u) bdry_vis.write_vtk_file("source-%s.vtu" % resolution, [ ("u", u), @@ -814,8 +831,8 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): try: solved_pot = bind( (qbx_tgt_tol, PointsTarget(fplot.points)), - op.representation(sym.var("u")) - )(queue, u=weighted_u, k=case.k) + op.representation(op.get_density_var("u")) + )(queue, u=weighted_u, **concrete_knl_kwargs) except QBXTargetAssociationFailedException as e: fplot.write_vtk_file( "failed-targets.vts", @@ -824,13 +841,12 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): ]) raise - from sumpy.kernel import LaplaceKernel ones_density = density_discr.zeros(queue) ones_density.fill(1) indicator = bind( (qbx_tgt_tol, PointsTarget(fplot.points)), -sym.D(LaplaceKernel(density_discr.ambient_dim), - sym.var("sigma"), + op.get_density_var("sigma"), qbx_forced_limit=None))( queue, sigma=ones_density).get() @@ -875,17 +891,27 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): # }}} -# {{{ test frontend - -@pytest.mark.parametrize("case", [ - EllipseIntEqTestCase(helmholtz_k=helmholtz_k, bc_type=bc_type, +cases = [ + EllipseIntEqTestCase(helmholtz_k, bc_type=bc_type, prob_side=prob_side) for helmholtz_k in [0, 1.2] for bc_type in ["dirichlet", "neumann"] for prob_side in [-1, +1] - ]) +] + +cases += [ + EllipseIntEqTestCase(BiharmonicKernel, bc_type="clamped_plate", + prob_side=-1, fmm_backend=None), + EllipseIntEqTestCase(BiharmonicKernel, bc_type="clamped_plate", + prob_side=-1, fmm_backend="sumpy", fmm_order=15), +] + +# {{{ test frontend + + +@pytest.mark.parametrize("case", cases) # Sample test run: -# 'test_integral_equation(cl._csc, EllipseIntEqTestCase(0, "dirichlet", +1), visualize=True)' # noqa: E501 +# 'test_integral_equation(cl._csc, EllipseIntEqTestCase(LaplaceKernel, "dirichlet", +1), visualize=True)' # noqa: E501 def test_integral_equation(ctx_factory, case, visualize=False): logging.basicConfig(level=logging.INFO) @@ -921,6 +947,8 @@ def test_integral_equation(ctx_factory, case, visualize=False): tgt_order = case.qbx_order elif case.bc_type == "neumann": tgt_order = case.qbx_order-1 + elif case.bc_type == "clamped_plate": + tgt_order = case.qbx_order else: assert False