Merge branch 'main' into stokes_biharmonic2

doc/conf.py

@@ -24,7 +24,6 @@
 
 intersphinx_mapping = {
     "https://docs.python.org/3/": None,
-    "https://documen.tician.de/boxtree/": None,
     "https://numpy.org/doc/stable/": None,
     "https://scipy.github.io/devdocs/": None,
     "https://documen.tician.de/arraycontext/": None,
@@ -34,5 +33,6 @@
     "https://documen.tician.de/pytools/": None,
     "https://documen.tician.de/pymbolic/": None,
     "https://documen.tician.de/loopy/": None,
+    "https://documen.tician.de/boxtree/": None,
     "https://documen.tician.de/sumpy/": None,
     }

doc/index.rst

@@ -20,7 +20,6 @@ Contents
     discretization
     symbolic
     linalg
-    tools
     misc
     qbx
     🚀 Github <https://github.com/inducer/pytential>

doc/linalg.rst

@@ -10,6 +10,11 @@ scheme is used:
   proxy-based skeletonization of the direct solver algorithms. Clusters
   are represented by a :class:`~pytential.linalg.TargetAndSourceClusterList`.
 
+GMRES
+-----
+
+.. automodule:: pytential.linalg.gmres
+
 .. _direct_solver:
 
 Hierarchical Direct Solver
@@ -28,4 +33,4 @@ Internal Functionality
 
 .. automodule:: pytential.linalg.direct_solver_symbolic
 
-.. vim: sw=4:tw=75
+.. vim: sw=4:tw=75:fdm=marker

examples/helmholtz-dirichlet.py

@@ -133,7 +133,7 @@ def u_incoming_func(x):
 
     bvp_rhs = bind(places, sqrt_w*sym.var("bc"))(actx, bc=bc)
 
-    from pytential.solve import gmres
+    from pytential.linalg.gmres import gmres
     gmres_result = gmres(
             bound_op.scipy_op(actx, sigma_sym.name, dtype=np.complex128, k=k),
             bvp_rhs, tol=1e-8, progress=True,

examples/laplace-dirichlet-3d.py

@@ -119,7 +119,7 @@ def u_incoming_func(x):
     bc = u_incoming_func(nodes)
     bvp_rhs = bind(places, sqrt_w*sym.var("bc"))(actx, bc=bc)
 
-    from pytential.solve import gmres
+    from pytential.linalg.gmres import gmres
     gmres_result = gmres(
             bound_op.scipy_op(actx, "sigma", dtype=np.float64),
             bvp_rhs, tol=1e-14, progress=True,

experiments/cahn-hilliard.py

@@ -73,7 +73,7 @@ def g(xvec):
         -g(nodes),
         ])
 
-    from pytential.solve import gmres
+    from pytential.linalg.gmres import gmres
     gmres_result = gmres(
             bound_op.scipy_op(queue, "sigma", dtype=np.complex128),
             bc, tol=1e-8, progress=True,

experiments/find-photonic-mode-sk.py

@@ -107,7 +107,7 @@ def find_mode():
     bound_op = bind(qbx, op)
 
     def muller_solve_func(ne):
-        from pytential.solve import gmres
+        from pytential.linalg.gmres import gmres
         gmres_result = gmres(
                 bound_op.scipy_op(queue, "u",
                     np.complex128, ne=ne, **base_context),

experiments/layerpot-coverage.py

@@ -82,7 +82,7 @@ def reference_solu(rvec):
             op.operator(sym.var('sigma')))
     rhs = bind(qbx.density_discr, op.prepare_rhs(sym.var("bc")))(queue, bc=bvals)
 
-    from pytential.solve import gmres
+    from pytential.linalg.gmres import gmres
     gmres_result = gmres(
             bound_op.scipy_op(queue, "sigma", dtype=np.float64),
             rhs,

experiments/maxwell.py

@@ -130,7 +130,7 @@ def green3m(x,y,z,source,strength,k):
         # magnetic field corresponding to dyadic green's function
         # due to monochromatic electric dipole located at "source".
         # "strength" is the the intensity of the dipole.
-        #  H = curl((I + Hess)(exp(ikr)/r) dot (strength)) = 
+        #  H = curl((I + Hess)(exp(ikr)/r) dot (strength)) =
         #  strength \cross \grad (exp(ikr)/r)
         #
             dx = x - source[0]
@@ -158,7 +158,7 @@ def dipole3e(x,y,z,source,strength,k):
             evec = evec*1j*k
             hvec = green3m(x,y,z,source,strength,k)
             return evec,hvec
-            
+
         def dipole3m(x,y,z,source,strength,k):
         #
         #  evalaute electric and magnetic field due
@@ -168,7 +168,7 @@ def dipole3m(x,y,z,source,strength,k):
             hvec = green3e(x,y,z,source,strength,k)
             hvec = -hvec*1j*k
             return evec,hvec
-            
+
 
         def dipole3eall(x,y,z,sources,strengths,k):
             ns = len(strengths)
@@ -187,7 +187,7 @@ def dipole3eall(x,y,z,sources,strengths,k):
 #        source[1] =-0.03
 #        source[2] = 0.02
         strength = np.ones(3)
-       
+
 #        evec = cl.array.to_device(queue,np.zeros((3,len(nodes[0])),dtype=np.complex128))
 #        hvec = cl.array.to_device(queue,np.zeros((3,len(nodes[0])),dtype=np.complex128))
 
@@ -210,7 +210,7 @@ def dipole3eall(x,y,z,sources,strengths,k):
 
         bound_op = bind(qbx, sym_operator)
 
-        from pytential.solve import gmres
+        from pytential.linalg.gmres import gmres
         if 0:
             gmres_result = gmres(
                 bound_op.scipy_op(queue, "sigma", dtype=np.complex128, k=k),

experiments/maxwell_sphere.py

@@ -208,7 +208,7 @@ def dipole3eall(x,y,z,sources,strengths,k):
 
         bound_op = bind(qbx, sym_operator)
 
-        from pytential.solve import gmres
+        from pytential.linalg.gmres import gmres
         if 1:
             gmres_result = gmres(
                 bound_op.scipy_op(queue, "sigma", dtype=np.complex128, k=k),

experiments/poisson.py

@@ -248,7 +248,7 @@ def get_kernel():
 
     bvp_rhs = bind(bdry_discr, op.prepare_rhs(sym.var("bc")))(queue, bc=bvp_bc)
 
-    from pytential.solve import gmres
+    from pytential.linalg.gmres import gmres
     gmres_result = gmres(
             bound_op.scipy_op(queue, "sigma", dtype=np.float64),
             bvp_rhs, tol=1e-14, progress=True,

experiments/stokes-2d-interior.py

@@ -127,7 +127,7 @@ def couette_soln(x, y, dp, h):
     # Get rhs vector
     bvp_rhs = bind(qbx, sqrt_w*sym.make_sym_vector("bc",dim))(queue, bc=bc)
 
-    from pytential.solve import gmres
+    from pytential.linalg.gmres import gmres
     gmres_result = gmres(
              bound_op.scipy_op(queue, "sigma", np.float64, mu=mu, normal=normal),
              bvp_rhs, tol=1e-9, progress=True,

experiments/two-domain-helmholtz.py

@@ -136,7 +136,7 @@ def main():
 
         bvp_rhs[i_bc] *= sqrt_w
 
-    from pytential.solve import gmres
+    from pytential.linalg.gmres import gmres
     gmres_result = gmres(
             bound_pde_op.scipy_op(queue, "unknown", dtype=np.complex128,
                 domains=[sym.DEFAULT_TARGET]*2, K0=K0, K1=K1),

pytential/linalg/__init__.py

@@ -23,18 +23,27 @@
 from pytential.linalg.utils import (
         IndexList, TargetAndSourceClusterList,
         make_index_list, make_index_cluster_cartesian_product,
+        interp_decomp,
         )
 from pytential.linalg.proxy import (
-        ProxyClusterGeometryData,
+        ProxyClusterGeometryData, ProxyPointTarget, ProxyPointSource,
         ProxyGeneratorBase, ProxyGenerator, QBXProxyGenerator,
         partition_by_nodes, gather_cluster_neighbor_points,
         )
+from pytential.linalg.skeletonization import (
+        SkeletonizationWrangler, make_skeletonization_wrangler,
+        SkeletonizationResult, skeletonize_by_proxy,
+        )
 
 __all__ = (
     "IndexList", "TargetAndSourceClusterList",
     "make_index_list", "make_index_cluster_cartesian_product",
+    "interp_decomp",
 
-    "ProxyClusterGeometryData",
+    "ProxyClusterGeometryData", "ProxyPointTarget", "ProxyPointSource",
     "ProxyGeneratorBase", "ProxyGenerator", "QBXProxyGenerator",
     "partition_by_nodes", "gather_cluster_neighbor_points",
+
+    "SkeletonizationWrangler", "make_skeletonization_wrangler",
+    "SkeletonizationResult", "skeletonize_by_proxy",
 )

pytential/linalg/direct_solver_symbolic.py

@@ -20,6 +20,8 @@
 THE SOFTWARE.
 """
 
+from pytools.obj_array import make_obj_array
+
 from pytential.symbolic.mappers import (
         IdentityMapper, OperatorCollector, LocationTagger)
 
@@ -30,6 +32,39 @@
 """
 
 
+# {{{ utils
+
+class PROXY_SKELETONIZATION_SOURCE:             # noqa: N801
+    pass
+
+
+class PROXY_SKELETONIZATION_TARGET:             # noqa: N801
+    pass
+
+
+def prepare_expr(places, exprs, auto_where=None):
+    from pytential.symbolic.execution import _prepare_expr
+    return make_obj_array([
+        _prepare_expr(places, expr, auto_where=auto_where)
+        for expr in exprs])
+
+
+def prepare_proxy_expr(places, exprs, auto_where=None):
+    def _prepare_expr(expr):
+        # remove all diagonal / non-operator terms in the expression
+        expr = IntGTermCollector()(expr)
+        # ensure all IntGs remove all the kernel derivatives
+        expr = KernelTransformationRemover()(expr)
+        # ensure all IntGs have their source and targets set
+        expr = DOFDescriptorReplacer(auto_where[0], auto_where[1])(expr)
+
+        return expr
+
+    return make_obj_array([_prepare_expr(expr) for expr in exprs])
+
+# }}}
+
+
 # {{{ KernelTransformationRemover
 
 class KernelTransformationRemover(IdentityMapper):