From dfad02861f1e893f2149c9e3c49fe5a32f63c4ea Mon Sep 17 00:00:00 2001 From: benSepanski Date: Thu, 18 Feb 2021 12:17:07 -0600 Subject: [PATCH 01/20] Set up connections, make interface for potentials --- firedrake/layer_potentials.py | 479 ++++++++++++++++++++++++++++++++++ 1 file changed, 479 insertions(+) diff --git a/firedrake/layer_potentials.py b/firedrake/layer_potentials.py index 979a40f82f..b7f1ac496a 100644 --- a/firedrake/layer_potentials.py +++ b/firedrake/layer_potentials.py @@ -1,6 +1,8 @@ import numpy as np from firedrake.pointwise_operators import AbstractExternalOperator from firedrake.utils import cached_property +from firedrake.functionspaceimpl import WithGeometry +import firedrake.function from firedrake import Function, FunctionSpace, interpolate, Interpolator, \ SpatialCoordinate, VectorFunctionSpace, TensorFunctionSpace, project @@ -12,6 +14,7 @@ from ufl.tensors import as_tensor from pyop2.datatypes import ScalarType +from warnings import warn __all__ = ("DoubleLayerPotential", "PytentialLayerOperation", @@ -20,6 +23,482 @@ "VolumePotential",) +class Potential(AbstractExternalOperator): + r""" + This is a function which represents a potential + computed on a firedrake function. + + For a firedrake function :math:`u:\Omega\to\mathbb C^n` + with :math:`\Omega \subseteq \mathbb R^m`, + the potential :math:`P(u):\Omega \to\mathbb C^n` is a + convolution of :math:`u` against some kernel function. + More concretely, given a source region :math:`\Gamma \subseteq \Omega` + a target region :math:`\Sigma \subseteq \Omega`, and + a kernel function :math:`K`, we define + + .. math:: + + P(u)(x) = \begin{cases} + \int_\Gamma K(x, y) u(y) \,dy & x \in \Sigma + 0 & x \notin \Sigma + \end{cases} + """ + _external_operator_type = 'GLOBAL' + + def __init__(self, density, **kwargs): + """ + :arg density: A :mod:`firedrake` :class:`firedrake.function.Function` + or UFL expression which represents the density :math:`u` + :kwarg connection: A :class:`PotentialEvaluationLibraryConnection` + :kwarg potential_operator: The external potential evaluation library + bound operator. + """ + # super + AbstractExternalOperator.__init__(self, density, **kwargs) + + # FIXME + for order in self.derivatives: + assert order == 1, "Assumes self.derivatives = (1,..,1)" + + # Get connection & bound op and validate + connection = kwargs.get("connection", None) + if connection is None: + raise ValueError("Missing kwarg 'connection'") + if not isinstance(connection, PotentialEvaluationLibraryConnection): + raise TypeError("connection must be of type " + "PotentialEvaluationLibraryConnection, not %s." + % type(connection)) + self.connection = connection + + self.potential_operator = kwargs.get("potential_operator", None) + if self.potential_operator is None: + raise ValueError("Missing kwarg 'potential_operator'") + + # Get function space and validate aginst bound op + function_space = self.ufl_function_space() + assert isinstance(function_space, WithGeometry) # sanity check + if function_space is not self.potential_operator.function_space(): + raise ValueError("function_space must be same object as " + "potential_operator.function_space().") + + # Make sure density is a member of our function space, if it is + if isinstance(density, Function): + if density.function_space() is not function_space: + raise ValueError("density.function_space() must be the " + "same as function_space") + # Make sure the shapes match, at least + elif density.shape != function_space.shape: + raise ValueError("Shape mismatch between function_space and " + "density. %s != %s." % + (density.shape, function_space.shape)) + + @cached_property + def _evaluator(self): + return PotentialEvaluator(self, + self.density, + self.connection, + self.potential_operator) + + def _evaluate(self): + raise NotImplementedError + + def _compute_derivatives(self, continuity_tolerance=None): + raise NotImplementedError + + def _evaluate_action(self, *args): + return self._evaluator._evaluate_action() + + def _evaluate_adjoint_action(self, *args): + return self._evaluator._evaluate_action() + + def evaluate_adj_component_control(self, x, idx): + derivatives = tuple(dj + int(idx == j) for j, dj in enumerate(self.derivatives)) + dN = self._ufl_expr_reconstruct_(*self.ufl_operands, derivatives=derivatives) + return dN._evaluate_adjoint_action((x.function,)).vector() + + def evaluate_adj_component(self, x, idx): + print(x, type(x)) + raise NotImplementedError + + +class PotentialEvaluator: + """ + Evaluates a potential + """ + def __init__(self, density, connection, potential_operator): + """ + :arg density: The UFL/firedrake density function + :arg connection: A :class:`PotentialEvaluationLibraryConnection` + :arg potential_operator: The external potential evaluation library + bound operator. + """ + self.density = density + self.connection = connection + self.potential_operator = potential_operator + + def _eval_potential_operator(self, action_coefficients, out=None): + """ + Evaluate the potential on the action coefficients and return. + If *out* is not *None*, stores the result in *out* + + :arg action_coefficients: A :class:`~firedrake.function.Function` + to evaluate the potential at + :arg out: If not *None*, then a :class:`~firedrake.function.Function` + in which to store the evaluated potential + :return: *out* if it is not *None*, otherwise a new + :class:`firedrake.function.Function` storing the evaluated + potential + """ + density = self.connection.from_firedrake(action_coefficients) + potential = self.potential_operator(density) + return self.connection.to_firedrake(potential, out=out) + + def _evaluate(self): + """ + Evaluate P(self.density) into self + """ + return self._eval_potential_operator(self.density, out=self) + + def _compute_derivatives(self): + """ + Return a function + Derivative(P, self.derivatives, self.density) + """ + # FIXME : Assumes derivatives are Jacobian + return self._eval_potential_operator + + def _evaluate_action(self, action_coefficients): + """ + Evaluate derivatives of layer potential at action coefficients. + i.e. Derivative(P, self.derivatives, self.density) evaluated at + the action coefficients and store into self + """ + operator = self._compute_derivatives() + return operator(action_coefficients, out=self) + + +class PotentialEvaluationLibraryConnection: + """ + A connection to an external library for potential evaluation + """ + def __init__(self, function_space, warn_if_cg=True): + """ + Initialize self to work on the function space + + :arg function_space: The :mod:`firedrake` function space + (of type :class:`~firedrake.functionspaceimpl.WithGeometry`) on + which to convert to/from. Must be a 'Lagrange' or + 'Discontinuous Lagrange' space. + :arg warn_if_cg: If *True*, warn if the space is "CG" + """ + # validate function space + if not isinstance(function_space, WithGeometry): + raise TypeError("function_space must be of type WithGeometry, not " + "%s" % type(function_space)) + + family = function_space.ufl_element().family() + acceptable_families = ['Discontinuous Lagrange', 'Lagrange'] + if family not in acceptable_families: + raise ValueError("function_space.ufl_element().family() must be " + "one of %s, not '%s'" % + (acceptable_families, family)) + if family == 'Lagrange' and warn_if_cg: + warn("Functions in continuous function space will be projected " + "to/from a 'Discontinuous Lagrange' space. Make sure " + "any operators evaluated are continuous. " + "Pass warn_if_cg=False to suppress this warning.") + + # build DG space if necessary + family = function_space.ufl_element().family() + if family == 'Discontinuous Lagrange': + dg_function_space = function_space + elif family == 'Lagrange': + mesh = function_space.mesh() + degree = function_space.ufl_element().degree() + shape = function_space.shape + if shape is None or len(shape) == 0: + dg_function_space = FunctionSpace(mesh, "DG", degree) + elif len(shape) == 1: + dg_function_space = VectorFunctionSpace(mesh, "DG", degree, + dim=shape) + else: + dg_function_space = TensorFunctionSpace(mesh, "DG", degree, + shape=shape) + else: + acceptable_families = ['Discontinuous Lagrange', 'Lagrange'] + raise ValueError("function_space.ufl_element().family() must be " + "one of %s, not '%s'" % + (acceptable_families, family)) + + # store function space and dg function space + self.function_space = function_space + self.dg_function_space = dg_function_space + self.is_dg = function_space == dg_function_space + + def from_firedrake(self, density): + """ + Convert the density into a form acceptable by an bound operation + in an external library + + :arg density: A :class:`~firedrake.function.Function` holding the + density. + :returns: The converted density + """ + raise NotImplementedError + + def to_firedrake(self, evaluated_potential, out=None): + """ + Convert the evaluated potential from an external library + into a firedrake function + + :arg evaluated_potential: the evaluated potential + :arg out: If not *None*, store the converted potential into this + :class:`firedrake.function.Function`. + :return: *out* if it is not *None*, otherwise a new + :class:`firedrake.function.Function` storing the evaluated + potential + """ + raise NotImplementedError + + +class MeshmodeConnection(PotentialEvaluationLibraryConnection): + """ + Build a connection into :mod:`meshmode` + """ + def __init__(self, actx, function_space, warn_if_cg=True, + meshmode_connection_kwargs=None): + """ + For other args see :class:`PotentialEvaluationLibraryConnection` + + :arg actx: a :class:`meshmode.array_context.PyOpenCLArrayContext` + used for :func:`meshmode.interop.build_connection_from_firedrake` + and for conversion by the + :class:`meshmode.interop.FiredrakeConnection`. + :arg meshmode_connection_kwargs: A dict passes as kwargs + to :func:`meshmode.interop.build_connection_from_firedrake` + """ + PotentialEvaluationLibraryConnection.__init__(self, function_space, + warn_if_cg=warn_if_cg) + # validate actx + from meshmode.array_context import PyOpenCLArrayContext + if not isinstance(actx, PyOpenCLArrayContext): + raise TypeError("actx must be of type PyOpenCLArrayContext, not " + "%s." % type(actx)) + self.actx = actx + # validate meshmode_connection_kwargs + if meshmode_connection_kwargs is None: + meshmode_connection_kwargs = {} + if not isinstance(meshmode_connection_kwargs, dict): + raise TypeError("meshmode_connection_kwargs must be *None* or of " + "type dict, not '%s'." % + type(meshmode_connection_kwargs)) + # build a meshmode connection + from meshmode.interop.firedrake import build_connection_from_firedrake + mm_conn = build_connection_from_firedrake(actx, + self.dg_function_space, + **meshmode_connection_kwargs) + self.meshmode_connection = mm_conn + + def from_firedrake(self, density): + return self.meshmode_connection.from_firedrake(density, actx=self.actx) + + def to_firedrake(self, evaluated_potential, out=None): + return self.meshmode_connection.from_meshmode(evaluated_potential, + out=out) + + +class VolumentialConnection(PotentialEvaluationLibraryConnection): + """ + Build a connection into :mod:`volumential` + """ + def __init__(self, cl_ctx, queue, + function_space, + box_fmm_geometry_kwargs, + reorder_potentials, + meshmode_connection_kwargs, + warn_if_cg=True): + """ + For other args see :class:`PotentialEvaluationLibraryConnection` + + :arg cl_ctx: A :class:`pyopencl.Context` + :arg queue: A :class:`pyopencl.CommandQueue` + :arg box_fmm_geometry_kwargs: kwargs passed to + :class:`volumential.geometry.BoxFMMGeometryFactory`. + Must include the key *'nlevels'* + :arg reorder_potentials: *True* iff the volumential FMM + is reordering potentials. *False* otherwise. + :arg meshmode_connection_kwargs: A dict passes as kwargs + to :func:`meshmode.interop.build_connection_from_firedrake`. + :mod:`meshmode` is used an intermediate level between + :mod:`firedrake` and :mod:`volumential`. + """ + # super + PotentialEvaluationLibraryConnection.__init__(self, function_space, + warn_if_cg=warn_if_cg) + # validate cl_ctx argument + from pyopencl import Context + if not isinstance(cl_ctx, Context): + raise TypeError("cl_ctx must be of type Context, not " + "%s." % type(cl_ctx)) + + # validate queue argument + from pyopencl import CommandQueue + if not isinstance(queue, CommandQueue): + raise TypeError("queue must be of type CommandQueue, not " + "%s." % type(queue)) + + # validate reorder_potentials + if not isinstance(reorder_potentials, bool): + raise TypeError("reorder_potentials must be of type bool, not " + "%s." % type(reorder_potentials)) + + # validate box_fmm_geometry_kwargs + if not isinstance(box_fmm_geometry_kwargs, dict): + raise TypeError("box_fmm_geometry_kwargs must be of " + "type dict, not '%s'." % + type(box_fmm_geometry_kwargs)) + if 'nlevels' not in box_fmm_geometry_kwargs: + raise ValueError("box_fmm_geometry_kwargs missing required keyword" + " 'nlevels'.") + + # validate meshmode_connection_kwargs + if meshmode_connection_kwargs is None: + meshmode_connection_kwargs = {} + if not isinstance(meshmode_connection_kwargs, dict): + raise TypeError("meshmode_connection_kwargs must be *None* or of " + "type dict, not '%s'." % + type(meshmode_connection_kwargs)) + + # Build intermediate connection into meshmode + from meshmode.interop.firedrake import build_connection_from_firedrake + from meshmode.array_context import PyOpenCLArrayContext + actx = PyOpenCLArrayContext(queue) + meshmode_connection = build_connection_from_firedrake( + actx, function_space, **meshmode_connection_kwargs) + + # Build connection from meshmode into volumential + # (following https://gitlab.tiker.net/xywei/volumential/-/blob/fe2c3e7af355d5c527060e783237c124c95397b5/test/test_interpolation.py#L72 ) # noqa : E501 + from volumential.geometry import ( + BoundingBoxFactory, BoxFMMGeometryFactory) + from volumential.interpolation import ElementsToSourcesLookupBuilder + dim = function_space.mesh().geometric_dimension() + bboxfmm_fac = BoxFMMGeometryFactory(cl_ctx, dim=dim, + **box_fmm_geometry_kwargs) + boxgeo = bboxfmm_fac(queue) + elt_to_src_lookup_fac = ElementsToSourcesLookupBuilder( + cl_ctx, tree=boxgeo.tree, discr=meshmode_connection.discr) + elt_to_src_lookup, evt = elt_to_src_lookup_fac(queue) + + # Build connection from volumential into meshmode + from volumential.interpolation import LeavesToNodesLookupBuilder + leaves_to_node_lookup_fac = LeavesToNodesLookupBuilder( + cl_ctx, trav=boxgeo.trav, discr=meshmode_connection.discr) + leaves_to_node_lookup, evt = leaves_to_node_lookup_fac(queue) + + # Figure out whether or not conversion from volumential -> meshmode + # should user order or tree order + order = "user" if reorder_potentials else "tree" + + # Store maps for conversion to/from volumential + self.meshmode_connection = meshmode_connection + self.to_volumential_lookup = elt_to_src_lookup + self.from_volumential_lookup = leaves_to_node_lookup + self.from_volumential_order = order + # store necessary computing contextx + self.queue = queue + self.actx = actx + + def from_firedrake(self, density): + # pass operand into a meshmode DOFArray + from meshmode.dof_array import flatten + meshmode_density = \ + self.meshmode_connection.from_firedrake(density, actx=self.actx) + meshmode_density = flatten(meshmode_density) + + # pass flattened operand into volumential interpolator + from volumential.interpolation import interpolate_from_meshmode + volumential_density = \ + interpolate_from_meshmode(self.queue, + meshmode_density, + self.to_volumential_lookup, + order="user") # user order is more intuitive + return volumential_density + + def to_firedrake(self, evaluated_potential, out=None): + # pass volumential back to meshmode DOFArray + from volumential.interpolation import interpolate_to_meshmode + from meshmode.dof_array import unflatten + meshmode_pot_vals = \ + interpolate_to_meshmode(self.queue, + evaluated_potential, + self.from_volumential_lookup, + order=self.from_volumential_order) + meshmode_pot_vals = unflatten(self.actx, + self.meshmode_connection.discr, + meshmode_pot_vals) + # get meshmode data back as firedrake function + return self.meshmode_connection.from_meshmode(meshmode_pot_vals, + out=self) + + +""" +functions + PytentialOperation + SingleLayerPotential + DoubleLayerPotential + VolumentialOperation + VolumePotential +""" + + +class PotentialSourceAndTarget: + """ + Holds the source and target for a layer or volume potential + """ + def __init__(self, mesh, + source_region_dim=None, + source_region_id=None, + target_region_dim=None, + target_region_id=None): + """ + Source and target of a layer or volume potential + + Region dims must be the geometric dimension or the geometric + dimension minus 1. *None* indicates the geometric dimension. + + Region ids must be either a valid mesh subdomain id (an + integer or tuple) or *None*. *None* indicates either + the entire mesh, or the entire exterior boundary + (as determined by the value of region). + + Currently, if both the source and target are of dimension + geometric dimension - 1, they must be disjoint. + """ + pass + + +def PytentialOperation(density, unbound_op, potential_src_and_tgt, + function_space=None): + pass + + +def VolumentialOperation(density, unbound_op, potential_src_and_tgt, + function_space=None): + pass + + +def SingleLayerPotential(density, kernel, potential_src_and_tgt, + function_space=None, + operator_kwargs=None): + pass + + +def DoubleLayerPotential(density, kernel, potential_src_and_tgt, + function_space=None, + operator_kwargs=None): + pass + + class PytentialLayerOperation(AbstractExternalOperator): r""" Evaluates a pytential layer operation on a 2 or 3D mesh From 83c1ac6a1b36147f99a86c6102bb0cb93cb4f332 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Thu, 18 Feb 2021 14:26:02 -0600 Subject: [PATCH 02/20] finished connections and pytential generic op --- firedrake/layer_potentials.py | 951 +++++++++++++--------------------- 1 file changed, 373 insertions(+), 578 deletions(-) diff --git a/firedrake/layer_potentials.py b/firedrake/layer_potentials.py index b7f1ac496a..43ac1e0b41 100644 --- a/firedrake/layer_potentials.py +++ b/firedrake/layer_potentials.py @@ -1,8 +1,10 @@ import numpy as np +import firedrake.function + +from firedrake.functionspaceimpl import WithGeometry +from firedrake.mesh import MeshGeometry from firedrake.pointwise_operators import AbstractExternalOperator from firedrake.utils import cached_property -from firedrake.functionspaceimpl import WithGeometry -import firedrake.function from firedrake import Function, FunctionSpace, interpolate, Interpolator, \ SpatialCoordinate, VectorFunctionSpace, TensorFunctionSpace, project @@ -181,7 +183,8 @@ class PotentialEvaluationLibraryConnection: """ A connection to an external library for potential evaluation """ - def __init__(self, function_space, warn_if_cg=True): + def __init__(self, function_space, potential_source_and_target, + warn_if_cg=True): """ Initialize self to work on the function space @@ -189,6 +192,8 @@ def __init__(self, function_space, warn_if_cg=True): (of type :class:`~firedrake.functionspaceimpl.WithGeometry`) on which to convert to/from. Must be a 'Lagrange' or 'Discontinuous Lagrange' space. + :arg potential_source_and_target: A :class:`PotentialSourceAndTarget`. + mesh must match that of function_space. :arg warn_if_cg: If *True*, warn if the space is "CG" """ # validate function space @@ -208,6 +213,19 @@ def __init__(self, function_space, warn_if_cg=True): "any operators evaluated are continuous. " "Pass warn_if_cg=False to suppress this warning.") + # validate potential_source_and_targets + if not isinstance(potential_source_and_target, + PotentialSourceAndTarget): + raise TypeError("potential_source_and_targets must be of type " + "PotentialSourceAndTarget, not '%s'." + % type(potential_source_and_target)) + + # make sure meshes match + if potential_source_and_target.mesh is not function_space.mesh(): + raise ValueError("function_space.mesh() and " + "potential_source_and_target.mesh must be the same" + " obejct.") + # build DG space if necessary family = function_space.ufl_element().family() if family == 'Discontinuous Lagrange': @@ -234,6 +252,8 @@ def __init__(self, function_space, warn_if_cg=True): self.function_space = function_space self.dg_function_space = dg_function_space self.is_dg = function_space == dg_function_space + # store source and targets + self.source_and_target = potential_source_and_target def from_firedrake(self, density): """ @@ -265,7 +285,8 @@ class MeshmodeConnection(PotentialEvaluationLibraryConnection): """ Build a connection into :mod:`meshmode` """ - def __init__(self, actx, function_space, warn_if_cg=True, + def __init__(self, function_space, potential_source_and_target, actx, + warn_if_cg=True, meshmode_connection_kwargs=None): """ For other args see :class:`PotentialEvaluationLibraryConnection` @@ -274,17 +295,31 @@ def __init__(self, actx, function_space, warn_if_cg=True, used for :func:`meshmode.interop.build_connection_from_firedrake` and for conversion by the :class:`meshmode.interop.FiredrakeConnection`. - :arg meshmode_connection_kwargs: A dict passes as kwargs - to :func:`meshmode.interop.build_connection_from_firedrake` + :arg meshmode_connection_kwargs: A dict passed as kwargs + to :func:`meshmode.interop.build_connection_from_firedrake`. + Must not include 'restrict_to_boundary' """ - PotentialEvaluationLibraryConnection.__init__(self, function_space, - warn_if_cg=warn_if_cg) + PotentialEvaluationLibraryConnection.__init__( + self, + function_space, + potential_source_and_target, + warn_if_cg=warn_if_cg) + + # FIXME : allow subdomain regions + tdim = potential_source_and_target.mesh.topological_dimension() + if potential_source_and_target.get_source_dim() != tdim: + if potential_source_and_target.get_source_id() != "everywhere": + raise NotImplementedError("subdomain sources not implemented") + if potential_source_and_target.get_target_dim() != tdim: + if potential_source_and_target.get_target_id() != "everywhere": + raise NotImplementedError("subdomain targets not implemented") + # validate actx from meshmode.array_context import PyOpenCLArrayContext if not isinstance(actx, PyOpenCLArrayContext): raise TypeError("actx must be of type PyOpenCLArrayContext, not " "%s." % type(actx)) - self.actx = actx + # validate meshmode_connection_kwargs if meshmode_connection_kwargs is None: meshmode_connection_kwargs = {} @@ -292,33 +327,135 @@ def __init__(self, actx, function_space, warn_if_cg=True, raise TypeError("meshmode_connection_kwargs must be *None* or of " "type dict, not '%s'." % type(meshmode_connection_kwargs)) - # build a meshmode connection + if 'restrict_to_boundary' in meshmode_connection_kwargs: + raise ValueError("'restrict_to_boundary' must not be set by " + "meshmode_connection_kwargs") + + # build a meshmode connection for the source + src_bdy = None + if potential_source_and_target.get_source_dim() != tdim: + # sanity check + assert potential_source_and_target.get_source_dim() == tdim - 1 + src_bdy = potential_source_and_target.get_source_id() + from meshmode.interop.firedrake import build_connection_from_firedrake - mm_conn = build_connection_from_firedrake(actx, - self.dg_function_space, - **meshmode_connection_kwargs) - self.meshmode_connection = mm_conn + src_conn = build_connection_from_firedrake(actx, + self.dg_function_space, + restrict_to_boundary=src_bdy, + **meshmode_connection_kwargs) + # If source is a boundary, build a connection to restrict it to the + # boundary + restrict_conn = None + if src_bdy is not None: + # convert "everywhere" to meshmode BTAG_ALL + meshmode_src_bdy = src_bdy + if meshmode_src_bdy == "everywhere": + from meshmode.mesh import BTAG_ALL + meshmode_src_bdy = BTAG_ALL + # get group factory + order = self.dg_function_space.degree() + from meshmode.discretization.poly_element import \ + PolynomialRecursiveNodesGroupFactory + default_grp_factory = \ + PolynomialRecursiveNodesGroupFactory(order, 'lgl') + grp_factory = meshmode_connection_kwargs.get('grp_factory', + default_grp_factory) + from meshmode.discretization.connection import make_face_restriction + restrict_conn = make_face_restriction(actx, + src_conn.discr, + grp_factory, + meshmode_src_bdy) + + # Build a meshmode connection for the target + tgt_bdy = None + if potential_source_and_target.get_target_dim() != tdim: + # sanity check + assert potential_source_and_target.get_target_dim() == tdim - 1 + tgt_bdy = potential_source_and_target.get_target_id() + + # Can we re-use the source connection? + src_dim = potential_source_and_target.get_source_dim() + tgt_dim = potential_source_and_target.get_target_dim() + if tgt_bdy == src_bdy and src_dim == tgt_dim: + tgt_conn = src_conn + else: + # If not, build a new one + tgt_conn = build_connection_from_firedrake( + actx, + self.dg_function_space, + restrict_to_boundary=tgt_bdy, + **meshmode_connection_kwargs) + + # store computing context + self.actx = actx + # store connections + self.source_to_meshmode_connection = src_conn + self.restrict_source_to_boundary = restrict_conn + self.target_to_meshmode_connection = tgt_conn def from_firedrake(self, density): - return self.meshmode_connection.from_firedrake(density, actx=self.actx) + # make sure we are in the dg function space + if not self.is_dg: + density = project(density, self.dg_function_space) + # Convert to meshmode. + density = \ + self.source_to_meshmode_connection.from_firedrake(density, + actx=self.actx) + # restrict to boundary if necessary + if self.restrict_source_to_boundary is not None: + density = self.restrict_source_to_boundary(density) + return density def to_firedrake(self, evaluated_potential, out=None): - return self.meshmode_connection.from_meshmode(evaluated_potential, - out=out) + # if we are dg, it's simple + if self.is_dg: + return self.target_to_meshmode_connection.from_meshmode( + evaluated_potential, + out=out) + else: + # Otherwise, we have to project back to our function space + pot = \ + self.target_to_meshmode_connection.from_meshmode(evaluated_potential) + pot = project(evaluated_potential, self.function_space) + if out is not None: + out.dat.data[:] = pot.dat.data[:] + pot = out + return pot + + def get_source_discretization(self): + """ + Get the :class:`meshmode.discretization.Discretization` + of the source + """ + if self.restrict_source_to_boundary is None: + return self.source_to_meshmode_connection.discr + return self.restrict_source_to_boundary.to_discr + + def get_target_discretization(self): + """ + Get the :class:`meshmode.discretization.Discretization` + of the target + """ + return self.target_to_meshmode_connection.discr class VolumentialConnection(PotentialEvaluationLibraryConnection): """ Build a connection into :mod:`volumential` """ - def __init__(self, cl_ctx, queue, + def __init__(self, + potential_source_and_target, function_space, + cl_ctx, + queue, box_fmm_geometry_kwargs, reorder_potentials, meshmode_connection_kwargs, warn_if_cg=True): """ - For other args see :class:`PotentialEvaluationLibraryConnection` + For other args see :class:`PotentialEvaluationLibraryConnection`. + The *potential_source_and_target* must have a 3D + source, 3D target, and 3D mesh of co-dimension 0. :arg cl_ctx: A :class:`pyopencl.Context` :arg queue: A :class:`pyopencl.CommandQueue` @@ -333,8 +470,28 @@ def __init__(self, cl_ctx, queue, :mod:`firedrake` and :mod:`volumential`. """ # super - PotentialEvaluationLibraryConnection.__init__(self, function_space, - warn_if_cg=warn_if_cg) + PotentialEvaluationLibraryConnection.__init__( + self, + function_space, + potential_source_and_target, + warn_if_cg=warn_if_cg) + + # validate that sources, target, and mesh are of right dimension + if potential_source_and_target.get_source_dimension() != 3: + raise ValueError("potential source must have dimension 3") + if potential_source_and_target.get_target_dimension() != 3: + raise ValueError("potential target must have dimension 3") + if potential_source_and_target.mesh().geometric_dimension() != 3: + raise ValueError("mesh must have geometric dimension 3") + if potential_source_and_target.mesh().topological_dimension() != 3: + raise ValueError("mesh must have topological dimension 3") + + # FIXME : allow subdomain regions + if potential_source_and_target.get_source_id() != "everywhere": + raise NotImplementedError("subdomain sources not implemented") + if potential_source_and_target.get_target_id() != "everywhere": + raise NotImplementedError("subdomain targets not implemented") + # validate cl_ctx argument from pyopencl import Context if not isinstance(cl_ctx, Context): @@ -373,8 +530,12 @@ def __init__(self, cl_ctx, queue, from meshmode.interop.firedrake import build_connection_from_firedrake from meshmode.array_context import PyOpenCLArrayContext actx = PyOpenCLArrayContext(queue) - meshmode_connection = build_connection_from_firedrake( - actx, function_space, **meshmode_connection_kwargs) + meshmode_connection = MeshmodeConnection( + function_space, + potential_source_and_target, + actx, + warn_if_cg=warn_if_cg, + meshmode_connection_kwargs=meshmode_connection_kwargs) # Build connection from meshmode into volumential # (following https://gitlab.tiker.net/xywei/volumential/-/blob/fe2c3e7af355d5c527060e783237c124c95397b5/test/test_interpolation.py#L72 ) # noqa : E501 @@ -411,8 +572,7 @@ def __init__(self, cl_ctx, queue, def from_firedrake(self, density): # pass operand into a meshmode DOFArray from meshmode.dof_array import flatten - meshmode_density = \ - self.meshmode_connection.from_firedrake(density, actx=self.actx) + meshmode_density = self.meshmode_connection.from_firedrake(density) meshmode_density = flatten(meshmode_density) # pass flattened operand into volumential interpolator @@ -437,13 +597,11 @@ def to_firedrake(self, evaluated_potential, out=None): self.meshmode_connection.discr, meshmode_pot_vals) # get meshmode data back as firedrake function - return self.meshmode_connection.from_meshmode(meshmode_pot_vals, - out=self) + return self.meshmode_connection.to_firedrake(meshmode_pot_vals, out=self) """ functions - PytentialOperation SingleLayerPotential DoubleLayerPotential VolumentialOperation @@ -461,577 +619,214 @@ def __init__(self, mesh, target_region_dim=None, target_region_id=None): """ - Source and target of a layer or volume potential + Source and target of a layer or volume potential. + The mesh must have co-dimension 0 or 1 - Region dims must be the geometric dimension or the geometric - dimension minus 1. *None* indicates the geometric dimension. + By region_dim, we mean the topological dimension of the + region. + Regions must have co-dimensions 0 or 1. + They cannot have a higher dimension than the mesh. + *None* indicates the topological dimension of the + mesh. Region ids must be either a valid mesh subdomain id (an integer or tuple) or *None*. *None* indicates either the entire mesh, or the entire exterior boundary (as determined by the value of region). + The string "everywhere" is also equivalent to a + value of *None* + + NOTE: For implementation reasons, if the target is + a boundary, instead of just evaluating the + target at DOFs on that boundary, it is instead + evaluated at any DOF of a cell which has at least + one vertex on the boundary. + """ + # validate mesh + if not isinstance(mesh, MeshGeometry): + raise TypeError("mesh must be of type MeshGeometry, not '%s'." % + type(mesh)) + + # mesh must have co-dimension 1 or 0. + valid_codims = [0, 1] + codim = mesh.geometric_dimension() - mesh.topological_dimension() + if codim not in valid_codims: + raise ValueError("mesh has invalid co-dimension of %s. " + "co-dimension must be one of %s." % + (codim, valid_codims)) + + # validate dims + for dim in [source_region_dim, target_region_dim]: + if dim is None: + continue + if not isinstance(dim, int): + raise TypeError("source and target dimensions must be *None*" + " or an integer, not of type '%s'." % type(dim)) + valid_dims = set([mesh.topological_dimension(), + mesh.geometric_dimension()-1]) + if dim < mesh.geometric_dimension() - 1: + raise ValueError("source and target dimensions must be " + "one of %s or *None*, not '%s'." % + (valid_dims, dim)) + # Get dims if *None* + if source_region_dim is None: + source_region_dim = mesh.topological_dimension() + if target_region_dim is None: + target_region_dim = mesh.topological_dimension() + + # Validate region ids type and validity + for id_, dim in zip([source_region_id, target_region_id], + [source_region_dim, target_region_dim]): + if id_ is None or id_ == "everywhere": + continue + # standardize id_ to a tuple + if isinstance(id_, int): + id_ = tuple(id_) + # check type + if not isinstance(id_, tuple): + raise TypeError("source and target region ids must be " + "*None*, an int, or a tuple, not of type") + # boundary case: + if dim == mesh.topological_dimension() - 1: + if not set(id_) <= set(mesh.exterior_facets.unique_markers): + raise ValueError(("boundary region ids %s are not a " + + "subset of mesh.exterior_facets." + + "unique_markers = %s.") % (id_, + mesh.exterior_facets.unique_markers)) + else: + # sanity check + assert dim == mesh.topological_dimension() + # this caches the cell subset on the mesh, which we'll + # probably be doing anyway if we're targeting the region + # + # It also throws a ValueError if the id_ is invalid + mesh.cell_subset(id_) + + # handle None case + if source_region_id is None: + source_region_id = "everywhere" + if target_region_id is None: + target_region_id = "everywhere" + + # store mesh, region dims, and region ids + self.mesh = mesh + self._source_region_dim = source_region_dim + self._target_region_dim = source_region_dim + self._source_region_id = source_region_id + self._target_region_id = source_region_id + + def get_source_dimension(self): + """ + Get the topological dimension of the source region + """ + return self._source_region_dim - Currently, if both the source and target are of dimension - geometric dimension - 1, they must be disjoint. + def get_target_dimension(self): """ - pass + Get the topological dimension of the target region + """ + return self._target_region_dim + def get_source_id(self): + """ + Get the subdomain id of the source, or the string + "everywhere" to represent the entire exterior boundary/mesh + """ + return self._source_region_id -def PytentialOperation(density, unbound_op, potential_src_and_tgt, - function_space=None): - pass + def get_target_id(self): + """ + Get the subdomain id of the target, or the string + "everywhere" to represent the entire exterior boundary/mesh + """ + return self._target_region_id -def VolumentialOperation(density, unbound_op, potential_src_and_tgt, - function_space=None): +def PytentialOperation(actx, + density, + unbound_op, + density_name, + potential_src_and_tgt, + **kwargs): + """ + :arg actx: A :class:`meshmode.dof_array.PyOpenCLArrayContext` + :arg density: the :mod:`firedrake` density function + :arg unbound_op: A :mod:`pytential` operation which has not been + bound (e.g. a :mod:`pymbolic` expression) + :arg density_name: A string, the name of the density function + in the unbound_op + :arg potential_src_and_tgt: A :class:`PotentialSourceAndTarget` + + :kwarg warn_if_cg: Pass *False* to suppress the warning if + a "CG" space is used + :kwarg meshmode_connection_kwargs: *None*, or + a dict with arguments to pass to + :func:`meshmode.interop.firedrake.build_connection_from_firedrake` + :kwarg qbx_kwargs: *None*, or a dict with arguments to pass to + :class:`pytential.qbx.QBXLayerPotentialSource`. + :kwarg op_kwargs: kwargs passed to the invocation of the + bound pytential operator (e.g. {'k'=0.5} for + a :class:`sumpy.kernel.HelmholtzKernel`). + Must not include *density_name* + """ + # make sure density name is a string + if not isinstance(density_name, str): + raise TypeError("density_name must be of type str, not '%s'." % + density_name) + + # get kwargs to build connection + warn_if_cg = kwargs.get('warn_if_cg', True) + meshmode_connection_kwargs = kwargs.get('meshmode_connection_kwargs', None) + + # Build meshmode connection + meshmode_connection = MeshmodeConnection( + density.function_space(), + potential_src_and_tgt, + actx, + warn_if_cg=warn_if_cg, + meshmode_connection_kwargs=meshmode_connection_kwargs) + + # Build QBX + src_discr = meshmode_connection.get_source_discretization() + qbx_kwargs = kwargs.get('qbx_kwargs', None) + from pytential.qbx import QBXLayerPotentialSource + qbx = QBXLayerPotentialSource(src_discr, **qbx_kwargs) + + # Bind pytential operator + tgt_discr = meshmode_connection.get_target_discretization() + from pytential import bind + pyt_op = bind((qbx, tgt_discr), unbound_op) + + # Get operator kwargs + op_kwargs = kwargs.get('op_kwargs', {}) + if density_name in op_kwargs: + raise ValueError(f"density_name '{density_name}' should not be included" + " in op_kwargs.") + + # build bound operator that just takes density as argument + def bound_op_with_kwargs(density_arg): + op_kwargs[density_name] = density_arg + return pyt_op(**op_kwargs) + + # Now build and return Potential object + return Potential(density, + connection=meshmode_connection, + potential_operator=bound_op_with_kwargs) + + +def SingleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): pass -def SingleLayerPotential(density, kernel, potential_src_and_tgt, - function_space=None, - operator_kwargs=None): +def DoubleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): pass -def DoubleLayerPotential(density, kernel, potential_src_and_tgt, - function_space=None, - operator_kwargs=None): +def VolumentialOperation(density, unbound_op, potential_src_and_tgt, + function_space=None): pass -class PytentialLayerOperation(AbstractExternalOperator): - r""" - Evaluates a pytential layer operation on a 2 or 3D mesh - (with geometric dim == topological dim) which has - a source and target of co-dimension 1 (the source - and target must be external boundaries in the mesh). - - IMPORTANT: The pytential normal direction is OPPOSITE that - of the firedrake normal direction. - This is not automatically accounted for. - (Although in :class:`SingleLayerPotential` - and :class:`DoubleLayerPotential` it is automatically - accounted for) - - Note: The function_space is the function space of the target, - this must be a "CG" space or "DG" space. - The function space of the operand (the density) - MUST be a "DG" space. - - :kwarg operator_data: A map with required keys - - * 'op': The pytential operation (see, for example, - :class:`SingleLayerPotential` or :class:`DoubleLayerPotential`). - Must be made up of :class:`pytential.sym.primitives` - as described in the *expr* argument to - :function:`pytential.bind`. This is not validated. - The output must have the same shape as the - function_space. - * 'density_name': The name of the density function in the - pytential operation *op* (see, for example, - :class:`SingleLayerPotential` or - :class:`DoubleLayerPotential`). - * 'actx': a :mod:`meshmode` :class:`meshmode.array_context.ArrayContext` - * 'source_bdy_id' A boundary id representing the source - * 'target_bdy_id' A boundary id representing the target - - And optional keys - - * 'op_kwargs': (optional) keyword arguments to be passed to the - evaluator of the pytential operation. For instance, - if a :class:`sumpy.kernel.HelmholtzKernel` is being - evaluated with a kappa of 0.5, you might pass - `operator_data['op_kwargs'] = {'k': 0.5}`. - DO NOT include the density functions in these kwargs, - it will be automatically inserted. - * 'project_to_dg': (optional) If *True*, when an operand in "CG" - space is received, it is interpolated into - a "DG" space of the same degree. - Default *False* - * 'grp_factory': (optional) An interpolatory group factory - inheriting from :class:`meshmode.discretization.ElementGroupFactory` - to be used in the intermediate :mod:`meshmode` representation - * 'qbx_order': As described in :class:`pytential.qbx.QBXLayerPotentialSource`. - Clear of FMM error and clear of quadrature area, should - see convergence order qbx_order + 1 (in layer approximation). - Default is function space degree + 2 - * 'fine_order': As described in :class:`pytential.qbx.QBXLayerPotentialSource` - Has to do with refinements used by QBX. - Default is 4 * function space degree - * 'fmm_order': As described in :class:`pytential.qbx.QBXLayerPotentialSource` - FMM error is bounded by c^{fmm_order + 1}, where - c = 0.5 for 2D and 0.75 for 3D. Default value is 6. - """ - - _external_operator_type = 'GLOBAL' - - def __init__(self, *operands, **kwargs): - self.operands = operands - self.kwargs = kwargs - - AbstractExternalOperator.__init__(self, operands[0], **kwargs) - - # Validate input - # print(len(self.operands), self.derivatives) - # assert self.derivatives == (0,), \ - # "Derivatives of single layer potential not currently supported: " + str(self.derivatives) - # from firedrake import Function - # - # if not isinstance(operand, Function): - # raise TypeError(":arg:`operand` must be of type firedrake.Function" - # ", not %s." % type(operand)) - # if operand.function_space().shape != tuple(): - # raise ValueError(":arg:`operand` must be a function with shape ()," - # " not %s." % operand.function_space().shape) - operator_data = kwargs["operator_data"] - assert isinstance(operator_data, dict) - required_keys = ('op', 'density_name', - 'actx', 'source_bdy_id', 'target_bdy_id') - optional_keys = ('op_kwargs', 'project_to_dg', - 'grp_factory', 'qbx_order', 'fmm_order', 'fine_order',) - permissible_keys = required_keys + optional_keys - if not all(key in operator_data for key in required_keys): - raise ValueError("operator_data is missing one of the required " - "keys: %s" % required_keys) - if not all(key in permissible_keys for key in operator_data): - raise ValueError("operator_data contains an unexpected key. All " - "keys must be one of %s." % (permissible_keys,)) - - @cached_property - def _evaluator(self): - return PytentialLayerOperationEvaluator(self, - *self.operands, - **self.kwargs) - - def _evaluate(self): - 1/0 - - def _compute_derivatives(self): - 1/0 - - def _evaluate_action(self, *args): - return self._evaluator._evaluate_action(*args) - - -class SingleOrDoubleLayerPotential(PytentialLayerOperation): - r""" - This is an abstract class to avoid code duplication between - single and double layer potentials. One should only - instantiate a :class:`SingleLayerPotential` or - :class:`DoubleLayerPotential` - - A single layer potential evaluates to - - .. math:: - - f(x)|_{x\in\Gamma} = \int_\Om K(x-y) op(y) \,dx - - A double layer potential evaluates to - - .. math:: - - f(x)|_{x\in\Gamma} = \int_\Om \partial_n K(x-y) op(y) \,dx - - - where \Gamma is the target boundary id and \Om is the - source bdy id - as described in :class:`~firedrake.layer_potentials.LayerPotential`, - and K is a :class:`sumpy.kernel.Kernel`. - The function space must have scalar shape. - - :kwarg operator_data: A map as described in - :class:`~firedrake.layer_potentials.LayerPotential` - except that - - * 'kernel' must be included, and map to a value of type - :class:`sumpy.kernel.Kernel` - * (Optional) 'kernel_kwargs': A map which tells the names of required - arguments to the kernel. For example, - if the kernel is a HelmohltzKernel, - (which requires an argument 'k'), - and you are using a :class:`pytential.sym.var` - of name 'kappa', one would pass - {'k': 'kappa'}. - NOTE: If kappa = 0.5, the - corresponding operator_data['op_kwargs'] - argument would be {'kappa': 0.5} - * 'op' must not be included - """ - def __init__(self, *operands, **kwargs): - operator_data = kwargs["operator_data"] - # Make sure invalid keys are not present and kernel is present - if 'op' in operator_data.keys(): - raise ValueError("operator_data must not contain key 'op'") - if 'kernel' not in operator_data: - raise ValueError("Missing 'kernel' in operator_data") - # Get kernel and validate it - kernel = operator_data['kernel'] - from sumpy.kernel import Kernel - if not isinstance(kernel, Kernel): - raise TypeError("operator_data['kernel'] must be of type " - "sumpy.kernel.Kernel, not %s." % type(kernel)) - # Make sure have valid density name - if 'density_name' not in operator_data: - raise ValueError("Missing 'density_name' in operator_data") - density_name = operator_data['density_name'] - if not isinstance(density_name, str): - raise TypeError("operator_data['density_name'] must be of type str" - ", not '%s'." % type(density_name)) - # Get kernel kwargs if any - kernel_kwargs = operator_data.get('kernel_kwargs', {}) - # Build single-layer potential - from pytential.sym import var - op = self._getOp(kernel, var(density_name), kernel_kwargs) - del operator_data['kernel'] - del operator_data['kernel_kwargs'] - operator_data['op'] = op - # Finish initialization - super().__init__(*operands, **kwargs) - - def _getOp(self, kernel, density_name, kernel_kwargs): - raise NotImplementedError("Must instantiate a SingleLayerPotential" - " or DoubleLayerPotential." - " SingleOrDoubleLayerPotential is abstract.") - - -class SingleLayerPotential(SingleOrDoubleLayerPotential): - r""" - Layer potential which evaluates to - - .. math:: - - f(x)|_{x\in\Gamma} = \int_\Om K(x-y) op(y) \,dx - - As described in :class:`SingleOrDoubleLayerPotential` - """ - def _getOp(self, kernel, density_name, kernel_kwargs): - from pytential import sym - return sym.S(kernel, sym.var(density_name), **kernel_kwargs) - - -class DoubleLayerPotential(SingleOrDoubleLayerPotential): - r""" - Layer potential which evaluates to - - .. math:: - - f(x)|_{x\in\Gamma} = \int_\Om \partial_n K(x-y) op(y) \,dx - - - As described in :class:`SingleOrDoubleLayerPotential` - """ - def _getOp(self, kernel, density_name, kernel_kwargs): - # Build double-layer potential (pytential normal points opposite - # direction, so need *= -1) - from pytential import sym - pyt_inner_normal_sign = -1 - return pyt_inner_normal_sign * \ - sym.D(kernel, sym.var(density_name), **kernel_kwargs) - - -def _get_target_points_and_indices(fspace, boundary_ids): - """ - Get the points from the function space which lie on the given boundary - id as a pytential PointsTarget, and their indices into the - firedrake function - - :arg fspace: The function space - :arg boundary_ids: the boundary ids (an int or tuple of ints, - not validated) - :return: (target_indices, target_points) - """ - if isinstance(boundary_ids, int): - boundary_ids = tuple(boundary_ids) - target_indices = set() - for marker in boundary_ids: - # Make sure geometric since discontinuous space - target_indices |= set( - fspace.boundary_nodes(marker, 'geometric')) - target_indices = np.array(list(target_indices), dtype=np.int32) - - target_indices = np.array(target_indices, dtype=np.int32) - # Get coordinates of nodes - coords = SpatialCoordinate(fspace.mesh()) - function_space_dim = VectorFunctionSpace( - fspace.mesh(), - fspace.ufl_element().family(), - degree=fspace.ufl_element().degree()) - - coords = Function(function_space_dim).interpolate(coords) - coords = np.real(coords.dat.data) - - target_pts = coords[target_indices] - # change from [nnodes][ambient_dim] to [ambient_dim][nnodes] - target_pts = np.transpose(target_pts).copy() - from pytential.target import PointsTarget - return (target_indices, PointsTarget(target_pts)) - - -class PytentialLayerOperationEvaluator: - def __init__(self, pyt_layer, *operands, **kwargs): - self.pyt_layer = pyt_layer - if len(operands) != 1: - raise NotImplementedError("I haven't thought this through") - operand = operands[0] - - # FIXME : Allow ufl - # Check operand is a function - if not isinstance(operand, Function): - raise TypeError("operand %s must be of type Function, not %s" % - operand, type(operand)) - - # Get funciton space - function_space = kwargs["function_space"] - if function_space.ufl_element().family() not in ['Lagrange', - 'Discontinuous Lagrange']: - raise ValueError("function_space must have family 'Lagrange' or" - "'Discontinuous Lagrange', not '%s'" % - operand.function_space().ufl_element().family()) - # validate mesh - valid_geo_dims = [2, 3] - mesh = function_space.mesh() - if mesh.geometric_dimension() not in valid_geo_dims: - raise ValueError("function_space.mesh().geometric_dimension() " - "%s is not in %s" % - (mesh.geometric_dimension(), valid_geo_dims)) - if mesh.geometric_dimension() != mesh.topological_dimension(): - raise ValueError("function_space.mesh().topological_dimension() of " - "%s does not equal " - "function_space.mesh().geometric_dimension() of %s" - % (mesh.topolocial_dimension(), - mesh.geometric_dimension())) - - # Create an interpolator from the original operand so that - # we can re-interpolate it into the space at each evaluation - # in case it has changed (e.g. cost functional in PDE constrained optimization. - - # Get operator data - operator_data = kwargs["operator_data"] - # project to dg? - project_to_dg = operator_data.get('project_to_dg', False) - if not isinstance(project_to_dg, bool): - raise TypeError("operator_data['project_to_dg'] must be of type " - "bool, not %s." % type(project_to_dg)) - # get operand fspace - operand_fspace = operand.function_space() - if operand_fspace.ufl_element().family() != 'Discontinuous Lagrange': - # If operand is not "DG" and not projecting, throw error - if not project_to_dg: - raise ValueError("operand must live in a function space with family " - "Discontinuous Lagrange, not '%s'. Look at the " - "optional operator_data['project_to_dg'] argument." % - operand.function_space().ufl_element().family()) - # otherwise build new function space if necessary - elif operand_fspace.shape == (): - operand_fspace = FunctionSpace(operand_fspace.mesh(), - "DG", - degree=operand_fspace.ufl_element().degree()) - elif len(operand_fspace.shape) == 1: - operand_fspace = VectorFunctionSpace(operand_fspace.mesh(), - "DG", - degree=operand_fspace.ufl_element().degree(), - dim=operand_fspace.shape[0]) - else: - operand_fspace = TensorFunctionSpace(operand_fspace.mesh(), - "DG", - degree=operand_fspace.ufl_element().degree(), - shape=operand_fspace.shape) - else: - # No need to project - project_to_dg = False - - - # get op and op-shape - op = operator_data["op"] - # Get density-name and validate - density_name = operator_data["density_name"] - if not isinstance(density_name, str): - raise TypeError("operator_data['density_name'] must be of type " - " str, not '%s'." % type(density_name)) - # Get op kwargs and validate - op_kwargs = operator_data.get('op_kwargs', {}) - if not isinstance(op_kwargs, dict): - raise TypeError("operator_data['op_kwargs'] must be of type " - " dict, not '%s'." % type(op_kwargs)) - for k in op_kwargs.keys(): - if not isinstance(k, str): - raise TypeError("Key '%s' in operator_data['op_kwargs'] must " - " be of type str, not '%s'." % type(k)) - - # Validate actx type - actx = operator_data['actx'] - from meshmode.array_context import PyOpenCLArrayContext - if not isinstance(actx, PyOpenCLArrayContext): - raise TypeError("operator_data['actx'] must be of type " - "PyOpenCLArrayContext, not %s." % type(actx)) - - source_bdy_id = operator_data["source_bdy_id"] - target_bdy_id = operator_data["target_bdy_id"] - # Make sure bdy ids are appropriate types - if not isinstance(source_bdy_id, int): - raise TypeError("operator_data['source_bdy_id'] must be of type int," - " not type %s" % type(source_bdy_id)) - if isinstance(target_bdy_id, int): - target_bdy_id = tuple(target_bdy_id) - if not isinstance(target_bdy_id, tuple): - raise TypeError("operator_data['target_bdy_id'] must be an int " - " or a tuple of ints, not of type %s" % - type(target_bdy_id)) - for bdy_id in target_bdy_id: - if not isinstance(bdy_id, int): - raise TypeError("non-integer value '%s' found in " - "operator_data['target_bdy_id']" % bdy_id) - # Make sure bdy ids are actually boundary ids - valid_ids = function_space.mesh().exterior_facets.unique_markers - if not set(valid_ids) >= set(target_bdy_id): - raise ValueError("Invalid target bdy id(s): %s." % - set(target_bdy_id) - set(valid_ids)) - if source_bdy_id not in valid_ids: - raise ValueError("Invalid source bdy id: %s" % source_bdy_id) - # Make sure bdy ids are disjoint - if source_bdy_id in target_bdy_id: - raise NotImplementedError("source and target boundaries must be " - "disjoint") - - degree = function_space.ufl_element().degree() - # Get group factory, if any - grp_factory = operator_data.get('grp_factory', None) - # validate grp_factory - from meshmode.discretization.poly_element import ElementGroupFactory - if grp_factory is not None: - if not isinstance(grp_factory, ElementGroupFactory): - raise TypeError("operator_data['grp_factory'] must be *None*" - " or of type ElementGroupFactory, not %s." % - type(grp_factory)) - # Set defaults for qbx kwargs - qbx_order = kwargs.get('qbx_order', degree+2) - fine_order = kwargs.get('fine_order', 4 * degree) - fmm_order = kwargs.get('fmm_order', 6) - # Validate qbx kwargs - for var, name in zip([qbx_order, fine_order, fmm_order], - ['qbx_order', 'fine_order', 'fmm_order']): - if not isinstance(var, int): - raise TypeError("operator_data['%s'] must be of type int, " - "not %s." % (name, type(var))) - if not var > 0: - raise ValueError("operator_data['%s'] = %s is not positive." - % (name, var)) - - qbx_kwargs = {'qbx_order': qbx_order, - 'fine_order': fine_order, - 'fmm_order': fmm_order, - 'fmm_backend': 'fmmlib', - } - # }}} - - # Build connection into meshmode - from meshmode.interop.firedrake import build_connection_from_firedrake - meshmode_connection = build_connection_from_firedrake( - actx, operand_fspace, grp_factory=grp_factory, - restrict_to_boundary=source_bdy_id) - - # Build recursive nodes group factory if none provided - if grp_factory is None: - from meshmode.discretization.poly_element import \ - PolynomialRecursiveNodesGroupFactory - grp_factory = PolynomialRecursiveNodesGroupFactory(degree, family='lgl') - - # build connection meshmode near src boundary -> src boundary inside meshmode - from meshmode.discretization.connection import make_face_restriction - src_bdy_connection = make_face_restriction(actx, - meshmode_connection.discr, - grp_factory, - source_bdy_id) - - # Build QBX - from pytential.qbx import QBXLayerPotentialSource - qbx = QBXLayerPotentialSource(src_bdy_connection.to_discr, **qbx_kwargs) - # Get target, and store the firedrake indices of the points - target_indices, target = _get_target_points_and_indices(function_space, - target_bdy_id) - - # Bind pytential operator - from pytential import bind - self.pyt_op = bind((qbx, target), op) - self.density_name = density_name - self.op_kwargs = op_kwargs - - # Store attributes that we may need later - self.project_to_dg = project_to_dg - self.operand_fspace = operand_fspace - self.actx = actx - self.meshmode_connection = meshmode_connection - self.src_bdy_connection = src_bdy_connection - self.target_indices = target_indices - self.function_space = function_space - # so we don't have to keep making a fd function during conversion - # from meshmode. AbstractExternalOperator s aren't functions, - # so can't use *self* - self.fd_pot = Function(function_space) - # initialize to zero so that only has values on target boundary - self.fd_pot.dat.data[:] = 0.0 - - def _evaluate(self, *args): - operand, = self.pyt_layer.ufl_operands - # TODO : Fix this by allowing for ufl arguments - assert isinstance(operand, Function) - operand_discrete = operand - if self.project_to_dg: - operand_discrete = project(operand_discrete, self.operand_fspace) - - # pass operand into a meshmode DOFArray - from meshmode.dof_array import flatten - meshmode_src_vals = \ - self.meshmode_connection.from_firedrake(operand_discrete, actx=self.actx) - # pass operand onto boundary - meshmode_src_vals_on_bdy = self.src_bdy_connection(meshmode_src_vals) - - # Evaluate pytential potential - self.op_kwargs[self.density_name] = meshmode_src_vals_on_bdy - target_values = self.pyt_op(**self.op_kwargs) - - # FIXME : Is this correct? - fspace_shape = self.function_space.shape - if target_values.shape != fspace_shape and not \ - (fspace_shape == () and len(target_values.shape) == 1): - raise ValueError("function_space shape of %s does not match" - "operator shape %s." % - (fspace_shape, target_values.shape)) - # Make sure conversion back to firedrake works for non-scalar types - # (copied and modified from - # https://github.com/inducer/meshmode/blob/be1bcc9d395ca51d6903993eb57acc865acec243/meshmode/interop/firedrake/connection.py#L544-L555 # noqa - # ) - # If scalar, just reorder and resample out - if fspace_shape == (): - self.fd_pot.dat.data[self.target_indices] = target_values.get(queue=self.actx.queue)[:] - else: - # otherwise, have to grab each dofarray and the corresponding - # data from *function_data* - for multi_index in np.ndindex(fspace_shape): - # have to be careful to take view and not copy - index = (np.s_[:],) + multi_index - fd_data = self.fd_pot.dat.data[index] - dof_array = target_values[multi_index] - fd_data[self.target_indices] = dof_array.get(queue=self.actx.queue) - - # Store in vp - self.pyt_layer.dat.data[:] = self.fd_pot.dat.data[:] - # Return evaluated potential - return self.fd_pot - - def _evaluate_action(self, *args): - assert len(args) == 1 - action_coeffs = args[0] - assert len(action_coeffs) in [0, 1] - if len(action_coeffs) == 1: - action_coeff, = action_coeffs - ufl_operand, = self.pyt_layer.ufl_operands - # Store new action coefficients into operand - ufl_operand.dat.data[:] = action_coeff.dat.data[:] - # Evaluate operation on current operands - return self._evaluate() - - class VolumePotential(AbstractExternalOperator): r""" Evaluates to From 4761f0457a564084bd8c2e205ce171eea90fdd78 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Thu, 18 Feb 2021 16:00:39 -0600 Subject: [PATCH 03/20] Split up interface into multiple files --- firedrake/layer_potentials.py | 1226 ----------------- firedrake/potential_evaluation/__init__.py | 416 ++++++ firedrake/potential_evaluation/potentials.py | 212 +++ firedrake/potential_evaluation/pytential.py | 233 ++++ firedrake/potential_evaluation/volumential.py | 568 ++++++++ 5 files changed, 1429 insertions(+), 1226 deletions(-) delete mode 100644 firedrake/layer_potentials.py create mode 100644 firedrake/potential_evaluation/__init__.py create mode 100644 firedrake/potential_evaluation/potentials.py create mode 100644 firedrake/potential_evaluation/pytential.py create mode 100644 firedrake/potential_evaluation/volumential.py diff --git a/firedrake/layer_potentials.py b/firedrake/layer_potentials.py deleted file mode 100644 index 43ac1e0b41..0000000000 --- a/firedrake/layer_potentials.py +++ /dev/null @@ -1,1226 +0,0 @@ -import numpy as np -import firedrake.function - -from firedrake.functionspaceimpl import WithGeometry -from firedrake.mesh import MeshGeometry -from firedrake.pointwise_operators import AbstractExternalOperator -from firedrake.utils import cached_property - -from firedrake import Function, FunctionSpace, interpolate, Interpolator, \ - SpatialCoordinate, VectorFunctionSpace, TensorFunctionSpace, project - -from ufl.algorithms.apply_derivatives import VariableRuleset -from ufl.algorithms import extract_coefficients -from ufl.constantvalue import as_ufl -from ufl.core.multiindex import indices -from ufl.tensors import as_tensor - -from pyop2.datatypes import ScalarType -from warnings import warn - - -__all__ = ("DoubleLayerPotential", "PytentialLayerOperation", - "SingleLayerPotential", - "SingleOrDoubleLayerPotential", # included for the docs - "VolumePotential",) - - -class Potential(AbstractExternalOperator): - r""" - This is a function which represents a potential - computed on a firedrake function. - - For a firedrake function :math:`u:\Omega\to\mathbb C^n` - with :math:`\Omega \subseteq \mathbb R^m`, - the potential :math:`P(u):\Omega \to\mathbb C^n` is a - convolution of :math:`u` against some kernel function. - More concretely, given a source region :math:`\Gamma \subseteq \Omega` - a target region :math:`\Sigma \subseteq \Omega`, and - a kernel function :math:`K`, we define - - .. math:: - - P(u)(x) = \begin{cases} - \int_\Gamma K(x, y) u(y) \,dy & x \in \Sigma - 0 & x \notin \Sigma - \end{cases} - """ - _external_operator_type = 'GLOBAL' - - def __init__(self, density, **kwargs): - """ - :arg density: A :mod:`firedrake` :class:`firedrake.function.Function` - or UFL expression which represents the density :math:`u` - :kwarg connection: A :class:`PotentialEvaluationLibraryConnection` - :kwarg potential_operator: The external potential evaluation library - bound operator. - """ - # super - AbstractExternalOperator.__init__(self, density, **kwargs) - - # FIXME - for order in self.derivatives: - assert order == 1, "Assumes self.derivatives = (1,..,1)" - - # Get connection & bound op and validate - connection = kwargs.get("connection", None) - if connection is None: - raise ValueError("Missing kwarg 'connection'") - if not isinstance(connection, PotentialEvaluationLibraryConnection): - raise TypeError("connection must be of type " - "PotentialEvaluationLibraryConnection, not %s." - % type(connection)) - self.connection = connection - - self.potential_operator = kwargs.get("potential_operator", None) - if self.potential_operator is None: - raise ValueError("Missing kwarg 'potential_operator'") - - # Get function space and validate aginst bound op - function_space = self.ufl_function_space() - assert isinstance(function_space, WithGeometry) # sanity check - if function_space is not self.potential_operator.function_space(): - raise ValueError("function_space must be same object as " - "potential_operator.function_space().") - - # Make sure density is a member of our function space, if it is - if isinstance(density, Function): - if density.function_space() is not function_space: - raise ValueError("density.function_space() must be the " - "same as function_space") - # Make sure the shapes match, at least - elif density.shape != function_space.shape: - raise ValueError("Shape mismatch between function_space and " - "density. %s != %s." % - (density.shape, function_space.shape)) - - @cached_property - def _evaluator(self): - return PotentialEvaluator(self, - self.density, - self.connection, - self.potential_operator) - - def _evaluate(self): - raise NotImplementedError - - def _compute_derivatives(self, continuity_tolerance=None): - raise NotImplementedError - - def _evaluate_action(self, *args): - return self._evaluator._evaluate_action() - - def _evaluate_adjoint_action(self, *args): - return self._evaluator._evaluate_action() - - def evaluate_adj_component_control(self, x, idx): - derivatives = tuple(dj + int(idx == j) for j, dj in enumerate(self.derivatives)) - dN = self._ufl_expr_reconstruct_(*self.ufl_operands, derivatives=derivatives) - return dN._evaluate_adjoint_action((x.function,)).vector() - - def evaluate_adj_component(self, x, idx): - print(x, type(x)) - raise NotImplementedError - - -class PotentialEvaluator: - """ - Evaluates a potential - """ - def __init__(self, density, connection, potential_operator): - """ - :arg density: The UFL/firedrake density function - :arg connection: A :class:`PotentialEvaluationLibraryConnection` - :arg potential_operator: The external potential evaluation library - bound operator. - """ - self.density = density - self.connection = connection - self.potential_operator = potential_operator - - def _eval_potential_operator(self, action_coefficients, out=None): - """ - Evaluate the potential on the action coefficients and return. - If *out* is not *None*, stores the result in *out* - - :arg action_coefficients: A :class:`~firedrake.function.Function` - to evaluate the potential at - :arg out: If not *None*, then a :class:`~firedrake.function.Function` - in which to store the evaluated potential - :return: *out* if it is not *None*, otherwise a new - :class:`firedrake.function.Function` storing the evaluated - potential - """ - density = self.connection.from_firedrake(action_coefficients) - potential = self.potential_operator(density) - return self.connection.to_firedrake(potential, out=out) - - def _evaluate(self): - """ - Evaluate P(self.density) into self - """ - return self._eval_potential_operator(self.density, out=self) - - def _compute_derivatives(self): - """ - Return a function - Derivative(P, self.derivatives, self.density) - """ - # FIXME : Assumes derivatives are Jacobian - return self._eval_potential_operator - - def _evaluate_action(self, action_coefficients): - """ - Evaluate derivatives of layer potential at action coefficients. - i.e. Derivative(P, self.derivatives, self.density) evaluated at - the action coefficients and store into self - """ - operator = self._compute_derivatives() - return operator(action_coefficients, out=self) - - -class PotentialEvaluationLibraryConnection: - """ - A connection to an external library for potential evaluation - """ - def __init__(self, function_space, potential_source_and_target, - warn_if_cg=True): - """ - Initialize self to work on the function space - - :arg function_space: The :mod:`firedrake` function space - (of type :class:`~firedrake.functionspaceimpl.WithGeometry`) on - which to convert to/from. Must be a 'Lagrange' or - 'Discontinuous Lagrange' space. - :arg potential_source_and_target: A :class:`PotentialSourceAndTarget`. - mesh must match that of function_space. - :arg warn_if_cg: If *True*, warn if the space is "CG" - """ - # validate function space - if not isinstance(function_space, WithGeometry): - raise TypeError("function_space must be of type WithGeometry, not " - "%s" % type(function_space)) - - family = function_space.ufl_element().family() - acceptable_families = ['Discontinuous Lagrange', 'Lagrange'] - if family not in acceptable_families: - raise ValueError("function_space.ufl_element().family() must be " - "one of %s, not '%s'" % - (acceptable_families, family)) - if family == 'Lagrange' and warn_if_cg: - warn("Functions in continuous function space will be projected " - "to/from a 'Discontinuous Lagrange' space. Make sure " - "any operators evaluated are continuous. " - "Pass warn_if_cg=False to suppress this warning.") - - # validate potential_source_and_targets - if not isinstance(potential_source_and_target, - PotentialSourceAndTarget): - raise TypeError("potential_source_and_targets must be of type " - "PotentialSourceAndTarget, not '%s'." - % type(potential_source_and_target)) - - # make sure meshes match - if potential_source_and_target.mesh is not function_space.mesh(): - raise ValueError("function_space.mesh() and " - "potential_source_and_target.mesh must be the same" - " obejct.") - - # build DG space if necessary - family = function_space.ufl_element().family() - if family == 'Discontinuous Lagrange': - dg_function_space = function_space - elif family == 'Lagrange': - mesh = function_space.mesh() - degree = function_space.ufl_element().degree() - shape = function_space.shape - if shape is None or len(shape) == 0: - dg_function_space = FunctionSpace(mesh, "DG", degree) - elif len(shape) == 1: - dg_function_space = VectorFunctionSpace(mesh, "DG", degree, - dim=shape) - else: - dg_function_space = TensorFunctionSpace(mesh, "DG", degree, - shape=shape) - else: - acceptable_families = ['Discontinuous Lagrange', 'Lagrange'] - raise ValueError("function_space.ufl_element().family() must be " - "one of %s, not '%s'" % - (acceptable_families, family)) - - # store function space and dg function space - self.function_space = function_space - self.dg_function_space = dg_function_space - self.is_dg = function_space == dg_function_space - # store source and targets - self.source_and_target = potential_source_and_target - - def from_firedrake(self, density): - """ - Convert the density into a form acceptable by an bound operation - in an external library - - :arg density: A :class:`~firedrake.function.Function` holding the - density. - :returns: The converted density - """ - raise NotImplementedError - - def to_firedrake(self, evaluated_potential, out=None): - """ - Convert the evaluated potential from an external library - into a firedrake function - - :arg evaluated_potential: the evaluated potential - :arg out: If not *None*, store the converted potential into this - :class:`firedrake.function.Function`. - :return: *out* if it is not *None*, otherwise a new - :class:`firedrake.function.Function` storing the evaluated - potential - """ - raise NotImplementedError - - -class MeshmodeConnection(PotentialEvaluationLibraryConnection): - """ - Build a connection into :mod:`meshmode` - """ - def __init__(self, function_space, potential_source_and_target, actx, - warn_if_cg=True, - meshmode_connection_kwargs=None): - """ - For other args see :class:`PotentialEvaluationLibraryConnection` - - :arg actx: a :class:`meshmode.array_context.PyOpenCLArrayContext` - used for :func:`meshmode.interop.build_connection_from_firedrake` - and for conversion by the - :class:`meshmode.interop.FiredrakeConnection`. - :arg meshmode_connection_kwargs: A dict passed as kwargs - to :func:`meshmode.interop.build_connection_from_firedrake`. - Must not include 'restrict_to_boundary' - """ - PotentialEvaluationLibraryConnection.__init__( - self, - function_space, - potential_source_and_target, - warn_if_cg=warn_if_cg) - - # FIXME : allow subdomain regions - tdim = potential_source_and_target.mesh.topological_dimension() - if potential_source_and_target.get_source_dim() != tdim: - if potential_source_and_target.get_source_id() != "everywhere": - raise NotImplementedError("subdomain sources not implemented") - if potential_source_and_target.get_target_dim() != tdim: - if potential_source_and_target.get_target_id() != "everywhere": - raise NotImplementedError("subdomain targets not implemented") - - # validate actx - from meshmode.array_context import PyOpenCLArrayContext - if not isinstance(actx, PyOpenCLArrayContext): - raise TypeError("actx must be of type PyOpenCLArrayContext, not " - "%s." % type(actx)) - - # validate meshmode_connection_kwargs - if meshmode_connection_kwargs is None: - meshmode_connection_kwargs = {} - if not isinstance(meshmode_connection_kwargs, dict): - raise TypeError("meshmode_connection_kwargs must be *None* or of " - "type dict, not '%s'." % - type(meshmode_connection_kwargs)) - if 'restrict_to_boundary' in meshmode_connection_kwargs: - raise ValueError("'restrict_to_boundary' must not be set by " - "meshmode_connection_kwargs") - - # build a meshmode connection for the source - src_bdy = None - if potential_source_and_target.get_source_dim() != tdim: - # sanity check - assert potential_source_and_target.get_source_dim() == tdim - 1 - src_bdy = potential_source_and_target.get_source_id() - - from meshmode.interop.firedrake import build_connection_from_firedrake - src_conn = build_connection_from_firedrake(actx, - self.dg_function_space, - restrict_to_boundary=src_bdy, - **meshmode_connection_kwargs) - # If source is a boundary, build a connection to restrict it to the - # boundary - restrict_conn = None - if src_bdy is not None: - # convert "everywhere" to meshmode BTAG_ALL - meshmode_src_bdy = src_bdy - if meshmode_src_bdy == "everywhere": - from meshmode.mesh import BTAG_ALL - meshmode_src_bdy = BTAG_ALL - # get group factory - order = self.dg_function_space.degree() - from meshmode.discretization.poly_element import \ - PolynomialRecursiveNodesGroupFactory - default_grp_factory = \ - PolynomialRecursiveNodesGroupFactory(order, 'lgl') - grp_factory = meshmode_connection_kwargs.get('grp_factory', - default_grp_factory) - from meshmode.discretization.connection import make_face_restriction - restrict_conn = make_face_restriction(actx, - src_conn.discr, - grp_factory, - meshmode_src_bdy) - - # Build a meshmode connection for the target - tgt_bdy = None - if potential_source_and_target.get_target_dim() != tdim: - # sanity check - assert potential_source_and_target.get_target_dim() == tdim - 1 - tgt_bdy = potential_source_and_target.get_target_id() - - # Can we re-use the source connection? - src_dim = potential_source_and_target.get_source_dim() - tgt_dim = potential_source_and_target.get_target_dim() - if tgt_bdy == src_bdy and src_dim == tgt_dim: - tgt_conn = src_conn - else: - # If not, build a new one - tgt_conn = build_connection_from_firedrake( - actx, - self.dg_function_space, - restrict_to_boundary=tgt_bdy, - **meshmode_connection_kwargs) - - # store computing context - self.actx = actx - # store connections - self.source_to_meshmode_connection = src_conn - self.restrict_source_to_boundary = restrict_conn - self.target_to_meshmode_connection = tgt_conn - - def from_firedrake(self, density): - # make sure we are in the dg function space - if not self.is_dg: - density = project(density, self.dg_function_space) - # Convert to meshmode. - density = \ - self.source_to_meshmode_connection.from_firedrake(density, - actx=self.actx) - # restrict to boundary if necessary - if self.restrict_source_to_boundary is not None: - density = self.restrict_source_to_boundary(density) - return density - - def to_firedrake(self, evaluated_potential, out=None): - # if we are dg, it's simple - if self.is_dg: - return self.target_to_meshmode_connection.from_meshmode( - evaluated_potential, - out=out) - else: - # Otherwise, we have to project back to our function space - pot = \ - self.target_to_meshmode_connection.from_meshmode(evaluated_potential) - pot = project(evaluated_potential, self.function_space) - if out is not None: - out.dat.data[:] = pot.dat.data[:] - pot = out - return pot - - def get_source_discretization(self): - """ - Get the :class:`meshmode.discretization.Discretization` - of the source - """ - if self.restrict_source_to_boundary is None: - return self.source_to_meshmode_connection.discr - return self.restrict_source_to_boundary.to_discr - - def get_target_discretization(self): - """ - Get the :class:`meshmode.discretization.Discretization` - of the target - """ - return self.target_to_meshmode_connection.discr - - -class VolumentialConnection(PotentialEvaluationLibraryConnection): - """ - Build a connection into :mod:`volumential` - """ - def __init__(self, - potential_source_and_target, - function_space, - cl_ctx, - queue, - box_fmm_geometry_kwargs, - reorder_potentials, - meshmode_connection_kwargs, - warn_if_cg=True): - """ - For other args see :class:`PotentialEvaluationLibraryConnection`. - The *potential_source_and_target* must have a 3D - source, 3D target, and 3D mesh of co-dimension 0. - - :arg cl_ctx: A :class:`pyopencl.Context` - :arg queue: A :class:`pyopencl.CommandQueue` - :arg box_fmm_geometry_kwargs: kwargs passed to - :class:`volumential.geometry.BoxFMMGeometryFactory`. - Must include the key *'nlevels'* - :arg reorder_potentials: *True* iff the volumential FMM - is reordering potentials. *False* otherwise. - :arg meshmode_connection_kwargs: A dict passes as kwargs - to :func:`meshmode.interop.build_connection_from_firedrake`. - :mod:`meshmode` is used an intermediate level between - :mod:`firedrake` and :mod:`volumential`. - """ - # super - PotentialEvaluationLibraryConnection.__init__( - self, - function_space, - potential_source_and_target, - warn_if_cg=warn_if_cg) - - # validate that sources, target, and mesh are of right dimension - if potential_source_and_target.get_source_dimension() != 3: - raise ValueError("potential source must have dimension 3") - if potential_source_and_target.get_target_dimension() != 3: - raise ValueError("potential target must have dimension 3") - if potential_source_and_target.mesh().geometric_dimension() != 3: - raise ValueError("mesh must have geometric dimension 3") - if potential_source_and_target.mesh().topological_dimension() != 3: - raise ValueError("mesh must have topological dimension 3") - - # FIXME : allow subdomain regions - if potential_source_and_target.get_source_id() != "everywhere": - raise NotImplementedError("subdomain sources not implemented") - if potential_source_and_target.get_target_id() != "everywhere": - raise NotImplementedError("subdomain targets not implemented") - - # validate cl_ctx argument - from pyopencl import Context - if not isinstance(cl_ctx, Context): - raise TypeError("cl_ctx must be of type Context, not " - "%s." % type(cl_ctx)) - - # validate queue argument - from pyopencl import CommandQueue - if not isinstance(queue, CommandQueue): - raise TypeError("queue must be of type CommandQueue, not " - "%s." % type(queue)) - - # validate reorder_potentials - if not isinstance(reorder_potentials, bool): - raise TypeError("reorder_potentials must be of type bool, not " - "%s." % type(reorder_potentials)) - - # validate box_fmm_geometry_kwargs - if not isinstance(box_fmm_geometry_kwargs, dict): - raise TypeError("box_fmm_geometry_kwargs must be of " - "type dict, not '%s'." % - type(box_fmm_geometry_kwargs)) - if 'nlevels' not in box_fmm_geometry_kwargs: - raise ValueError("box_fmm_geometry_kwargs missing required keyword" - " 'nlevels'.") - - # validate meshmode_connection_kwargs - if meshmode_connection_kwargs is None: - meshmode_connection_kwargs = {} - if not isinstance(meshmode_connection_kwargs, dict): - raise TypeError("meshmode_connection_kwargs must be *None* or of " - "type dict, not '%s'." % - type(meshmode_connection_kwargs)) - - # Build intermediate connection into meshmode - from meshmode.interop.firedrake import build_connection_from_firedrake - from meshmode.array_context import PyOpenCLArrayContext - actx = PyOpenCLArrayContext(queue) - meshmode_connection = MeshmodeConnection( - function_space, - potential_source_and_target, - actx, - warn_if_cg=warn_if_cg, - meshmode_connection_kwargs=meshmode_connection_kwargs) - - # Build connection from meshmode into volumential - # (following https://gitlab.tiker.net/xywei/volumential/-/blob/fe2c3e7af355d5c527060e783237c124c95397b5/test/test_interpolation.py#L72 ) # noqa : E501 - from volumential.geometry import ( - BoundingBoxFactory, BoxFMMGeometryFactory) - from volumential.interpolation import ElementsToSourcesLookupBuilder - dim = function_space.mesh().geometric_dimension() - bboxfmm_fac = BoxFMMGeometryFactory(cl_ctx, dim=dim, - **box_fmm_geometry_kwargs) - boxgeo = bboxfmm_fac(queue) - elt_to_src_lookup_fac = ElementsToSourcesLookupBuilder( - cl_ctx, tree=boxgeo.tree, discr=meshmode_connection.discr) - elt_to_src_lookup, evt = elt_to_src_lookup_fac(queue) - - # Build connection from volumential into meshmode - from volumential.interpolation import LeavesToNodesLookupBuilder - leaves_to_node_lookup_fac = LeavesToNodesLookupBuilder( - cl_ctx, trav=boxgeo.trav, discr=meshmode_connection.discr) - leaves_to_node_lookup, evt = leaves_to_node_lookup_fac(queue) - - # Figure out whether or not conversion from volumential -> meshmode - # should user order or tree order - order = "user" if reorder_potentials else "tree" - - # Store maps for conversion to/from volumential - self.meshmode_connection = meshmode_connection - self.to_volumential_lookup = elt_to_src_lookup - self.from_volumential_lookup = leaves_to_node_lookup - self.from_volumential_order = order - # store necessary computing contextx - self.queue = queue - self.actx = actx - - def from_firedrake(self, density): - # pass operand into a meshmode DOFArray - from meshmode.dof_array import flatten - meshmode_density = self.meshmode_connection.from_firedrake(density) - meshmode_density = flatten(meshmode_density) - - # pass flattened operand into volumential interpolator - from volumential.interpolation import interpolate_from_meshmode - volumential_density = \ - interpolate_from_meshmode(self.queue, - meshmode_density, - self.to_volumential_lookup, - order="user") # user order is more intuitive - return volumential_density - - def to_firedrake(self, evaluated_potential, out=None): - # pass volumential back to meshmode DOFArray - from volumential.interpolation import interpolate_to_meshmode - from meshmode.dof_array import unflatten - meshmode_pot_vals = \ - interpolate_to_meshmode(self.queue, - evaluated_potential, - self.from_volumential_lookup, - order=self.from_volumential_order) - meshmode_pot_vals = unflatten(self.actx, - self.meshmode_connection.discr, - meshmode_pot_vals) - # get meshmode data back as firedrake function - return self.meshmode_connection.to_firedrake(meshmode_pot_vals, out=self) - - -""" -functions - SingleLayerPotential - DoubleLayerPotential - VolumentialOperation - VolumePotential -""" - - -class PotentialSourceAndTarget: - """ - Holds the source and target for a layer or volume potential - """ - def __init__(self, mesh, - source_region_dim=None, - source_region_id=None, - target_region_dim=None, - target_region_id=None): - """ - Source and target of a layer or volume potential. - The mesh must have co-dimension 0 or 1 - - By region_dim, we mean the topological dimension of the - region. - Regions must have co-dimensions 0 or 1. - They cannot have a higher dimension than the mesh. - *None* indicates the topological dimension of the - mesh. - - Region ids must be either a valid mesh subdomain id (an - integer or tuple) or *None*. *None* indicates either - the entire mesh, or the entire exterior boundary - (as determined by the value of region). - The string "everywhere" is also equivalent to a - value of *None* - - NOTE: For implementation reasons, if the target is - a boundary, instead of just evaluating the - target at DOFs on that boundary, it is instead - evaluated at any DOF of a cell which has at least - one vertex on the boundary. - """ - # validate mesh - if not isinstance(mesh, MeshGeometry): - raise TypeError("mesh must be of type MeshGeometry, not '%s'." % - type(mesh)) - - # mesh must have co-dimension 1 or 0. - valid_codims = [0, 1] - codim = mesh.geometric_dimension() - mesh.topological_dimension() - if codim not in valid_codims: - raise ValueError("mesh has invalid co-dimension of %s. " - "co-dimension must be one of %s." % - (codim, valid_codims)) - - # validate dims - for dim in [source_region_dim, target_region_dim]: - if dim is None: - continue - if not isinstance(dim, int): - raise TypeError("source and target dimensions must be *None*" - " or an integer, not of type '%s'." % type(dim)) - valid_dims = set([mesh.topological_dimension(), - mesh.geometric_dimension()-1]) - if dim < mesh.geometric_dimension() - 1: - raise ValueError("source and target dimensions must be " - "one of %s or *None*, not '%s'." % - (valid_dims, dim)) - # Get dims if *None* - if source_region_dim is None: - source_region_dim = mesh.topological_dimension() - if target_region_dim is None: - target_region_dim = mesh.topological_dimension() - - # Validate region ids type and validity - for id_, dim in zip([source_region_id, target_region_id], - [source_region_dim, target_region_dim]): - if id_ is None or id_ == "everywhere": - continue - # standardize id_ to a tuple - if isinstance(id_, int): - id_ = tuple(id_) - # check type - if not isinstance(id_, tuple): - raise TypeError("source and target region ids must be " - "*None*, an int, or a tuple, not of type") - # boundary case: - if dim == mesh.topological_dimension() - 1: - if not set(id_) <= set(mesh.exterior_facets.unique_markers): - raise ValueError(("boundary region ids %s are not a " - + "subset of mesh.exterior_facets." - + "unique_markers = %s.") % (id_, - mesh.exterior_facets.unique_markers)) - else: - # sanity check - assert dim == mesh.topological_dimension() - # this caches the cell subset on the mesh, which we'll - # probably be doing anyway if we're targeting the region - # - # It also throws a ValueError if the id_ is invalid - mesh.cell_subset(id_) - - # handle None case - if source_region_id is None: - source_region_id = "everywhere" - if target_region_id is None: - target_region_id = "everywhere" - - # store mesh, region dims, and region ids - self.mesh = mesh - self._source_region_dim = source_region_dim - self._target_region_dim = source_region_dim - self._source_region_id = source_region_id - self._target_region_id = source_region_id - - def get_source_dimension(self): - """ - Get the topological dimension of the source region - """ - return self._source_region_dim - - def get_target_dimension(self): - """ - Get the topological dimension of the target region - """ - return self._target_region_dim - - def get_source_id(self): - """ - Get the subdomain id of the source, or the string - "everywhere" to represent the entire exterior boundary/mesh - """ - return self._source_region_id - - def get_target_id(self): - """ - Get the subdomain id of the target, or the string - "everywhere" to represent the entire exterior boundary/mesh - """ - return self._target_region_id - - -def PytentialOperation(actx, - density, - unbound_op, - density_name, - potential_src_and_tgt, - **kwargs): - """ - :arg actx: A :class:`meshmode.dof_array.PyOpenCLArrayContext` - :arg density: the :mod:`firedrake` density function - :arg unbound_op: A :mod:`pytential` operation which has not been - bound (e.g. a :mod:`pymbolic` expression) - :arg density_name: A string, the name of the density function - in the unbound_op - :arg potential_src_and_tgt: A :class:`PotentialSourceAndTarget` - - :kwarg warn_if_cg: Pass *False* to suppress the warning if - a "CG" space is used - :kwarg meshmode_connection_kwargs: *None*, or - a dict with arguments to pass to - :func:`meshmode.interop.firedrake.build_connection_from_firedrake` - :kwarg qbx_kwargs: *None*, or a dict with arguments to pass to - :class:`pytential.qbx.QBXLayerPotentialSource`. - :kwarg op_kwargs: kwargs passed to the invocation of the - bound pytential operator (e.g. {'k'=0.5} for - a :class:`sumpy.kernel.HelmholtzKernel`). - Must not include *density_name* - """ - # make sure density name is a string - if not isinstance(density_name, str): - raise TypeError("density_name must be of type str, not '%s'." % - density_name) - - # get kwargs to build connection - warn_if_cg = kwargs.get('warn_if_cg', True) - meshmode_connection_kwargs = kwargs.get('meshmode_connection_kwargs', None) - - # Build meshmode connection - meshmode_connection = MeshmodeConnection( - density.function_space(), - potential_src_and_tgt, - actx, - warn_if_cg=warn_if_cg, - meshmode_connection_kwargs=meshmode_connection_kwargs) - - # Build QBX - src_discr = meshmode_connection.get_source_discretization() - qbx_kwargs = kwargs.get('qbx_kwargs', None) - from pytential.qbx import QBXLayerPotentialSource - qbx = QBXLayerPotentialSource(src_discr, **qbx_kwargs) - - # Bind pytential operator - tgt_discr = meshmode_connection.get_target_discretization() - from pytential import bind - pyt_op = bind((qbx, tgt_discr), unbound_op) - - # Get operator kwargs - op_kwargs = kwargs.get('op_kwargs', {}) - if density_name in op_kwargs: - raise ValueError(f"density_name '{density_name}' should not be included" - " in op_kwargs.") - - # build bound operator that just takes density as argument - def bound_op_with_kwargs(density_arg): - op_kwargs[density_name] = density_arg - return pyt_op(**op_kwargs) - - # Now build and return Potential object - return Potential(density, - connection=meshmode_connection, - potential_operator=bound_op_with_kwargs) - - -def SingleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): - pass - - -def DoubleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): - pass - - -def VolumentialOperation(density, unbound_op, potential_src_and_tgt, - function_space=None): - pass - - -class VolumePotential(AbstractExternalOperator): - r""" - Evaluates to - - .. math:: - - f(x) = \int_\Om K(x-y) op(y) \,dx - - where op is a function-type living in *function_space*, \Om is the - 3-D mesh on which *function_space* lives, and K - a :mod:`sumpy` kernel. - - :arg operator_data: A map with required keys - - * 'kernel': a :class:`sumpy.kernel.Kernel` (*K* in the math above) - * 'kernel_type': A string describing the kernel, e.g. ``"Laplace"`` - * 'cl_ctx': a :mod:`pyopencl` computing context - * 'queue': a :mod:`pyopencl` command queue compatible with the - computing context - * 'nlevels': The desired number of levels to be used by - :mod:`volumential` - * 'm_order': multipole order used by the expansion - * 'dataset_filename': the filename to pass to - :class:`volumential.table_manager.NearFieldInteractionTableManager` - - And optional keys - - * 'grp_factory': (optional) An interpolatory group factory - inheriting from :class:`meshmode.discretization.ElementGroupFactory` - to be used in the intermediate :mod:`meshmode` representation - * 'q_order': (optional) The desired :mod:`volumential` quadrature - order, defaults to *function_space*'s degree - * 'force_direct_evaluation': (optional) As in - :func:`volumential.volume_fmm.drive_volume_fmm`. - Defaults to *False* - * 'fmm_kwargs': (optional) A dictionary of kwargs - to pass to :func:`volumential.volume_fmm.drive_fmm`. - * 'root_extent': (optional) the root extent to pass to - :class:`volumential.table_manager.NearFieldInteractionTableManager` - (defaults to 1) - * 'table_compute_method': (optional) used by - :class:`volumential.table_manager.NearFieldInteractionTableManager` - (defaults to "DrosteSum") - * 'table_kwargs': (optional) passed to - :class:`volumential.table_manager.NearFieldInteractionTableManager` - 's *get_table* method - - """ - - _external_operator_type = 'GLOBAL' - - def __init__(self, *operands, **kwargs): - self.operands = operands - self.kwargs = kwargs - AbstractExternalOperator.__init__(self, operands[0], **kwargs) - - - # Build a broken version of input function space - fs = kwargs["function_space"] - fsdg = FunctionSpace(fs.mesh(), "DG", fs.ufl_element().degree()) - - self.evaluator_kwargs = kwargs.copy() - self.evaluator_kwargs["function_space"] = fsdg - - # Validate input - # print(len(self.operands), self.derivatives) - # assert self.derivatives == (0,), \ - # "Derivatives of volume potential not currently supported: " + str(self.derivatives) - # from firedrake import Function - # - # if not isinstance(operand, Function): - # raise TypeError(":arg:`operand` must be of type firedrake.Function" - # ", not %s." % type(operand)) - # if operand.function_space().shape != tuple(): - # raise ValueError(":arg:`operand` must be a function with shape ()," - # " not %s." % operand.function_space().shape) - operator_data = kwargs["operator_data"] - assert isinstance(operator_data, dict) - required_keys = ('kernel', 'kernel_type', 'cl_ctx', 'queue', 'nlevels', - 'm_order', 'dataset_filename') - optional_keys = ('grp_factory', 'q_order', 'force_direct_evaluation', - 'fmm_kwargs', 'root_extent', - 'table_compute_method', 'table_kwargs', - ) - permissible_keys = required_keys + optional_keys - if not all(key in operator_data for key in required_keys): - raise ValueError("operator_data is missing one of the required " - "keys: %s" % required_keys) - if not all(key in permissible_keys for key in operator_data): - raise ValueError("operator_data contains an unexpected key. All " - "keys must be one of %s." % (permissible_keys,)) - - @cached_property - def _evaluator(self): - return VolumePotentialEvaluator(self, *self.operands, **self.evaluator_kwargs) - - def _evaluate(self, continuity_tolerance=None): - 1/0 - - def _compute_derivatives(self, continuity_tolerance=None): - 1/0 - - def _evaluate_action(self, *args, continuity_tolerance=None): - return self._evaluator._evaluate_action(continuity_tolerance=continuity_tolerance) - - def _evaluate_adjoint_action(self, *args, continuity_tolerance=None): - return self._evaluator._evaluate_action(continuity_tolerance=continuity_tolerance) - - def evaluate_adj_component_control(self, x, idx): - derivatives = tuple(dj + int(idx == j) for j, dj in enumerate(self.derivatives)) - dN = self._ufl_expr_reconstruct_(*self.ufl_operands, derivatives=derivatives) - return dN._evaluate_adjoint_action((x.function,)).vector() - - def evaluate_adj_component(self, x, idx): - print(x, type(x)) - 1/0 - - -class VolumePotentialEvaluator: - def __init__(self, vp, *operands, **kwargs): - self.vp = vp - operand = operands[0] - - cg_function_space = vp.function_space() - function_space = kwargs["function_space"] - self.dg_function_space = function_space - - # Create an interpolator from the original operand so that - # we can re-interpolate it into the space at each evaluation - # in case it has changed (e.g. cost functional in PDE constrained optimization. - - operator_data = kwargs["operator_data"] - kernel = operator_data['kernel'] - kernel_type = operator_data['kernel_type'] - cl_ctx = operator_data['cl_ctx'] - queue = operator_data['queue'] - nlevels = operator_data['nlevels'] - m_order = operator_data['m_order'] - dataset_filename = operator_data['dataset_filename'] - degree = function_space.ufl_element().degree() - grp_factory = operator_data.get('grp_factory', None) - q_order = operator_data.get('q_order', degree) - force_direct_evaluation = operator_data.get('force_direct_evaluation', - False) - fmm_kwargs = operator_data.get('fmm_kwargs', {}) - root_extent = operator_data.get("root_extent", 1) - table_compute_method = operator_data.get('table_compute_method', - 'DrosteSum') - table_kwargs = operator_data.get('table_kwargs', {}) - - from sumpy.kernel import Kernel - if not isinstance(kernel, Kernel): - raise TypeError("operator_data['kernel'] must be of type " - "sumpy.kernel.Kernel, not %s." % type(kernel)) - if not isinstance(nlevels, int): - raise TypeError("operator_data['nlevels'] must be of type int, " - "not %s." % type(nlevels)) - if not isinstance(m_order, int): - raise TypeError("operator_data['m_order'] must be of type int, " - "not %s." % type(m_order)) - if not isinstance(q_order, int): - raise TypeError("operator_data['q_order'] must be of type int, " - "not %s." % type(q_order)) - - # Build connection into meshmode - from meshmode.interop.firedrake import build_connection_from_firedrake - from meshmode.array_context import PyOpenCLArrayContext - actx = PyOpenCLArrayContext(queue) - meshmode_connection = build_connection_from_firedrake( - actx, function_space, grp_factory=grp_factory) - - # Build connection from meshmode into volumential - # (following https://gitlab.tiker.net/xywei/volumential/-/blob/fe2c3e7af355d5c527060e783237c124c95397b5/test/test_interpolation.py#L72 ) # noqa : E501 - from volumential.geometry import ( - BoundingBoxFactory, BoxFMMGeometryFactory) - from volumential.interpolation import ElementsToSourcesLookupBuilder - dim = function_space.mesh().geometric_dimension() - bbox_fac = BoundingBoxFactory(dim=dim) - bboxfmm_fac = BoxFMMGeometryFactory( - cl_ctx, dim=dim, order=q_order, - nlevels=nlevels, bbox_getter=bbox_fac, - expand_to_hold_mesh=meshmode_connection.discr.mesh, - mesh_padding_factor=0.0) - boxgeo = bboxfmm_fac(queue) - elt_to_src_lookup_fac = ElementsToSourcesLookupBuilder( - cl_ctx, tree=boxgeo.tree, discr=meshmode_connection.discr) - elt_to_src_lookup, evt = elt_to_src_lookup_fac(queue) - - # Build connection from volumential into meshmode - from volumential.interpolation import LeavesToNodesLookupBuilder - leaves_to_node_lookup_fac = LeavesToNodesLookupBuilder( - cl_ctx, trav=boxgeo.trav, discr=meshmode_connection.discr) - leaves_to_node_lookup, evt = leaves_to_node_lookup_fac(queue) - - # Create near-field table in volumential - from volumential.table_manager import NearFieldInteractionTableManager - table_manager = NearFieldInteractionTableManager( - dataset_filename, root_extent=root_extent, queue=queue) - - nftable, _ = table_manager.get_table( - dim, - kernel_type, - q_order, - compute_method=table_compute_method, - queue=queue, - **table_kwargs) - - # Create kernel wrangler in volumential - from sumpy.expansion import DefaultExpansionFactory - out_kernels = [kernel] - expn_factory = DefaultExpansionFactory() - local_expn_class = expn_factory.get_local_expansion_class(kernel) - mpole_expn_class = expn_factory.get_multipole_expansion_class(kernel) - - from functools import partial - from volumential.expansion_wrangler_fpnd import ( - FPNDExpansionWrangler, FPNDExpansionWranglerCodeContainer) - wrangler_cc = FPNDExpansionWranglerCodeContainer( - cl_ctx, - partial(mpole_expn_class, kernel), - partial(local_expn_class, kernel), - out_kernels) - wrangler = FPNDExpansionWrangler( - code_container=wrangler_cc, - queue=queue, - tree=boxgeo.tree, - near_field_table=nftable, - dtype=ScalarType, - fmm_level_to_order=lambda kernel, kernel_args, tree, lev: m_order, - quad_order=q_order) - - # Store attributes that we may need later - # self.ufl_operands = operands - self.actx = actx - self.queue = queue - self.sumpy_kernel = kernel - self.meshmode_connection = meshmode_connection - self.volumential_boxgeo = boxgeo - self.to_volumential_lookup = elt_to_src_lookup - self.from_volumential_lookup = leaves_to_node_lookup - self.force_direct_evaluation = force_direct_evaluation - self.fmm_kwargs = fmm_kwargs - self.expansion_wrangler = wrangler - # so we don't have to keep making a fd function during conversion - # from meshmode. AbstractExternalOperator s aren't functions, - # so can't use *self* - self.fd_pot = Function(self.dg_function_space) - self.fd_pot_cg = Function(vp.function_space()) - - def _evaluate(self, continuity_tolerance=None): - """ - :arg continuity_tolerance: If *None* then ignored. Otherwise, a floating - point value. If the function space associated to this operator - is a continuous one and *continuity_tolerance* is a float, it is - verified that, - during the conversion of the evaluated potential - from meshmode's discontinuous function space representation into - firedrake's continuous one, the potential's value at - a firedrake node is within *continuity_tolerance* of its - value at any duplicated firedrake node - """ - operand, = self.vp.ufl_operands - operand_discrete = interpolate(operand, self.dg_function_space) - - # pass operand into a meshmode DOFArray - from meshmode.dof_array import flatten - meshmode_src_vals = \ - self.meshmode_connection.from_firedrake(operand_discrete, actx=self.actx) - meshmode_src_vals = flatten(meshmode_src_vals) - - # if 1: - # alpha = 40 - # dim = 2 - - # if 1: - # mm_discr = self.meshmode_connection.discr - # from meshmode.dof_array import thaw, flat_norm - # x = flatten(thaw(self.actx, mm_discr.nodes()[0])) - # y = flatten(thaw(self.actx, mm_discr.nodes()[1])) - # import pyopencl.clmath as clmath - # import pyopencl.array as cla - # norm2 = (x-0.5)**2 + (y-0.5)**2 - # ref_src = -(4 * alpha ** 2 * norm2 - 2 * dim * alpha) * clmath.exp( - # -alpha * norm2) - # denom = np.max(np.abs(ref_src.get())) - # err = np.abs((ref_src - meshmode_src_vals).get()) - # print("Fdrake -> meshmode error:", np.max(err)/denom) - - # pass flattened operand into volumential interpolator - from volumential.interpolation import interpolate_from_meshmode - volumential_src_vals = \ - interpolate_from_meshmode(self.queue, - meshmode_src_vals, - self.to_volumential_lookup, - order="user") # user order is more intuitive - - # if 1: - # # check source density against known value - # # x = self.volumential_boxgeo.tree.sources[0].with_queue(self.queue) - # # y = self.volumential_boxgeo.tree.sources[1].with_queue(self.queue) - # x, y = self.volumential_boxgeo.nodes - # import pyopencl.clmath as clmath - # import pyopencl.array as cla - # norm2 = (x-0.5)**2 + (y-0.5)**2 - # ref_src = -(4 * alpha ** 2 * norm2 - 2 * dim * alpha) * clmath.exp( - # -alpha * norm2) - # denom = np.max(np.abs(ref_src.get())) - # err = np.abs((ref_src-volumential_src_vals).get()) - # print("Fdrake -> volumential error:", np.max(err)/denom) - - # evaluate volume fmm - from volumential.volume_fmm import drive_volume_fmm - pot, = drive_volume_fmm( - self.volumential_boxgeo.traversal, - self.expansion_wrangler, - volumential_src_vals * self.volumential_boxgeo.weights, - volumential_src_vals, - direct_evaluation=self.force_direct_evaluation, - reorder_sources=True, - **self.fmm_kwargs) - - # if 1: - # # check potential against known value - # ref_pot = clmath.exp(-alpha * norm2) - # denom = np.max(np.abs(ref_pot.get())) - # err = np.abs((ref_pot-pot).get()) - # print("volumential potential error:", np.max(err)/denom) - - # # evaluate volume fmm with ref_src - # pot_w_ref_src, = drive_volume_fmm( - # self.volumential_boxgeo.traversal, - # self.expansion_wrangler, - # ref_src * self.volumential_boxgeo.weights, - # ref_src, - # direct_evaluation=self.force_direct_evaluation, - # reorder_sources=True, - # **self.fmm_kwargs) - - # err = np.abs((ref_pot-pot_w_ref_src).get()) - # print("volumential potential w/t ref src error:", np.max(err)/denom) - - # pass volumential back to meshmode DOFArray - from volumential.interpolation import interpolate_to_meshmode - from meshmode.dof_array import unflatten - user_order = self.fmm_kwargs.get("reorder_potentials", True) - if user_order: - order = "user" - else: - order = "tree" - meshmode_pot_vals = interpolate_to_meshmode(self.queue, - pot, - self.from_volumential_lookup, - order=order) - meshmode_pot_vals = unflatten(self.actx, - self.meshmode_connection.discr, - meshmode_pot_vals) - # get meshmode data back as firedrake function - self.meshmode_connection.from_meshmode( - meshmode_pot_vals, - out=self.fd_pot) - - self.fd_pot_cg = project(self.fd_pot, self.vp.function_space()) - - self.vp.dat.data[:] = self.fd_pot_cg.dat.data[:] - - return self.fd_pot_cg - - def _compute_derivatives(self, continuity_tolerance=None): - # TODO : Support derivatives - return self._evaluate(continuity_tolerance=continuity_tolerance) - - def _evaluate_action(self, *args, continuity_tolerance=None): - # From tests/pointwiseoperator/test_point_expr.py - # https://github.com/firedrakeproject/firedrake/blob/c0d9b592f587fa8c7437f690da7a6595f6804c1b/tests/pointwiseoperator/test_point_expr.py # noqa - if len(args) == 0: - # Evaluate the operator - return self._evaluate(continuity_tolerance=continuity_tolerance) - - # Evaluate the Jacobian/Hessian action - operands = self.ufl_operands - operator = self._compute_derivatives(continuity_tolerance=continuity_tolerance) - expr = as_ufl(operator(*operands)) - if expr.ufl_shape == () and expr != 0: - var = VariableRuleset(self.ufl_operands[0]) - expr = expr*var._Id - elif expr == 0: - return self.assign(expr) - - for arg in args: - mi = indices(len(expr.ufl_shape)) - aa = mi - bb = mi[-len(arg.ufl_shape):] - expr = arg[bb] * expr[aa] - mi_tensor = tuple(e for e in mi if not (e in aa and e in bb)) - if len(expr.ufl_free_indices): - expr = as_tensor(expr, mi_tensor) - return self.interpolate(expr) diff --git a/firedrake/potential_evaluation/__init__.py b/firedrake/potential_evaluation/__init__.py new file mode 100644 index 0000000000..80d7465050 --- /dev/null +++ b/firedrake/potential_evaluation/__init__.py @@ -0,0 +1,416 @@ +import numpy as np +import firedrake.function + +from firedrake.functionspaceimpl import WithGeometry +from firedrake.mesh import MeshGeometry +from firedrake.pointwise_operators import AbstractExternalOperator +from firedrake.utils import cached_property + +from firedrake import Function, FunctionSpace, interpolate, Interpolator, \ + SpatialCoordinate, VectorFunctionSpace, TensorFunctionSpace, project + +from ufl.algorithms.apply_derivatives import VariableRuleset +from ufl.algorithms import extract_coefficients +from ufl.constantvalue import as_ufl +from ufl.core.multiindex import indices +from ufl.tensors import as_tensor + +from pyop2.datatypes import ScalarType +from warnings import warn + + +from firedrake.potential_evaluation.potentials import \ + SingleLayerPotential, DoubleLayerPotential, VolumePotential +__all__ = ("SingleLayerPotential", + "DoubleLayerPotential", + "VolumePotential", + "PotentialSourceAndTarget") + + +class Potential(AbstractExternalOperator): + r""" + This is a function which represents a potential + computed on a firedrake function. + + For a firedrake function :math:`u:\Omega\to\mathbb C^n` + with :math:`\Omega \subseteq \mathbb R^m`, + the potential :math:`P(u):\Omega \to\mathbb C^n` is a + convolution of :math:`u` against some kernel function. + More concretely, given a source region :math:`\Gamma \subseteq \Omega` + a target region :math:`\Sigma \subseteq \Omega`, and + a kernel function :math:`K`, we define + + .. math:: + + P(u)(x) = \begin{cases} + \int_\Gamma K(x, y) u(y) \,dy & x \in \Sigma + 0 & x \notin \Sigma + \end{cases} + """ + _external_operator_type = 'GLOBAL' + + def __init__(self, density, **kwargs): + """ + :arg density: A :mod:`firedrake` :class:`firedrake.function.Function` + or UFL expression which represents the density :math:`u` + :kwarg connection: A :class:`PotentialEvaluationLibraryConnection` + :kwarg potential_operator: The external potential evaluation library + bound operator. + """ + # super + AbstractExternalOperator.__init__(self, density, **kwargs) + + # FIXME + for order in self.derivatives: + assert order == 1, "Assumes self.derivatives = (1,..,1)" + + # Get connection & bound op and validate + connection = kwargs.get("connection", None) + if connection is None: + raise ValueError("Missing kwarg 'connection'") + if not isinstance(connection, PotentialEvaluationLibraryConnection): + raise TypeError("connection must be of type " + "PotentialEvaluationLibraryConnection, not %s." + % type(connection)) + self.connection = connection + + self.potential_operator = kwargs.get("potential_operator", None) + if self.potential_operator is None: + raise ValueError("Missing kwarg 'potential_operator'") + + # Get function space and validate aginst bound op + function_space = self.ufl_function_space() + assert isinstance(function_space, WithGeometry) # sanity check + if function_space is not self.potential_operator.function_space(): + raise ValueError("function_space must be same object as " + "potential_operator.function_space().") + + # Make sure density is a member of our function space, if it is + if isinstance(density, Function): + if density.function_space() is not function_space: + raise ValueError("density.function_space() must be the " + "same as function_space") + # Make sure the shapes match, at least + elif density.shape != function_space.shape: + raise ValueError("Shape mismatch between function_space and " + "density. %s != %s." % + (density.shape, function_space.shape)) + + @cached_property + def _evaluator(self): + return PotentialEvaluator(self, + self.density, + self.connection, + self.potential_operator) + + def _evaluate(self): + raise NotImplementedError + + def _compute_derivatives(self, continuity_tolerance=None): + raise NotImplementedError + + def _evaluate_action(self, *args): + return self._evaluator._evaluate_action() + + def _evaluate_adjoint_action(self, *args): + return self._evaluator._evaluate_action() + + def evaluate_adj_component_control(self, x, idx): + derivatives = tuple(dj + int(idx == j) for j, dj in enumerate(self.derivatives)) + dN = self._ufl_expr_reconstruct_(*self.ufl_operands, derivatives=derivatives) + return dN._evaluate_adjoint_action((x.function,)).vector() + + def evaluate_adj_component(self, x, idx): + print(x, type(x)) + raise NotImplementedError + + +class PotentialEvaluator: + """ + Evaluates a potential + """ + def __init__(self, density, connection, potential_operator): + """ + :arg density: The UFL/firedrake density function + :arg connection: A :class:`PotentialEvaluationLibraryConnection` + :arg potential_operator: The external potential evaluation library + bound operator. + """ + self.density = density + self.connection = connection + self.potential_operator = potential_operator + + def _eval_potential_operator(self, action_coefficients, out=None): + """ + Evaluate the potential on the action coefficients and return. + If *out* is not *None*, stores the result in *out* + + :arg action_coefficients: A :class:`~firedrake.function.Function` + to evaluate the potential at + :arg out: If not *None*, then a :class:`~firedrake.function.Function` + in which to store the evaluated potential + :return: *out* if it is not *None*, otherwise a new + :class:`firedrake.function.Function` storing the evaluated + potential + """ + density = self.connection.from_firedrake(action_coefficients) + potential = self.potential_operator(density) + return self.connection.to_firedrake(potential, out=out) + + def _evaluate(self): + """ + Evaluate P(self.density) into self + """ + return self._eval_potential_operator(self.density, out=self) + + def _compute_derivatives(self): + """ + Return a function + Derivative(P, self.derivatives, self.density) + """ + # FIXME : Assumes derivatives are Jacobian + return self._eval_potential_operator + + def _evaluate_action(self, action_coefficients): + """ + Evaluate derivatives of layer potential at action coefficients. + i.e. Derivative(P, self.derivatives, self.density) evaluated at + the action coefficients and store into self + """ + operator = self._compute_derivatives() + return operator(action_coefficients, out=self) + + +class PotentialEvaluationLibraryConnection: + """ + A connection to an external library for potential evaluation + """ + def __init__(self, function_space, potential_source_and_target, + warn_if_cg=True): + """ + Initialize self to work on the function space + + :arg function_space: The :mod:`firedrake` function space + (of type :class:`~firedrake.functionspaceimpl.WithGeometry`) on + which to convert to/from. Must be a 'Lagrange' or + 'Discontinuous Lagrange' space. + :arg potential_source_and_target: A :class:`PotentialSourceAndTarget`. + mesh must match that of function_space. + :arg warn_if_cg: If *True*, warn if the space is "CG" + """ + # validate function space + if not isinstance(function_space, WithGeometry): + raise TypeError("function_space must be of type WithGeometry, not " + "%s" % type(function_space)) + + family = function_space.ufl_element().family() + acceptable_families = ['Discontinuous Lagrange', 'Lagrange'] + if family not in acceptable_families: + raise ValueError("function_space.ufl_element().family() must be " + "one of %s, not '%s'" % + (acceptable_families, family)) + if family == 'Lagrange' and warn_if_cg: + warn("Functions in continuous function space will be projected " + "to/from a 'Discontinuous Lagrange' space. Make sure " + "any operators evaluated are continuous. " + "Pass warn_if_cg=False to suppress this warning.") + + # validate potential_source_and_targets + if not isinstance(potential_source_and_target, + PotentialSourceAndTarget): + raise TypeError("potential_source_and_targets must be of type " + "PotentialSourceAndTarget, not '%s'." + % type(potential_source_and_target)) + + # make sure meshes match + if potential_source_and_target.mesh is not function_space.mesh(): + raise ValueError("function_space.mesh() and " + "potential_source_and_target.mesh must be the same" + " obejct.") + + # build DG space if necessary + family = function_space.ufl_element().family() + if family == 'Discontinuous Lagrange': + dg_function_space = function_space + elif family == 'Lagrange': + mesh = function_space.mesh() + degree = function_space.ufl_element().degree() + shape = function_space.shape + if shape is None or len(shape) == 0: + dg_function_space = FunctionSpace(mesh, "DG", degree) + elif len(shape) == 1: + dg_function_space = VectorFunctionSpace(mesh, "DG", degree, + dim=shape) + else: + dg_function_space = TensorFunctionSpace(mesh, "DG", degree, + shape=shape) + else: + acceptable_families = ['Discontinuous Lagrange', 'Lagrange'] + raise ValueError("function_space.ufl_element().family() must be " + "one of %s, not '%s'" % + (acceptable_families, family)) + + # store function space and dg function space + self.function_space = function_space + self.dg_function_space = dg_function_space + self.is_dg = function_space == dg_function_space + # store source and targets + self.source_and_target = potential_source_and_target + + def from_firedrake(self, density): + """ + Convert the density into a form acceptable by an bound operation + in an external library + + :arg density: A :class:`~firedrake.function.Function` holding the + density. + :returns: The converted density + """ + raise NotImplementedError + + def to_firedrake(self, evaluated_potential, out=None): + """ + Convert the evaluated potential from an external library + into a firedrake function + + :arg evaluated_potential: the evaluated potential + :arg out: If not *None*, store the converted potential into this + :class:`firedrake.function.Function`. + :return: *out* if it is not *None*, otherwise a new + :class:`firedrake.function.Function` storing the evaluated + potential + """ + raise NotImplementedError + + +class PotentialSourceAndTarget: + """ + Holds the source and target for a layer or volume potential + """ + def __init__(self, mesh, + source_region_dim=None, + source_region_id=None, + target_region_dim=None, + target_region_id=None): + """ + Source and target of a layer or volume potential. + The mesh must have co-dimension 0 or 1 + + By region_dim, we mean the topological dimension of the + region. + Regions must have co-dimensions 0 or 1. + They cannot have a higher dimension than the mesh. + *None* indicates the topological dimension of the + mesh. + + Region ids must be either a valid mesh subdomain id (an + integer or tuple) or *None*. *None* indicates either + the entire mesh, or the entire exterior boundary + (as determined by the value of region). + The string "everywhere" is also equivalent to a + value of *None* + + NOTE: For implementation reasons, if the target is + a boundary, instead of just evaluating the + target at DOFs on that boundary, it is instead + evaluated at any DOF of a cell which has at least + one vertex on the boundary. + """ + # validate mesh + if not isinstance(mesh, MeshGeometry): + raise TypeError("mesh must be of type MeshGeometry, not '%s'." % + type(mesh)) + + # mesh must have co-dimension 1 or 0. + valid_codims = [0, 1] + codim = mesh.geometric_dimension() - mesh.topological_dimension() + if codim not in valid_codims: + raise ValueError("mesh has invalid co-dimension of %s. " + "co-dimension must be one of %s." % + (codim, valid_codims)) + + # validate dims + for dim in [source_region_dim, target_region_dim]: + if dim is None: + continue + if not isinstance(dim, int): + raise TypeError("source and target dimensions must be *None*" + " or an integer, not of type '%s'." % type(dim)) + valid_dims = set([mesh.topological_dimension(), + mesh.geometric_dimension()-1]) + if dim < mesh.geometric_dimension() - 1: + raise ValueError("source and target dimensions must be " + "one of %s or *None*, not '%s'." % + (valid_dims, dim)) + # Get dims if *None* + if source_region_dim is None: + source_region_dim = mesh.topological_dimension() + if target_region_dim is None: + target_region_dim = mesh.topological_dimension() + + # Validate region ids type and validity + for id_, dim in zip([source_region_id, target_region_id], + [source_region_dim, target_region_dim]): + if id_ is None or id_ == "everywhere": + continue + # standardize id_ to a tuple + if isinstance(id_, int): + id_ = tuple(id_) + # check type + if not isinstance(id_, tuple): + raise TypeError("source and target region ids must be " + "*None*, an int, or a tuple, not of type") + # boundary case: + if dim == mesh.topological_dimension() - 1: + if not set(id_) <= set(mesh.exterior_facets.unique_markers): + raise ValueError(("boundary region ids %s are not a " + + "subset of mesh.exterior_facets." + + "unique_markers = %s.") % (id_, + mesh.exterior_facets.unique_markers)) + else: + # sanity check + assert dim == mesh.topological_dimension() + # this caches the cell subset on the mesh, which we'll + # probably be doing anyway if we're targeting the region + # + # It also throws a ValueError if the id_ is invalid + mesh.cell_subset(id_) + + # handle None case + if source_region_id is None: + source_region_id = "everywhere" + if target_region_id is None: + target_region_id = "everywhere" + + # store mesh, region dims, and region ids + self.mesh = mesh + self._source_region_dim = source_region_dim + self._target_region_dim = source_region_dim + self._source_region_id = source_region_id + self._target_region_id = source_region_id + + def get_source_dimension(self): + """ + Get the topological dimension of the source region + """ + return self._source_region_dim + + def get_target_dimension(self): + """ + Get the topological dimension of the target region + """ + return self._target_region_dim + + def get_source_id(self): + """ + Get the subdomain id of the source, or the string + "everywhere" to represent the entire exterior boundary/mesh + """ + return self._source_region_id + + def get_target_id(self): + """ + Get the subdomain id of the target, or the string + "everywhere" to represent the entire exterior boundary/mesh + """ + return self._target_region_id diff --git a/firedrake/potential_evaluation/potentials.py b/firedrake/potential_evaluation/potentials.py new file mode 100644 index 0000000000..275dffb710 --- /dev/null +++ b/firedrake/potential_evaluation/potentials.py @@ -0,0 +1,212 @@ +from firedrake.functionspaceimpl import WithGeometry +from firedrake.potential_evaluation.pytential import PytentialOperation + + +def SingleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): + """ + Evaluate the single layer potential of the density + function convoluted against the kernel on the source + at each target point. + + :arg density: The :mod:`firedrake` density function + :arg kernel: the :class:`sumpy.kernel.Kernel` + :arg potential_src_and_tgt: A + :class:`firedrake.potential_evaluation.PotentialSourceAndTarget` + + :kwarg cl_ctx: (Optional) a :class:`pyopencl.Context` used + to create a command queue. + At most one of *cl_ctx*, *queue*, *actx* should be included. + If none is included, a *cl_ctx* will be created. + :kwarg queue: (Optional) a :class:`pyopencl.CommandQueue` used + to create an array context. + At most one of *cl_ctx*, *queue*, *actx* should be included. + If none is included, a *cl_ctx* will be created. + :kwarg actx: (Optional) a :class:`meshmode.array_context.PyOpenCLArrayContext` + used for :mod:`pytential` and :mod:`meshmode` computations. + At most one of *cl_ctx*, *queue*, *actx* should be included. + If none is included, a *cl_ctx* will be created. + :kwarg qbx_forced_limit: (Optional) As described in + https://documen.tician.de/pytential/_modules/pytential/symbolic/primitives.html#IntG + DEFAULT: *None* + :kwarg op_kwargs: (Optional) kwargs passed to the invocation of the + bound pytential operator (e.g. {'k': 0.5} for + a :class:`sumpy.kernel.HelmholtzKernel`). + :kwarg qbx_kwargs: (Optional) *None*, or a dict with arguments to pass to + :class:`pytential.qbx.QBXLayerPotentialSource`. + DEFAULTS: + - fine order: 4 8 function space degree + - qbx_order: function space degree + - fmm_order: *None* + - fmm_level_to_order: A + :class:`sumpy.level_to_order.FMMLibExpansionOrderFinder` set + to provide tolerance of double machine epsilon + 2**-53 with an extra order of 1. + This allows for taking 1 derivative, but should be + increased for taking more derivatives. + - fmm_backend: "fmmlib" + :kwarg meshmode_connection_kwargs: (Optional) + *None*, or a dict with arguments to pass to + :func:`meshmode.interop.firedrake.build_connection_from_firedrake` + :kwarg warn_if_cg: (Optional) Pass *False* to suppress the warning if + a "CG" space is used. *True* by default + """ + # TODO: Make fmm_level_to_order depend on number of derivatives somehow? + from pytential.symbolic.primitives import S + _layer_potential(S, density, kernel, potential_src_and_tgt, **kwargs) + + +def DoubleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): + """ + Same as SingleLayerPotential, but evaluates the double layer potential + instead of the single + """ + from pytential.symbolic.primitives import D + _layer_potential(D, density, kernel, potential_src_and_tgt, **kwargs) + + +def _layer_potential(layer_potential_sym, + density, kernel, potential_src_and_tgt, **kwargs): + """ + Build a layer potential. For usage, see + :func:`SingleLayerPotential` or :func:`DoubleLayerPotential`. + """ + kwargs = _validate_layer_potential_args_and_set_defaults( + density, kernel, potential_src_and_tgt, **kwargs) + qbx_forced_limit = kwargs.get('qbx_forced_limit', None) + # build our unbound operator + from pytential import sym + unbound_op = layer_potential_sym(kernel, + sym.var("density"), + qbx_forced_limit=qbx_forced_limit) + # make sure we got an actx + assert 'actx' in kwargs + actx = kwargs['actx'] + # extract Pytential operation kwargs + pyt_op_kwargs = _extract_pytential_operation_kwargs(**kwargs) + # now return the pytential operation as an external operator + return PytentialOperation(actx, density, unbound_op, "density", + potential_src_and_tgt, **pyt_op_kwargs) + + +def _validate_layer_potential_args_and_set_defaults(density, + kernel, + places, + **kwargs): + """ + Validate the arguments for single/double layer potential. + + Returns dictionary updated with default arguments + """ + # validate density function space + if not isinstance(density.function_space(), WithGeometry): + raise TypeError("density.function_space() must be of type " + f"WithGeometry, not {type(density.function_space())}") + + # validate kernel + from sumpy.kernel import Kernel + if not isinstance(kernel, Kernel): + raise TypeError("kernel must be a sumpy.kernel.Kernel, not of " + f"type '{type(kernel)}'.") + # validate src and tgts + from firedrake.potential_evaluation import PotentialSourceAndTarget + if not isinstance(places, PotentialSourceAndTarget): + raise TypeError("potential_src_and_tgt must be a sumpy.kernel.Kernel, " + f"not of type '{type(places)}'.") + # Make sure src is of right dimension + mesh_gdim = places.mesh.geometric_dimension() + mesh_tdim = places.mesh.topological_dimension() + src_tdim = places.get_source_dimension() + # sanity check + assert mesh_gdim - mesh_tdim in [0, 1] + assert mesh_gdim - src_tdim in [0, 1] + assert mesh_tdim in (mesh_gdim, src_tdim) + # now do the real user-input check + if mesh_gdim - src_tdim != 1: + raise ValueError("source of a layer potential must have co-dimension 1," + f" not {mesh_gdim - src_tdim}.") + + # Make sure all kwargs are recognized + allowed_kwargs = ("cl_ctx", "queue", "actx", "op_kwargs", "qbx_kwargs", + "meshmode_connection_kwargs") + for key in kwargs: + if key not in allowed_kwargs: + raise ValueError(f"Unrecognized kwarg {key}") + + # Now handle pyopencl computing contexts and build + # a PyOpenCLArrayContext + from meshmode.array_context import PyOpenCLArrayContext + cl_ctx = None + queue = None + actx = None + if 'cl_ctx' in kwargs: + if 'actx' in kwargs or 'queue' in kwargs: + raise ValueError("At most one of 'actx', 'queue', 'cl_ctx' should " + "be supplied") + cl_ctx = kwargs['cl_ctx'] + from pyopencl import Context + if not isinstance(cl_ctx, Context): + raise TypeError("cl_ctx must be of type Context, not " + f"{type(cl_ctx)}") + queue = None + actx = None + elif 'queue' in kwargs: + if 'actx' in kwargs: + raise ValueError("At most one of 'actx', 'queue' should " + "be supplied") + queue = kwargs['queue'] + from pyopencl import CommandQueue + if not isinstance(queue, CommandQueue): + raise TypeError("queue must be of type CommandQueue, not " + f"{type(queue)}") + actx = None + elif 'actx' in kwargs: + actx = kwargs['actx'] + if not isinstance(actx, PyOpenCLArrayContext): + raise TypeError("actx must be of type PyOpenCLArrayContext, not " + f"{type(actx)}") + + # now make sure we actually get an actx in kwargs + if actx is None: + if queue is None: + if cl_ctx is None: + from pyopencl import create_some_context + cl_ctx = create_some_context() + from pyopencl import CommandQueue + queue = CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + kwargs['actx'] = actx + + # Set qbx_kwargs defaults + if 'qbx_kwargs' not in kwargs: + degree = density.function_space().ufl_element().degree() + from sumpy.expansion.level_to_order import FMMLibExpansionOrderFinder + qbx_kwargs = {'fine order': 4 * degree, + 'qbx_order': degree, + 'fmm_order': None, + 'fmm_level_to_order': + FMMLibExpansionOrderFinder(2**-53, extra_order=1), + 'fmm_backend': "fmmlib", + } + kwargs['qbx_kwargs'] = qbx_kwargs + + return kwargs + + +def _extract_pytential_operation_kwargs(**kwargs): + """ + Extract kwargs to be passed to :func:`PytentialOperation` + """ + pyt_op_kwargs = {} + pyt_op_possible_keywords = ("warn_if_cg", + "meshmode_connection_kwargs", + "qbx_kwargs", + "op_kwargs") + for key in pyt_op_possible_keywords: + if key in kwargs: + pyt_op_kwargs[key] = kwargs[key] + + return pyt_op_kwargs + + +def VolumePotential(density, kernel, potential_src_and_tgt, **kwargs): + raise NotImplementedError diff --git a/firedrake/potential_evaluation/pytential.py b/firedrake/potential_evaluation/pytential.py new file mode 100644 index 0000000000..ad6b14b291 --- /dev/null +++ b/firedrake/potential_evaluation/pytential.py @@ -0,0 +1,233 @@ +from firedrake.potential_evaluation import \ + Potential, PotentialEvaluationLibraryConnection +from firedrake import project + + +class MeshmodeConnection(PotentialEvaluationLibraryConnection): + """ + Build a connection into :mod:`meshmode` + """ + def __init__(self, function_space, potential_source_and_target, actx, + warn_if_cg=True, + meshmode_connection_kwargs=None): + """ + For other args see :class:`PotentialEvaluationLibraryConnection` + + :arg actx: a :class:`meshmode.array_context.PyOpenCLArrayContext` + used for :func:`meshmode.interop.build_connection_from_firedrake` + and for conversion by the + :class:`meshmode.interop.FiredrakeConnection`. + :arg meshmode_connection_kwargs: A dict passed as kwargs + to :func:`meshmode.interop.build_connection_from_firedrake`. + Must not include 'restrict_to_boundary' + """ + PotentialEvaluationLibraryConnection.__init__( + self, + function_space, + potential_source_and_target, + warn_if_cg=warn_if_cg) + + # FIXME : allow subdomain regions + tdim = potential_source_and_target.mesh.topological_dimension() + if potential_source_and_target.get_source_dim() != tdim: + if potential_source_and_target.get_source_id() != "everywhere": + raise NotImplementedError("subdomain sources not implemented") + if potential_source_and_target.get_target_dim() != tdim: + if potential_source_and_target.get_target_id() != "everywhere": + raise NotImplementedError("subdomain targets not implemented") + + # validate actx + from meshmode.array_context import PyOpenCLArrayContext + if not isinstance(actx, PyOpenCLArrayContext): + raise TypeError("actx must be of type PyOpenCLArrayContext, not " + "%s." % type(actx)) + + # validate meshmode_connection_kwargs + if meshmode_connection_kwargs is None: + meshmode_connection_kwargs = {} + if not isinstance(meshmode_connection_kwargs, dict): + raise TypeError("meshmode_connection_kwargs must be *None* or of " + "type dict, not '%s'." % + type(meshmode_connection_kwargs)) + if 'restrict_to_boundary' in meshmode_connection_kwargs: + raise ValueError("'restrict_to_boundary' must not be set by " + "meshmode_connection_kwargs") + + # build a meshmode connection for the source + src_bdy = None + if potential_source_and_target.get_source_dim() != tdim: + # sanity check + assert potential_source_and_target.get_source_dim() == tdim - 1 + src_bdy = potential_source_and_target.get_source_id() + + from meshmode.interop.firedrake import build_connection_from_firedrake + src_conn = build_connection_from_firedrake(actx, + self.dg_function_space, + restrict_to_boundary=src_bdy, + **meshmode_connection_kwargs) + # If source is a boundary, build a connection to restrict it to the + # boundary + restrict_conn = None + if src_bdy is not None: + # convert "everywhere" to meshmode BTAG_ALL + meshmode_src_bdy = src_bdy + if meshmode_src_bdy == "everywhere": + from meshmode.mesh import BTAG_ALL + meshmode_src_bdy = BTAG_ALL + # get group factory + order = self.dg_function_space.degree() + from meshmode.discretization.poly_element import \ + PolynomialRecursiveNodesGroupFactory + default_grp_factory = \ + PolynomialRecursiveNodesGroupFactory(order, 'lgl') + grp_factory = meshmode_connection_kwargs.get('grp_factory', + default_grp_factory) + from meshmode.discretization.connection import make_face_restriction + restrict_conn = make_face_restriction(actx, + src_conn.discr, + grp_factory, + meshmode_src_bdy) + + # Build a meshmode connection for the target + tgt_bdy = None + if potential_source_and_target.get_target_dim() != tdim: + # sanity check + assert potential_source_and_target.get_target_dim() == tdim - 1 + tgt_bdy = potential_source_and_target.get_target_id() + + # Can we re-use the source connection? + src_dim = potential_source_and_target.get_source_dim() + tgt_dim = potential_source_and_target.get_target_dim() + if tgt_bdy == src_bdy and src_dim == tgt_dim: + tgt_conn = src_conn + else: + # If not, build a new one + tgt_conn = build_connection_from_firedrake( + actx, + self.dg_function_space, + restrict_to_boundary=tgt_bdy, + **meshmode_connection_kwargs) + + # store computing context + self.actx = actx + # store connections + self.source_to_meshmode_connection = src_conn + self.restrict_source_to_boundary = restrict_conn + self.target_to_meshmode_connection = tgt_conn + + def from_firedrake(self, density): + # make sure we are in the dg function space + if not self.is_dg: + density = project(density, self.dg_function_space) + # Convert to meshmode. + density = \ + self.source_to_meshmode_connection.from_firedrake(density, + actx=self.actx) + # restrict to boundary if necessary + if self.restrict_source_to_boundary is not None: + density = self.restrict_source_to_boundary(density) + return density + + def to_firedrake(self, evaluated_potential, out=None): + # if we are dg, it's simple + if self.is_dg: + return self.target_to_meshmode_connection.from_meshmode( + evaluated_potential, + out=out) + else: + # Otherwise, we have to project back to our function space + pot = \ + self.target_to_meshmode_connection.from_meshmode(evaluated_potential) + pot = project(evaluated_potential, self.function_space) + if out is not None: + out.dat.data[:] = pot.dat.data[:] + pot = out + return pot + + def get_source_discretization(self): + """ + Get the :class:`meshmode.discretization.Discretization` + of the source + """ + if self.restrict_source_to_boundary is None: + return self.source_to_meshmode_connection.discr + return self.restrict_source_to_boundary.to_discr + + def get_target_discretization(self): + """ + Get the :class:`meshmode.discretization.Discretization` + of the target + """ + return self.target_to_meshmode_connection.discr + + +def PytentialOperation(actx, + density, + unbound_op, + density_name, + potential_src_and_tgt, + **kwargs): + """ + :arg actx: A :class:`meshmode.dof_array.PyOpenCLArrayContext` + :arg density: the :mod:`firedrake` density function + :arg unbound_op: A :mod:`pytential` operation which has not been + bound (e.g. a :mod:`pymbolic` expression) + :arg density_name: A string, the name of the density function + in the unbound_op + :arg potential_src_and_tgt: A :class:`PotentialSourceAndTarget` + + :kwarg warn_if_cg: Pass *False* to suppress the warning if + a "CG" space is used + :kwarg meshmode_connection_kwargs: *None*, or + a dict with arguments to pass to + :func:`meshmode.interop.firedrake.build_connection_from_firedrake` + :kwarg qbx_kwargs: *None*, or a dict with arguments to pass to + :class:`pytential.qbx.QBXLayerPotentialSource`. + :kwarg op_kwargs: kwargs passed to the invocation of the + bound pytential operator (e.g. {'k': 0.5} for + a :class:`sumpy.kernel.HelmholtzKernel`). + Must not include *density_name* + """ + # make sure density name is a string + if not isinstance(density_name, str): + raise TypeError("density_name must be of type str, not '%s'." % + density_name) + + # get kwargs to build connection + warn_if_cg = kwargs.get('warn_if_cg', True) + meshmode_connection_kwargs = kwargs.get('meshmode_connection_kwargs', None) + + # Build meshmode connection + meshmode_connection = MeshmodeConnection( + density.function_space(), + potential_src_and_tgt, + actx, + warn_if_cg=warn_if_cg, + meshmode_connection_kwargs=meshmode_connection_kwargs) + + # Build QBX + src_discr = meshmode_connection.get_source_discretization() + qbx_kwargs = kwargs.get('qbx_kwargs', None) + from pytential.qbx import QBXLayerPotentialSource + qbx = QBXLayerPotentialSource(src_discr, **qbx_kwargs) + + # Bind pytential operator + tgt_discr = meshmode_connection.get_target_discretization() + from pytential import bind + pyt_op = bind((qbx, tgt_discr), unbound_op) + + # Get operator kwargs + op_kwargs = kwargs.get('op_kwargs', {}) + if density_name in op_kwargs: + raise ValueError(f"density_name '{density_name}' should not be included" + " in op_kwargs.") + + # build bound operator that just takes density as argument + def bound_op_with_kwargs(density_arg): + op_kwargs[density_name] = density_arg + return pyt_op(**op_kwargs) + + # Now build and return Potential object + return Potential(density, + connection=meshmode_connection, + potential_operator=bound_op_with_kwargs) diff --git a/firedrake/potential_evaluation/volumential.py b/firedrake/potential_evaluation/volumential.py new file mode 100644 index 0000000000..829d375f80 --- /dev/null +++ b/firedrake/potential_evaluation/volumential.py @@ -0,0 +1,568 @@ +from firedrake.potential_evaluation import PotentialEvaluationLibraryConnection +from firedrake.potential_evaluation.pytential import \ + MeshmodeConnection + + +class VolumentialConnection(PotentialEvaluationLibraryConnection): + """ + Build a connection into :mod:`volumential` + """ + def __init__(self, + potential_source_and_target, + function_space, + cl_ctx, + queue, + box_fmm_geometry_kwargs, + reorder_potentials, + meshmode_connection_kwargs, + warn_if_cg=True): + """ + For other args see :class:`PotentialEvaluationLibraryConnection`. + The *potential_source_and_target* must have a 3D + source, 3D target, and 3D mesh of co-dimension 0. + + :arg cl_ctx: A :class:`pyopencl.Context` + :arg queue: A :class:`pyopencl.CommandQueue` + :arg box_fmm_geometry_kwargs: kwargs passed to + :class:`volumential.geometry.BoxFMMGeometryFactory`. + Must include the key *'nlevels'* + :arg reorder_potentials: *True* iff the volumential FMM + is reordering potentials. *False* otherwise. + :arg meshmode_connection_kwargs: A dict passes as kwargs + to :func:`meshmode.interop.build_connection_from_firedrake`. + :mod:`meshmode` is used an intermediate level between + :mod:`firedrake` and :mod:`volumential`. + """ + # super + PotentialEvaluationLibraryConnection.__init__( + self, + function_space, + potential_source_and_target, + warn_if_cg=warn_if_cg) + + # validate that sources, target, and mesh are of right dimension + if potential_source_and_target.get_source_dimension() != 3: + raise ValueError("potential source must have dimension 3") + if potential_source_and_target.get_target_dimension() != 3: + raise ValueError("potential target must have dimension 3") + if potential_source_and_target.mesh().geometric_dimension() != 3: + raise ValueError("mesh must have geometric dimension 3") + if potential_source_and_target.mesh().topological_dimension() != 3: + raise ValueError("mesh must have topological dimension 3") + + # FIXME : allow subdomain regions + if potential_source_and_target.get_source_id() != "everywhere": + raise NotImplementedError("subdomain sources not implemented") + if potential_source_and_target.get_target_id() != "everywhere": + raise NotImplementedError("subdomain targets not implemented") + + # validate cl_ctx argument + from pyopencl import Context + if not isinstance(cl_ctx, Context): + raise TypeError("cl_ctx must be of type Context, not " + "%s." % type(cl_ctx)) + + # validate queue argument + from pyopencl import CommandQueue + if not isinstance(queue, CommandQueue): + raise TypeError("queue must be of type CommandQueue, not " + "%s." % type(queue)) + + # validate reorder_potentials + if not isinstance(reorder_potentials, bool): + raise TypeError("reorder_potentials must be of type bool, not " + "%s." % type(reorder_potentials)) + + # validate box_fmm_geometry_kwargs + if not isinstance(box_fmm_geometry_kwargs, dict): + raise TypeError("box_fmm_geometry_kwargs must be of " + "type dict, not '%s'." % + type(box_fmm_geometry_kwargs)) + if 'nlevels' not in box_fmm_geometry_kwargs: + raise ValueError("box_fmm_geometry_kwargs missing required keyword" + " 'nlevels'.") + + # validate meshmode_connection_kwargs + if meshmode_connection_kwargs is None: + meshmode_connection_kwargs = {} + if not isinstance(meshmode_connection_kwargs, dict): + raise TypeError("meshmode_connection_kwargs must be *None* or of " + "type dict, not '%s'." % + type(meshmode_connection_kwargs)) + + # Build intermediate connection into meshmode + from meshmode.interop.firedrake import build_connection_from_firedrake + from meshmode.array_context import PyOpenCLArrayContext + actx = PyOpenCLArrayContext(queue) + meshmode_connection = MeshmodeConnection( + function_space, + potential_source_and_target, + actx, + warn_if_cg=warn_if_cg, + meshmode_connection_kwargs=meshmode_connection_kwargs) + + # Build connection from meshmode into volumential + # (following https://gitlab.tiker.net/xywei/volumential/-/blob/fe2c3e7af355d5c527060e783237c124c95397b5/test/test_interpolation.py#L72 ) # noqa : E501 + from volumential.geometry import ( + BoundingBoxFactory, BoxFMMGeometryFactory) + from volumential.interpolation import ElementsToSourcesLookupBuilder + dim = function_space.mesh().geometric_dimension() + bboxfmm_fac = BoxFMMGeometryFactory(cl_ctx, dim=dim, + **box_fmm_geometry_kwargs) + boxgeo = bboxfmm_fac(queue) + elt_to_src_lookup_fac = ElementsToSourcesLookupBuilder( + cl_ctx, tree=boxgeo.tree, discr=meshmode_connection.discr) + elt_to_src_lookup, evt = elt_to_src_lookup_fac(queue) + + # Build connection from volumential into meshmode + from volumential.interpolation import LeavesToNodesLookupBuilder + leaves_to_node_lookup_fac = LeavesToNodesLookupBuilder( + cl_ctx, trav=boxgeo.trav, discr=meshmode_connection.discr) + leaves_to_node_lookup, evt = leaves_to_node_lookup_fac(queue) + + # Figure out whether or not conversion from volumential -> meshmode + # should user order or tree order + order = "user" if reorder_potentials else "tree" + + # Store maps for conversion to/from volumential + self.meshmode_connection = meshmode_connection + self.to_volumential_lookup = elt_to_src_lookup + self.from_volumential_lookup = leaves_to_node_lookup + self.from_volumential_order = order + # store necessary computing contextx + self.queue = queue + self.actx = actx + + def from_firedrake(self, density): + # pass operand into a meshmode DOFArray + from meshmode.dof_array import flatten + meshmode_density = self.meshmode_connection.from_firedrake(density) + meshmode_density = flatten(meshmode_density) + + # pass flattened operand into volumential interpolator + from volumential.interpolation import interpolate_from_meshmode + volumential_density = \ + interpolate_from_meshmode(self.queue, + meshmode_density, + self.to_volumential_lookup, + order="user") # user order is more intuitive + return volumential_density + + def to_firedrake(self, evaluated_potential, out=None): + # pass volumential back to meshmode DOFArray + from volumential.interpolation import interpolate_to_meshmode + from meshmode.dof_array import unflatten + meshmode_pot_vals = \ + interpolate_to_meshmode(self.queue, + evaluated_potential, + self.from_volumential_lookup, + order=self.from_volumential_order) + meshmode_pot_vals = unflatten(self.actx, + self.meshmode_connection.discr, + meshmode_pot_vals) + # get meshmode data back as firedrake function + return self.meshmode_connection.to_firedrake(meshmode_pot_vals, out=self) + +def VolumentialOperation(density, unbound_op, potential_src_and_tgt, + function_space=None): + pass + + +# TODO: Merge this into VolumentialOperation function +class VolumePotential(AbstractExternalOperator): + r""" + Evaluates to + + .. math:: + + f(x) = \int_\Om K(x-y) op(y) \,dx + + where op is a function-type living in *function_space*, \Om is the + 3-D mesh on which *function_space* lives, and K + a :mod:`sumpy` kernel. + + :arg operator_data: A map with required keys + + * 'kernel': a :class:`sumpy.kernel.Kernel` (*K* in the math above) + * 'kernel_type': A string describing the kernel, e.g. ``"Laplace"`` + * 'cl_ctx': a :mod:`pyopencl` computing context + * 'queue': a :mod:`pyopencl` command queue compatible with the + computing context + * 'nlevels': The desired number of levels to be used by + :mod:`volumential` + * 'm_order': multipole order used by the expansion + * 'dataset_filename': the filename to pass to + :class:`volumential.table_manager.NearFieldInteractionTableManager` + + And optional keys + + * 'grp_factory': (optional) An interpolatory group factory + inheriting from :class:`meshmode.discretization.ElementGroupFactory` + to be used in the intermediate :mod:`meshmode` representation + * 'q_order': (optional) The desired :mod:`volumential` quadrature + order, defaults to *function_space*'s degree + * 'force_direct_evaluation': (optional) As in + :func:`volumential.volume_fmm.drive_volume_fmm`. + Defaults to *False* + * 'fmm_kwargs': (optional) A dictionary of kwargs + to pass to :func:`volumential.volume_fmm.drive_fmm`. + * 'root_extent': (optional) the root extent to pass to + :class:`volumential.table_manager.NearFieldInteractionTableManager` + (defaults to 1) + * 'table_compute_method': (optional) used by + :class:`volumential.table_manager.NearFieldInteractionTableManager` + (defaults to "DrosteSum") + * 'table_kwargs': (optional) passed to + :class:`volumential.table_manager.NearFieldInteractionTableManager` + 's *get_table* method + + """ + + _external_operator_type = 'GLOBAL' + + def __init__(self, *operands, **kwargs): + self.operands = operands + self.kwargs = kwargs + AbstractExternalOperator.__init__(self, operands[0], **kwargs) + + + # Build a broken version of input function space + fs = kwargs["function_space"] + fsdg = FunctionSpace(fs.mesh(), "DG", fs.ufl_element().degree()) + + self.evaluator_kwargs = kwargs.copy() + self.evaluator_kwargs["function_space"] = fsdg + + # Validate input + # print(len(self.operands), self.derivatives) + # assert self.derivatives == (0,), \ + # "Derivatives of volume potential not currently supported: " + str(self.derivatives) + # from firedrake import Function + # + # if not isinstance(operand, Function): + # raise TypeError(":arg:`operand` must be of type firedrake.Function" + # ", not %s." % type(operand)) + # if operand.function_space().shape != tuple(): + # raise ValueError(":arg:`operand` must be a function with shape ()," + # " not %s." % operand.function_space().shape) + operator_data = kwargs["operator_data"] + assert isinstance(operator_data, dict) + required_keys = ('kernel', 'kernel_type', 'cl_ctx', 'queue', 'nlevels', + 'm_order', 'dataset_filename') + optional_keys = ('grp_factory', 'q_order', 'force_direct_evaluation', + 'fmm_kwargs', 'root_extent', + 'table_compute_method', 'table_kwargs', + ) + permissible_keys = required_keys + optional_keys + if not all(key in operator_data for key in required_keys): + raise ValueError("operator_data is missing one of the required " + "keys: %s" % required_keys) + if not all(key in permissible_keys for key in operator_data): + raise ValueError("operator_data contains an unexpected key. All " + "keys must be one of %s." % (permissible_keys,)) + + @cached_property + def _evaluator(self): + return VolumePotentialEvaluator(self, *self.operands, **self.evaluator_kwargs) + + def _evaluate(self, continuity_tolerance=None): + 1/0 + + def _compute_derivatives(self, continuity_tolerance=None): + 1/0 + + def _evaluate_action(self, *args, continuity_tolerance=None): + return self._evaluator._evaluate_action(continuity_tolerance=continuity_tolerance) + + def _evaluate_adjoint_action(self, *args, continuity_tolerance=None): + return self._evaluator._evaluate_action(continuity_tolerance=continuity_tolerance) + + def evaluate_adj_component_control(self, x, idx): + derivatives = tuple(dj + int(idx == j) for j, dj in enumerate(self.derivatives)) + dN = self._ufl_expr_reconstruct_(*self.ufl_operands, derivatives=derivatives) + return dN._evaluate_adjoint_action((x.function,)).vector() + + def evaluate_adj_component(self, x, idx): + print(x, type(x)) + 1/0 + + +class VolumePotentialEvaluator: + def __init__(self, vp, *operands, **kwargs): + self.vp = vp + operand = operands[0] + + cg_function_space = vp.function_space() + function_space = kwargs["function_space"] + self.dg_function_space = function_space + + # Create an interpolator from the original operand so that + # we can re-interpolate it into the space at each evaluation + # in case it has changed (e.g. cost functional in PDE constrained optimization. + + operator_data = kwargs["operator_data"] + kernel = operator_data['kernel'] + kernel_type = operator_data['kernel_type'] + cl_ctx = operator_data['cl_ctx'] + queue = operator_data['queue'] + nlevels = operator_data['nlevels'] + m_order = operator_data['m_order'] + dataset_filename = operator_data['dataset_filename'] + degree = function_space.ufl_element().degree() + grp_factory = operator_data.get('grp_factory', None) + q_order = operator_data.get('q_order', degree) + force_direct_evaluation = operator_data.get('force_direct_evaluation', + False) + fmm_kwargs = operator_data.get('fmm_kwargs', {}) + root_extent = operator_data.get("root_extent", 1) + table_compute_method = operator_data.get('table_compute_method', + 'DrosteSum') + table_kwargs = operator_data.get('table_kwargs', {}) + + from sumpy.kernel import Kernel + if not isinstance(kernel, Kernel): + raise TypeError("operator_data['kernel'] must be of type " + "sumpy.kernel.Kernel, not %s." % type(kernel)) + if not isinstance(nlevels, int): + raise TypeError("operator_data['nlevels'] must be of type int, " + "not %s." % type(nlevels)) + if not isinstance(m_order, int): + raise TypeError("operator_data['m_order'] must be of type int, " + "not %s." % type(m_order)) + if not isinstance(q_order, int): + raise TypeError("operator_data['q_order'] must be of type int, " + "not %s." % type(q_order)) + + # Build connection into meshmode + from meshmode.interop.firedrake import build_connection_from_firedrake + from meshmode.array_context import PyOpenCLArrayContext + actx = PyOpenCLArrayContext(queue) + meshmode_connection = build_connection_from_firedrake( + actx, function_space, grp_factory=grp_factory) + + # Build connection from meshmode into volumential + # (following https://gitlab.tiker.net/xywei/volumential/-/blob/fe2c3e7af355d5c527060e783237c124c95397b5/test/test_interpolation.py#L72 ) # noqa : E501 + from volumential.geometry import ( + BoundingBoxFactory, BoxFMMGeometryFactory) + from volumential.interpolation import ElementsToSourcesLookupBuilder + dim = function_space.mesh().geometric_dimension() + bbox_fac = BoundingBoxFactory(dim=dim) + bboxfmm_fac = BoxFMMGeometryFactory( + cl_ctx, dim=dim, order=q_order, + nlevels=nlevels, bbox_getter=bbox_fac, + expand_to_hold_mesh=meshmode_connection.discr.mesh, + mesh_padding_factor=0.0) + boxgeo = bboxfmm_fac(queue) + elt_to_src_lookup_fac = ElementsToSourcesLookupBuilder( + cl_ctx, tree=boxgeo.tree, discr=meshmode_connection.discr) + elt_to_src_lookup, evt = elt_to_src_lookup_fac(queue) + + # Build connection from volumential into meshmode + from volumential.interpolation import LeavesToNodesLookupBuilder + leaves_to_node_lookup_fac = LeavesToNodesLookupBuilder( + cl_ctx, trav=boxgeo.trav, discr=meshmode_connection.discr) + leaves_to_node_lookup, evt = leaves_to_node_lookup_fac(queue) + + # Create near-field table in volumential + from volumential.table_manager import NearFieldInteractionTableManager + table_manager = NearFieldInteractionTableManager( + dataset_filename, root_extent=root_extent, queue=queue) + + nftable, _ = table_manager.get_table( + dim, + kernel_type, + q_order, + compute_method=table_compute_method, + queue=queue, + **table_kwargs) + + # Create kernel wrangler in volumential + from sumpy.expansion import DefaultExpansionFactory + out_kernels = [kernel] + expn_factory = DefaultExpansionFactory() + local_expn_class = expn_factory.get_local_expansion_class(kernel) + mpole_expn_class = expn_factory.get_multipole_expansion_class(kernel) + + from functools import partial + from volumential.expansion_wrangler_fpnd import ( + FPNDExpansionWrangler, FPNDExpansionWranglerCodeContainer) + wrangler_cc = FPNDExpansionWranglerCodeContainer( + cl_ctx, + partial(mpole_expn_class, kernel), + partial(local_expn_class, kernel), + out_kernels) + wrangler = FPNDExpansionWrangler( + code_container=wrangler_cc, + queue=queue, + tree=boxgeo.tree, + near_field_table=nftable, + dtype=ScalarType, + fmm_level_to_order=lambda kernel, kernel_args, tree, lev: m_order, + quad_order=q_order) + + # Store attributes that we may need later + # self.ufl_operands = operands + self.actx = actx + self.queue = queue + self.sumpy_kernel = kernel + self.meshmode_connection = meshmode_connection + self.volumential_boxgeo = boxgeo + self.to_volumential_lookup = elt_to_src_lookup + self.from_volumential_lookup = leaves_to_node_lookup + self.force_direct_evaluation = force_direct_evaluation + self.fmm_kwargs = fmm_kwargs + self.expansion_wrangler = wrangler + # so we don't have to keep making a fd function during conversion + # from meshmode. AbstractExternalOperator s aren't functions, + # so can't use *self* + self.fd_pot = Function(self.dg_function_space) + self.fd_pot_cg = Function(vp.function_space()) + + def _evaluate(self, continuity_tolerance=None): + """ + :arg continuity_tolerance: If *None* then ignored. Otherwise, a floating + point value. If the function space associated to this operator + is a continuous one and *continuity_tolerance* is a float, it is + verified that, + during the conversion of the evaluated potential + from meshmode's discontinuous function space representation into + firedrake's continuous one, the potential's value at + a firedrake node is within *continuity_tolerance* of its + value at any duplicated firedrake node + """ + operand, = self.vp.ufl_operands + operand_discrete = interpolate(operand, self.dg_function_space) + + # pass operand into a meshmode DOFArray + from meshmode.dof_array import flatten + meshmode_src_vals = \ + self.meshmode_connection.from_firedrake(operand_discrete, actx=self.actx) + meshmode_src_vals = flatten(meshmode_src_vals) + + # if 1: + # alpha = 40 + # dim = 2 + + # if 1: + # mm_discr = self.meshmode_connection.discr + # from meshmode.dof_array import thaw, flat_norm + # x = flatten(thaw(self.actx, mm_discr.nodes()[0])) + # y = flatten(thaw(self.actx, mm_discr.nodes()[1])) + # import pyopencl.clmath as clmath + # import pyopencl.array as cla + # norm2 = (x-0.5)**2 + (y-0.5)**2 + # ref_src = -(4 * alpha ** 2 * norm2 - 2 * dim * alpha) * clmath.exp( + # -alpha * norm2) + # denom = np.max(np.abs(ref_src.get())) + # err = np.abs((ref_src - meshmode_src_vals).get()) + # print("Fdrake -> meshmode error:", np.max(err)/denom) + + # pass flattened operand into volumential interpolator + from volumential.interpolation import interpolate_from_meshmode + volumential_src_vals = \ + interpolate_from_meshmode(self.queue, + meshmode_src_vals, + self.to_volumential_lookup, + order="user") # user order is more intuitive + + # if 1: + # # check source density against known value + # # x = self.volumential_boxgeo.tree.sources[0].with_queue(self.queue) + # # y = self.volumential_boxgeo.tree.sources[1].with_queue(self.queue) + # x, y = self.volumential_boxgeo.nodes + # import pyopencl.clmath as clmath + # import pyopencl.array as cla + # norm2 = (x-0.5)**2 + (y-0.5)**2 + # ref_src = -(4 * alpha ** 2 * norm2 - 2 * dim * alpha) * clmath.exp( + # -alpha * norm2) + # denom = np.max(np.abs(ref_src.get())) + # err = np.abs((ref_src-volumential_src_vals).get()) + # print("Fdrake -> volumential error:", np.max(err)/denom) + + # evaluate volume fmm + from volumential.volume_fmm import drive_volume_fmm + pot, = drive_volume_fmm( + self.volumential_boxgeo.traversal, + self.expansion_wrangler, + volumential_src_vals * self.volumential_boxgeo.weights, + volumential_src_vals, + direct_evaluation=self.force_direct_evaluation, + reorder_sources=True, + **self.fmm_kwargs) + + # if 1: + # # check potential against known value + # ref_pot = clmath.exp(-alpha * norm2) + # denom = np.max(np.abs(ref_pot.get())) + # err = np.abs((ref_pot-pot).get()) + # print("volumential potential error:", np.max(err)/denom) + + # # evaluate volume fmm with ref_src + # pot_w_ref_src, = drive_volume_fmm( + # self.volumential_boxgeo.traversal, + # self.expansion_wrangler, + # ref_src * self.volumential_boxgeo.weights, + # ref_src, + # direct_evaluation=self.force_direct_evaluation, + # reorder_sources=True, + # **self.fmm_kwargs) + + # err = np.abs((ref_pot-pot_w_ref_src).get()) + # print("volumential potential w/t ref src error:", np.max(err)/denom) + + # pass volumential back to meshmode DOFArray + from volumential.interpolation import interpolate_to_meshmode + from meshmode.dof_array import unflatten + user_order = self.fmm_kwargs.get("reorder_potentials", True) + if user_order: + order = "user" + else: + order = "tree" + meshmode_pot_vals = interpolate_to_meshmode(self.queue, + pot, + self.from_volumential_lookup, + order=order) + meshmode_pot_vals = unflatten(self.actx, + self.meshmode_connection.discr, + meshmode_pot_vals) + # get meshmode data back as firedrake function + self.meshmode_connection.from_meshmode( + meshmode_pot_vals, + out=self.fd_pot) + + self.fd_pot_cg = project(self.fd_pot, self.vp.function_space()) + + self.vp.dat.data[:] = self.fd_pot_cg.dat.data[:] + + return self.fd_pot_cg + + def _compute_derivatives(self, continuity_tolerance=None): + # TODO : Support derivatives + return self._evaluate(continuity_tolerance=continuity_tolerance) + + def _evaluate_action(self, *args, continuity_tolerance=None): + # From tests/pointwiseoperator/test_point_expr.py + # https://github.com/firedrakeproject/firedrake/blob/c0d9b592f587fa8c7437f690da7a6595f6804c1b/tests/pointwiseoperator/test_point_expr.py # noqa + if len(args) == 0: + # Evaluate the operator + return self._evaluate(continuity_tolerance=continuity_tolerance) + + # Evaluate the Jacobian/Hessian action + operands = self.ufl_operands + operator = self._compute_derivatives(continuity_tolerance=continuity_tolerance) + expr = as_ufl(operator(*operands)) + if expr.ufl_shape == () and expr != 0: + var = VariableRuleset(self.ufl_operands[0]) + expr = expr*var._Id + elif expr == 0: + return self.assign(expr) + + for arg in args: + mi = indices(len(expr.ufl_shape)) + aa = mi + bb = mi[-len(arg.ufl_shape):] + expr = arg[bb] * expr[aa] + mi_tensor = tuple(e for e in mi if not (e in aa and e in bb)) + if len(expr.ufl_free_indices): + expr = as_tensor(expr, mi_tensor) + return self.interpolate(expr) From 64bf5b58d10cd0bc0b55da0bdda3c2b6487d8a7c Mon Sep 17 00:00:00 2001 From: benSepanski Date: Thu, 18 Feb 2021 16:10:48 -0600 Subject: [PATCH 04/20] fix imports and begin fixing test --- firedrake/__init__.py | 2 +- .../test_layer_potentials.py | 75 +++++-------------- 2 files changed, 21 insertions(+), 56 deletions(-) diff --git a/firedrake/__init__.py b/firedrake/__init__.py index 5c4a6d3ab3..9561018191 100644 --- a/firedrake/__init__.py +++ b/firedrake/__init__.py @@ -89,7 +89,7 @@ from firedrake.ensemble import * from firedrake.randomfunctiongen import * from firedrake.pointwise_operators import * -from firedrake.layer_potentials import * +from firedrake.potential_evaluation import * from firedrake.logging import * # Set default log level diff --git a/tests/pointwiseoperator/test_layer_potentials.py b/tests/pointwiseoperator/test_layer_potentials.py index 5fa3978b02..81c2455a96 100644 --- a/tests/pointwiseoperator/test_layer_potentials.py +++ b/tests/pointwiseoperator/test_layer_potentials.py @@ -7,7 +7,8 @@ ln, pi, SpatialCoordinate, sqrt, PytentialLayerOperation, grad, \ FunctionSpace, VectorFunctionSpace, FacetNormal, inner, assemble, \ TestFunction, ds, Function, exp, TrialFunction, dx, project, solve, \ - utils, OpenCascadeMeshHierarchy, dot + utils, OpenCascadeMeshHierarchy, dot, SingleLayerPotential, + DoubleLayerPotential, PotentialSourceAndTarget from math import factorial from warnings import warn @@ -138,64 +139,28 @@ def get_true_sol_expr(spatial_coord): # Build function spaces cgfspace = FunctionSpace(mesh, "CG", fspace_degree) cgvfspace = VectorFunctionSpace(mesh, "CG", fspace_degree) - dgvfspace = VectorFunctionSpace(mesh, "DG", fspace_degree) - - # pytential normals point opposite direction of firedrake - pyt_inner_normal_sign = -1 - ambient_dim = mesh.geometric_dimension() - # Build rhs pytential operations - from pytential import sym - from sumpy.kernel import HelmholtzKernel - sigma = sym.make_sym_vector("density", ambient_dim) - r""" - ..math: - - x \in \Sigma - grad_op(x) = \nabla( \int_\Gamma( f(y) H_0^{(1)}(\kappa |x - y|) )d\gamma(y)) - """ - rhs_pyt_grad_layer = pyt_inner_normal_sign * \ - sym.grad(ambient_dim, sym.S(HelmholtzKernel(ambient_dim), - sym.n_dot(sigma), - k=sym.var("k"), qbx_forced_limit=None)) - r""" - ..math: - x \in \Sigma - op(x) = i \kappa \cdot \int_\Gamma( f(y) H_0^{(1)}(\kappa |x - y|) )d\gamma(y) - """ - rhs_pyt_layer = 1j * sym.var("k") * pyt_inner_normal_sign * \ - sym.S(HelmholtzKernel(ambient_dim), - sym.n_dot(sigma), - k=sym.var("k"), - qbx_forced_limit=None) - dg_true_sol_grad = Function(dgvfspace).interpolate(grad(true_sol_expr)) - # general operator data settings, missing 'op' - operator_data = {'actx': actx, - 'density_name': 'density', - 'source_bdy_id': source_bdy_id, - 'target_bdy_id': target_bdy_id, - 'op_kwargs': {'k': kappa}, - 'qbx_order': fspace_degree + 2, - 'fine_order': 4 * fspace_degree, - 'fmm_order': 50, - } - # Bind rhs pytential operations into firedrake - rhs_grad_layer_operator_data = dict(operator_data) - rhs_grad_layer_operator_data['op'] = rhs_pyt_grad_layer - rhs_grad_layer = PytentialLayerOperation(dg_true_sol_grad, - function_space=cgvfspace, - operator_data=rhs_grad_layer_operator_data) - rhs_layer_operator_data = dict(operator_data) - rhs_layer_operator_data['op'] = rhs_pyt_layer - rhs_layer = PytentialLayerOperation(dg_true_sol_grad, - function_space=cgfspace, - operator_data=rhs_layer_operator_data) + # TODO: Fix RHS + # TODO: Fix true solution expr + true_sol = Function(cgfspace).interpolate(true_sol_expr) + + # FIXME: Handle normal signs + places = PotentialSourceAndTarget(mesh, + source_region_dim=2, + source_region_id=source_bdy_id, + target_region_dim=2, + target_region_id=target_bdy_id) + + Dtrue_sol = DoubleLayerPotential(dot(true_sol, n), actx=actx) + Strue_sol = SingleLayerPotential(dot(true_sol, n), actx=actx) # get rhs form v = TestFunction(cgfspace) - rhs_form = inner(dot(grad(true_sol_expr), FacetNormal(mesh)), + rhs_form = inner(dot(grad(true_sol), n), v) * ds(source_bdy_id, metadata={'quadrature_degree': 2 * fspace_degree}) \ - + inner(rhs_layer, v) * ds(target_bdy_id) \ - - inner(dot(rhs_grad_layer, FacetNormal(mesh)), v) * ds(target_bdy_id) + - 1j * inner(Strue_sol, v) * ds(target_bdy_id) \ + - inner(dot(grad(Dtrue_sol), n), v) * ds(target_bdy_id) + + # TODO: Continue fixing test # local helmholtz operator r""" From d3b5800eb7f77fb85e329320bc87438a37b741d1 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Thu, 18 Feb 2021 16:38:24 -0600 Subject: [PATCH 05/20] fixed circular import --- firedrake/potential_evaluation/potentials.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/potential_evaluation/potentials.py b/firedrake/potential_evaluation/potentials.py index 275dffb710..1410729688 100644 --- a/firedrake/potential_evaluation/potentials.py +++ b/firedrake/potential_evaluation/potentials.py @@ -1,5 +1,4 @@ from firedrake.functionspaceimpl import WithGeometry -from firedrake.potential_evaluation.pytential import PytentialOperation def SingleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): @@ -84,6 +83,7 @@ def _layer_potential(layer_potential_sym, # extract Pytential operation kwargs pyt_op_kwargs = _extract_pytential_operation_kwargs(**kwargs) # now return the pytential operation as an external operator + from firedrake.potential_evaluation.pytential import PytentialOperation return PytentialOperation(actx, density, unbound_op, "density", potential_src_and_tgt, **pyt_op_kwargs) From 36741173dc9743c94f3f95f14662fa88b949ab00 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Thu, 18 Feb 2021 16:47:57 -0600 Subject: [PATCH 06/20] Make sure to ask for function space argument if density is not a function --- firedrake/potential_evaluation/potentials.py | 35 +++++++++++++++++--- firedrake/potential_evaluation/pytential.py | 13 +++++--- 2 files changed, 38 insertions(+), 10 deletions(-) diff --git a/firedrake/potential_evaluation/potentials.py b/firedrake/potential_evaluation/potentials.py index 1410729688..b0577b9665 100644 --- a/firedrake/potential_evaluation/potentials.py +++ b/firedrake/potential_evaluation/potentials.py @@ -12,6 +12,10 @@ def SingleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): :arg potential_src_and_tgt: A :class:`firedrake.potential_evaluation.PotentialSourceAndTarget` + :kwarg function_space: If density is a firedrake function, + this is optional. Otherwise it must be supplied. + If supplied and density is a function, they must + be the same object. :kwarg cl_ctx: (Optional) a :class:`pyopencl.Context` used to create a command queue. At most one of *cl_ctx*, *queue*, *actx* should be included. @@ -77,15 +81,20 @@ def _layer_potential(layer_potential_sym, unbound_op = layer_potential_sym(kernel, sym.var("density"), qbx_forced_limit=qbx_forced_limit) - # make sure we got an actx + # make sure we got an actx during validation assert 'actx' in kwargs actx = kwargs['actx'] + # make sure we got a function_space during validation + assert 'function_space' in kwargs + function_space = kwargs['function_space'] # extract Pytential operation kwargs pyt_op_kwargs = _extract_pytential_operation_kwargs(**kwargs) # now return the pytential operation as an external operator from firedrake.potential_evaluation.pytential import PytentialOperation return PytentialOperation(actx, density, unbound_op, "density", - potential_src_and_tgt, **pyt_op_kwargs) + potential_src_and_tgt, + function_space=function_space, + **pyt_op_kwargs) def _validate_layer_potential_args_and_set_defaults(density, @@ -98,9 +107,25 @@ def _validate_layer_potential_args_and_set_defaults(density, Returns dictionary updated with default arguments """ # validate density function space - if not isinstance(density.function_space(), WithGeometry): - raise TypeError("density.function_space() must be of type " - f"WithGeometry, not {type(density.function_space())}") + if hasattr(density, 'function_space'): + if not isinstance(density.function_space(), WithGeometry): + raise TypeError("density.function_space() must be of type " + f"WithGeometry, not {type(density.function_space())}") + function_space = kwargs.get('function_space', None) + if function_space is not None: + if function_space is not density.function_space(): + raise ValueError("density.function_space() and function_space" + " must be the same object") + else: + kwargs['function_space'] = density.function_space() + else: + function_space = kwargs.get('function_space', None) + if function_space is None: + raise ValueError("density has no function_space method, so" + "function_space kwarg must be supplied") + if not isinstance(function_space, WithGeometry): + raise TypeError("function_space must be of type " + f"WithGeometry, not {type(function_space)}") # validate kernel from sumpy.kernel import Kernel diff --git a/firedrake/potential_evaluation/pytential.py b/firedrake/potential_evaluation/pytential.py index ad6b14b291..86e8a70645 100644 --- a/firedrake/potential_evaluation/pytential.py +++ b/firedrake/potential_evaluation/pytential.py @@ -187,6 +187,8 @@ def PytentialOperation(actx, bound pytential operator (e.g. {'k': 0.5} for a :class:`sumpy.kernel.HelmholtzKernel`). Must not include *density_name* + + Remaining kwargs are passed to Potential.__init__ """ # make sure density name is a string if not isinstance(density_name, str): @@ -194,8 +196,8 @@ def PytentialOperation(actx, density_name) # get kwargs to build connection - warn_if_cg = kwargs.get('warn_if_cg', True) - meshmode_connection_kwargs = kwargs.get('meshmode_connection_kwargs', None) + warn_if_cg = kwargs.pop('warn_if_cg', True) + meshmode_connection_kwargs = kwargs.pop('meshmode_connection_kwargs', None) # Build meshmode connection meshmode_connection = MeshmodeConnection( @@ -207,7 +209,7 @@ def PytentialOperation(actx, # Build QBX src_discr = meshmode_connection.get_source_discretization() - qbx_kwargs = kwargs.get('qbx_kwargs', None) + qbx_kwargs = kwargs.pop('qbx_kwargs', None) from pytential.qbx import QBXLayerPotentialSource qbx = QBXLayerPotentialSource(src_discr, **qbx_kwargs) @@ -217,7 +219,7 @@ def PytentialOperation(actx, pyt_op = bind((qbx, tgt_discr), unbound_op) # Get operator kwargs - op_kwargs = kwargs.get('op_kwargs', {}) + op_kwargs = kwargs.pop('op_kwargs', {}) if density_name in op_kwargs: raise ValueError(f"density_name '{density_name}' should not be included" " in op_kwargs.") @@ -230,4 +232,5 @@ def bound_op_with_kwargs(density_arg): # Now build and return Potential object return Potential(density, connection=meshmode_connection, - potential_operator=bound_op_with_kwargs) + potential_operator=bound_op_with_kwargs, + **kwargs) From 960475e441a367e3f7c65eb8f6f262636890ca0e Mon Sep 17 00:00:00 2001 From: benSepanski Date: Thu, 18 Feb 2021 19:33:53 -0600 Subject: [PATCH 07/20] Chasing down bugs --- firedrake/potential_evaluation/__init__.py | 30 +++------ firedrake/potential_evaluation/potentials.py | 37 ++++++----- firedrake/potential_evaluation/pytential.py | 34 ++++------ .../test_layer_potentials.py | 64 +++++++++++-------- 4 files changed, 81 insertions(+), 84 deletions(-) diff --git a/firedrake/potential_evaluation/__init__.py b/firedrake/potential_evaluation/__init__.py index 80d7465050..2ab4ff0241 100644 --- a/firedrake/potential_evaluation/__init__.py +++ b/firedrake/potential_evaluation/__init__.py @@ -49,7 +49,7 @@ class Potential(AbstractExternalOperator): """ _external_operator_type = 'GLOBAL' - def __init__(self, density, **kwargs): + def __init__(self, density, connection, potential_operator, **kwargs): """ :arg density: A :mod:`firedrake` :class:`firedrake.function.Function` or UFL expression which represents the density :math:`u` @@ -60,41 +60,23 @@ def __init__(self, density, **kwargs): # super AbstractExternalOperator.__init__(self, density, **kwargs) - # FIXME - for order in self.derivatives: - assert order == 1, "Assumes self.derivatives = (1,..,1)" - # Get connection & bound op and validate - connection = kwargs.get("connection", None) - if connection is None: - raise ValueError("Missing kwarg 'connection'") if not isinstance(connection, PotentialEvaluationLibraryConnection): raise TypeError("connection must be of type " "PotentialEvaluationLibraryConnection, not %s." % type(connection)) self.connection = connection - - self.potential_operator = kwargs.get("potential_operator", None) - if self.potential_operator is None: - raise ValueError("Missing kwarg 'potential_operator'") + self.potential_operator = potential_operator # Get function space and validate aginst bound op function_space = self.ufl_function_space() assert isinstance(function_space, WithGeometry) # sanity check - if function_space is not self.potential_operator.function_space(): - raise ValueError("function_space must be same object as " - "potential_operator.function_space().") # Make sure density is a member of our function space, if it is if isinstance(density, Function): if density.function_space() is not function_space: raise ValueError("density.function_space() must be the " "same as function_space") - # Make sure the shapes match, at least - elif density.shape != function_space.shape: - raise ValueError("Shape mismatch between function_space and " - "density. %s != %s." % - (density.shape, function_space.shape)) @cached_property def _evaluator(self): @@ -355,7 +337,7 @@ def __init__(self, mesh, continue # standardize id_ to a tuple if isinstance(id_, int): - id_ = tuple(id_) + id_ = tuple([id_]) # check type if not isinstance(id_, tuple): raise TypeError("source and target region ids must be " @@ -374,6 +356,12 @@ def __init__(self, mesh, # probably be doing anyway if we're targeting the region # # It also throws a ValueError if the id_ is invalid + # FIXME: Figure out how to support multiple ids + if len(id_) > 1: + raise NotImplementedError("Currently, any tuple of " + "cell subdomain ids must be of " + "length 1.") + id_, = id_ mesh.cell_subset(id_) # handle None case diff --git a/firedrake/potential_evaluation/potentials.py b/firedrake/potential_evaluation/potentials.py index b0577b9665..fabe683a02 100644 --- a/firedrake/potential_evaluation/potentials.py +++ b/firedrake/potential_evaluation/potentials.py @@ -28,16 +28,16 @@ def SingleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): used for :mod:`pytential` and :mod:`meshmode` computations. At most one of *cl_ctx*, *queue*, *actx* should be included. If none is included, a *cl_ctx* will be created. - :kwarg qbx_forced_limit: (Optional) As described in - https://documen.tician.de/pytential/_modules/pytential/symbolic/primitives.html#IntG - DEFAULT: *None* - :kwarg op_kwargs: (Optional) kwargs passed to the invocation of the - bound pytential operator (e.g. {'k': 0.5} for - a :class:`sumpy.kernel.HelmholtzKernel`). + :kwarg op_kwargs: (Optional) kwargs passed to the construction + of the unbound pytential operator (e.g. {'k': 0.5} for + a :class:`sumpy.kernel.HelmholtzKernel` must be passed + to :func:`pytential.symbolic.primitives.S`), + DEFAULT: + * {'qbx_forced_limit': None} :kwarg qbx_kwargs: (Optional) *None*, or a dict with arguments to pass to :class:`pytential.qbx.QBXLayerPotentialSource`. DEFAULTS: - - fine order: 4 8 function space degree + - fine_order: 4 8 function space degree - qbx_order: function space degree - fmm_order: *None* - fmm_level_to_order: A @@ -55,7 +55,7 @@ def SingleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): """ # TODO: Make fmm_level_to_order depend on number of derivatives somehow? from pytential.symbolic.primitives import S - _layer_potential(S, density, kernel, potential_src_and_tgt, **kwargs) + return _layer_potential(S, density, kernel, potential_src_and_tgt, **kwargs) def DoubleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): @@ -64,7 +64,7 @@ def DoubleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): instead of the single """ from pytential.symbolic.primitives import D - _layer_potential(D, density, kernel, potential_src_and_tgt, **kwargs) + return _layer_potential(D, density, kernel, potential_src_and_tgt, **kwargs) def _layer_potential(layer_potential_sym, @@ -75,12 +75,13 @@ def _layer_potential(layer_potential_sym, """ kwargs = _validate_layer_potential_args_and_set_defaults( density, kernel, potential_src_and_tgt, **kwargs) - qbx_forced_limit = kwargs.get('qbx_forced_limit', None) # build our unbound operator + assert 'op_kwargs' in kwargs # make sure we got op kwargs + op_kwargs = kwargs['op_kwargs'] from pytential import sym unbound_op = layer_potential_sym(kernel, sym.var("density"), - qbx_forced_limit=qbx_forced_limit) + **op_kwargs) # make sure we got an actx during validation assert 'actx' in kwargs actx = kwargs['actx'] @@ -93,7 +94,6 @@ def _layer_potential(layer_potential_sym, from firedrake.potential_evaluation.pytential import PytentialOperation return PytentialOperation(actx, density, unbound_op, "density", potential_src_and_tgt, - function_space=function_space, **pyt_op_kwargs) @@ -126,6 +126,7 @@ def _validate_layer_potential_args_and_set_defaults(density, if not isinstance(function_space, WithGeometry): raise TypeError("function_space must be of type " f"WithGeometry, not {type(function_space)}") + function_space = kwargs['function_space'] # validate kernel from sumpy.kernel import Kernel @@ -152,7 +153,7 @@ def _validate_layer_potential_args_and_set_defaults(density, # Make sure all kwargs are recognized allowed_kwargs = ("cl_ctx", "queue", "actx", "op_kwargs", "qbx_kwargs", - "meshmode_connection_kwargs") + "meshmode_connection_kwargs", "function_space") for key in kwargs: if key not in allowed_kwargs: raise ValueError(f"Unrecognized kwarg {key}") @@ -203,9 +204,9 @@ def _validate_layer_potential_args_and_set_defaults(density, # Set qbx_kwargs defaults if 'qbx_kwargs' not in kwargs: - degree = density.function_space().ufl_element().degree() + degree = function_space.ufl_element().degree() from sumpy.expansion.level_to_order import FMMLibExpansionOrderFinder - qbx_kwargs = {'fine order': 4 * degree, + qbx_kwargs = {'fine_order': 4 * degree, 'qbx_order': degree, 'fmm_order': None, 'fmm_level_to_order': @@ -214,6 +215,10 @@ def _validate_layer_potential_args_and_set_defaults(density, } kwargs['qbx_kwargs'] = qbx_kwargs + # Set op_kwargs defaults + if 'op_kwargs' not in kwargs: + kwargs['op_kwargs'] = {'qbx_forced_limit': None} + return kwargs @@ -225,7 +230,7 @@ def _extract_pytential_operation_kwargs(**kwargs): pyt_op_possible_keywords = ("warn_if_cg", "meshmode_connection_kwargs", "qbx_kwargs", - "op_kwargs") + "function_space") for key in pyt_op_possible_keywords: if key in kwargs: pyt_op_kwargs[key] = kwargs[key] diff --git a/firedrake/potential_evaluation/pytential.py b/firedrake/potential_evaluation/pytential.py index 86e8a70645..06affba5e7 100644 --- a/firedrake/potential_evaluation/pytential.py +++ b/firedrake/potential_evaluation/pytential.py @@ -29,10 +29,10 @@ def __init__(self, function_space, potential_source_and_target, actx, # FIXME : allow subdomain regions tdim = potential_source_and_target.mesh.topological_dimension() - if potential_source_and_target.get_source_dim() != tdim: + if potential_source_and_target.get_source_dimension() == tdim: if potential_source_and_target.get_source_id() != "everywhere": raise NotImplementedError("subdomain sources not implemented") - if potential_source_and_target.get_target_dim() != tdim: + if potential_source_and_target.get_target_dimension() == tdim: if potential_source_and_target.get_target_id() != "everywhere": raise NotImplementedError("subdomain targets not implemented") @@ -55,9 +55,9 @@ def __init__(self, function_space, potential_source_and_target, actx, # build a meshmode connection for the source src_bdy = None - if potential_source_and_target.get_source_dim() != tdim: + if potential_source_and_target.get_source_dimension() != tdim: # sanity check - assert potential_source_and_target.get_source_dim() == tdim - 1 + assert potential_source_and_target.get_source_dimension() == tdim - 1 src_bdy = potential_source_and_target.get_source_id() from meshmode.interop.firedrake import build_connection_from_firedrake @@ -75,7 +75,7 @@ def __init__(self, function_space, potential_source_and_target, actx, from meshmode.mesh import BTAG_ALL meshmode_src_bdy = BTAG_ALL # get group factory - order = self.dg_function_space.degree() + order = self.dg_function_space.ufl_element().degree() from meshmode.discretization.poly_element import \ PolynomialRecursiveNodesGroupFactory default_grp_factory = \ @@ -90,14 +90,14 @@ def __init__(self, function_space, potential_source_and_target, actx, # Build a meshmode connection for the target tgt_bdy = None - if potential_source_and_target.get_target_dim() != tdim: + if potential_source_and_target.get_target_dimension() != tdim: # sanity check - assert potential_source_and_target.get_target_dim() == tdim - 1 + assert potential_source_and_target.get_target_dimension() == tdim - 1 tgt_bdy = potential_source_and_target.get_target_id() # Can we re-use the source connection? - src_dim = potential_source_and_target.get_source_dim() - tgt_dim = potential_source_and_target.get_target_dim() + src_dim = potential_source_and_target.get_source_dimension() + tgt_dim = potential_source_and_target.get_target_dimension() if tgt_bdy == src_bdy and src_dim == tgt_dim: tgt_conn = src_conn else: @@ -176,6 +176,7 @@ def PytentialOperation(actx, in the unbound_op :arg potential_src_and_tgt: A :class:`PotentialSourceAndTarget` + :kwarg function_space: The function space of the operator :kwarg warn_if_cg: Pass *False* to suppress the warning if a "CG" space is used :kwarg meshmode_connection_kwargs: *None*, or @@ -183,10 +184,6 @@ def PytentialOperation(actx, :func:`meshmode.interop.firedrake.build_connection_from_firedrake` :kwarg qbx_kwargs: *None*, or a dict with arguments to pass to :class:`pytential.qbx.QBXLayerPotentialSource`. - :kwarg op_kwargs: kwargs passed to the invocation of the - bound pytential operator (e.g. {'k': 0.5} for - a :class:`sumpy.kernel.HelmholtzKernel`). - Must not include *density_name* Remaining kwargs are passed to Potential.__init__ """ @@ -201,7 +198,7 @@ def PytentialOperation(actx, # Build meshmode connection meshmode_connection = MeshmodeConnection( - density.function_space(), + kwargs['function_space'], potential_src_and_tgt, actx, warn_if_cg=warn_if_cg, @@ -218,16 +215,9 @@ def PytentialOperation(actx, from pytential import bind pyt_op = bind((qbx, tgt_discr), unbound_op) - # Get operator kwargs - op_kwargs = kwargs.pop('op_kwargs', {}) - if density_name in op_kwargs: - raise ValueError(f"density_name '{density_name}' should not be included" - " in op_kwargs.") - # build bound operator that just takes density as argument def bound_op_with_kwargs(density_arg): - op_kwargs[density_name] = density_arg - return pyt_op(**op_kwargs) + return pyt_op(**{density_name: density_arg}) # Now build and return Potential object return Potential(density, diff --git a/tests/pointwiseoperator/test_layer_potentials.py b/tests/pointwiseoperator/test_layer_potentials.py index 81c2455a96..f5e29dd271 100644 --- a/tests/pointwiseoperator/test_layer_potentials.py +++ b/tests/pointwiseoperator/test_layer_potentials.py @@ -4,10 +4,10 @@ import pytest from firedrake import MeshHierarchy, norms, Constant, \ - ln, pi, SpatialCoordinate, sqrt, PytentialLayerOperation, grad, \ + ln, pi, SpatialCoordinate, sqrt, grad, \ FunctionSpace, VectorFunctionSpace, FacetNormal, inner, assemble, \ TestFunction, ds, Function, exp, TrialFunction, dx, project, solve, \ - utils, OpenCascadeMeshHierarchy, dot, SingleLayerPotential, + utils, OpenCascadeMeshHierarchy, dot, SingleLayerPotential, \ DoubleLayerPotential, PotentialSourceAndTarget from math import factorial from warnings import warn @@ -128,37 +128,45 @@ def get_true_sol_expr(spatial_coord): levels=3, order=2, project_refinements_to_cad=False) - source_bdy_id = 5 # (inner boundary) the circle - target_bdy_id = (1, 2, 3, 4) # (outer boundary) the square + scatterer_bdy = 5 # (inner boundary) the circle + truncated_bdy = (1, 2, 3, 4) # (outer boundary) the square # Solve for each mesh in hierarchy for h, mesh in zip([0.5 * 2**i for i in range(len(mesh_hierarchy))], mesh_hierarchy): + # Build function spaces + V = FunctionSpace(mesh, "CG", fspace_degree) + # Get true solution spatial_coord = SpatialCoordinate(mesh) true_sol_expr = get_true_sol_expr(spatial_coord) + # TODO: Fix true sol + true_sol = Function(V).interpolate(true_sol_expr) + n = FacetNormal(mesh) + f = dot(grad(true_sol), n) - # Build function spaces - cgfspace = FunctionSpace(mesh, "CG", fspace_degree) - cgvfspace = VectorFunctionSpace(mesh, "CG", fspace_degree) - - # TODO: Fix RHS - # TODO: Fix true solution expr - true_sol = Function(cgfspace).interpolate(true_sol_expr) - # FIXME: Handle normal signs + # places has source of inner boundary and target of outer boundary + from ufl import derivative + import petsc4py.PETSc + petsc4py.PETSc.Sys.popErrorHandler() places = PotentialSourceAndTarget(mesh, - source_region_dim=2, - source_region_id=source_bdy_id, - target_region_dim=2, - target_region_id=target_bdy_id) - - Dtrue_sol = DoubleLayerPotential(dot(true_sol, n), actx=actx) - Strue_sol = SingleLayerPotential(dot(true_sol, n), actx=actx) - # get rhs form - v = TestFunction(cgfspace) - rhs_form = inner(dot(grad(true_sol), n), - v) * ds(source_bdy_id, metadata={'quadrature_degree': 2 * fspace_degree}) \ - - 1j * inner(Strue_sol, v) * ds(target_bdy_id) \ - - inner(dot(grad(Dtrue_sol), n), v) * ds(target_bdy_id) + source_region_dim=1, + source_region_id=scatterer_bdy, + target_region_dim=1, + target_region_id=truncated_bdy) + # TODO: Fix RHS + from sumpy.kernel import HelmholtzKernel + Sf = SingleLayerPotential(dot(grad(true_sol), n), + HelmholtzKernel(dim=2), + places, + actx=actx, + function_space=V, + op_kwargs={'k': kappa, 'qbx_forced_limit': None}) + v = TestFunction(V) + rhs = inner(f, v) * ds(scatterer_bdy) + \ + 1j * kappa * inner(Sf, v) * ds(truncated_bdy) - \ + inner(dot(grad(Sf), n), v) * ds(truncated_bdy) + + assemble(rhs) # TODO: Continue fixing test @@ -269,3 +277,9 @@ def get_true_sol_expr(spatial_coord): assert(eoc_recorder.order_estimate() >= fspace_degree or eoc_recorder.max_error() < 2e-14) + + +fspace_degree = 1 +kappa = 1.0 +from pyopencl import create_some_context +test_sommerfeld_helmholtz(create_some_context, fspace_degree, kappa) From 07f5ce25d933c44a5d67cac6ff20f370ffc4d23c Mon Sep 17 00:00:00 2001 From: benSepanski Date: Thu, 18 Feb 2021 19:51:25 -0600 Subject: [PATCH 08/20] Chasing down more bugs --- firedrake/potential_evaluation/__init__.py | 43 +++++++++++++++------ firedrake/potential_evaluation/pytential.py | 9 +++-- 2 files changed, 36 insertions(+), 16 deletions(-) diff --git a/firedrake/potential_evaluation/__init__.py b/firedrake/potential_evaluation/__init__.py index 2ab4ff0241..90f9a9579b 100644 --- a/firedrake/potential_evaluation/__init__.py +++ b/firedrake/potential_evaluation/__init__.py @@ -49,39 +49,58 @@ class Potential(AbstractExternalOperator): """ _external_operator_type = 'GLOBAL' - def __init__(self, density, connection, potential_operator, **kwargs): + def __init__(self, density, **kwargs): """ :arg density: A :mod:`firedrake` :class:`firedrake.function.Function` or UFL expression which represents the density :math:`u` - :kwarg connection: A :class:`PotentialEvaluationLibraryConnection` - :kwarg potential_operator: The external potential evaluation library - bound operator. + + :kwarg operator_data: Should have keys + + :kwarg connection: A :class:`PotentialEvaluationLibraryConnection` + :kwarg potential_operator: The external potential evaluation library + bound operator. + :kwarg function_space: the function space """ # super AbstractExternalOperator.__init__(self, density, **kwargs) + if 'operator_data' not in kwargs: + raise ValueError("Missing kwarg 'operator_data'") + operator_data = kwargs['operator_data'] + # Get connection & bound op and validate + if 'connection' not in operator_data: + raise ValueError("Missing operator_data 'connection'") + connection = operator_data['connection'] if not isinstance(connection, PotentialEvaluationLibraryConnection): raise TypeError("connection must be of type " "PotentialEvaluationLibraryConnection, not %s." % type(connection)) - self.connection = connection - self.potential_operator = potential_operator + if 'potential_operator' not in operator_data: + raise ValueError("Missing operator_data 'potential_operator'") # Get function space and validate aginst bound op - function_space = self.ufl_function_space() - assert isinstance(function_space, WithGeometry) # sanity check + if 'function_space' not in kwargs: + raise ValueError("Missing kwarg 'function_space'") + + function_space = kwargs['function_space'] + if not isinstance(function_space, WithGeometry): + raise TypeError("function_space must be of type WithGeometry, " + f"not {type(function_space)}") # Make sure density is a member of our function space, if it is if isinstance(density, Function): if density.function_space() is not function_space: raise ValueError("density.function_space() must be the " "same as function_space") + # save attributes + self.density = density + self.connection = connection + self.potential_operator = operator_data['potential_operator'] @cached_property def _evaluator(self): - return PotentialEvaluator(self, - self.density, + return PotentialEvaluator(self.density, self.connection, self.potential_operator) @@ -92,10 +111,10 @@ def _compute_derivatives(self, continuity_tolerance=None): raise NotImplementedError def _evaluate_action(self, *args): - return self._evaluator._evaluate_action() + return self._evaluator._evaluate_action(*args) def _evaluate_adjoint_action(self, *args): - return self._evaluator._evaluate_action() + return self._evaluator._evaluate_action(*args) def evaluate_adj_component_control(self, x, idx): derivatives = tuple(dj + int(idx == j) for j, dj in enumerate(self.derivatives)) diff --git a/firedrake/potential_evaluation/pytential.py b/firedrake/potential_evaluation/pytential.py index 06affba5e7..90eac8b82e 100644 --- a/firedrake/potential_evaluation/pytential.py +++ b/firedrake/potential_evaluation/pytential.py @@ -220,7 +220,8 @@ def bound_op_with_kwargs(density_arg): return pyt_op(**{density_name: density_arg}) # Now build and return Potential object - return Potential(density, - connection=meshmode_connection, - potential_operator=bound_op_with_kwargs, - **kwargs) + operator_data = {'connection': meshmode_connection, + 'potential_operator': bound_op_with_kwargs, + } + kwargs['operator_data'] = operator_data + return Potential(density, **kwargs) From 10f6c8ddedc30a2af1a3774d3614257dc70b98d5 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Fri, 19 Feb 2021 12:04:16 -0600 Subject: [PATCH 09/20] Trudging through UFL using test --- firedrake/potential_evaluation/__init__.py | 31 +++++++---- firedrake/potential_evaluation/pytential.py | 2 +- .../test_layer_potentials.py | 54 ++++++++++++++----- 3 files changed, 64 insertions(+), 23 deletions(-) diff --git a/firedrake/potential_evaluation/__init__.py b/firedrake/potential_evaluation/__init__.py index 90f9a9579b..4851d4a059 100644 --- a/firedrake/potential_evaluation/__init__.py +++ b/firedrake/potential_evaluation/__init__.py @@ -100,9 +100,11 @@ def __init__(self, density, **kwargs): @cached_property def _evaluator(self): - return PotentialEvaluator(self.density, + return PotentialEvaluator(self, + self.density, self.connection, - self.potential_operator) + self.potential_operator, + self.function_space()) def _evaluate(self): raise NotImplementedError @@ -130,31 +132,35 @@ class PotentialEvaluator: """ Evaluates a potential """ - def __init__(self, density, connection, potential_operator): + def __init__(self, potential, density, connection, potential_operator, function_space): """ :arg density: The UFL/firedrake density function :arg connection: A :class:`PotentialEvaluationLibraryConnection` :arg potential_operator: The external potential evaluation library bound operator. + :arg function_space: the :mod:`firedrake` function space of this + operator """ + self.potential = potential self.density = density self.connection = connection self.potential_operator = potential_operator + self.function_space = function_space - def _eval_potential_operator(self, action_coefficients, out=None): + def _eval_potential_operator(self, density, out=None): """ Evaluate the potential on the action coefficients and return. If *out* is not *None*, stores the result in *out* - :arg action_coefficients: A :class:`~firedrake.function.Function` - to evaluate the potential at + :arg density: the density :arg out: If not *None*, then a :class:`~firedrake.function.Function` in which to store the evaluated potential :return: *out* if it is not *None*, otherwise a new :class:`firedrake.function.Function` storing the evaluated potential """ - density = self.connection.from_firedrake(action_coefficients) + density_discrete = interpolate(density, self.function_space) + density = self.connection.from_firedrake(density_discrete) potential = self.potential_operator(density) return self.connection.to_firedrake(potential, out=out) @@ -162,7 +168,7 @@ def _evaluate(self): """ Evaluate P(self.density) into self """ - return self._eval_potential_operator(self.density, out=self) + return self._eval_potential_operator(self.density, out=self.potential) def _compute_derivatives(self): """ @@ -172,14 +178,19 @@ def _compute_derivatives(self): # FIXME : Assumes derivatives are Jacobian return self._eval_potential_operator - def _evaluate_action(self, action_coefficients): + def _evaluate_action(self, args): """ Evaluate derivatives of layer potential at action coefficients. i.e. Derivative(P, self.derivatives, self.density) evaluated at the action coefficients and store into self """ + if len(args) == 0: + return self._evaluate() + + # FIXME: Assumes just taking Jacobian + assert len(args) == 1 operator = self._compute_derivatives() - return operator(action_coefficients, out=self) + return operator(*args, out=self.potential) class PotentialEvaluationLibraryConnection: diff --git a/firedrake/potential_evaluation/pytential.py b/firedrake/potential_evaluation/pytential.py index 90eac8b82e..9711b904e4 100644 --- a/firedrake/potential_evaluation/pytential.py +++ b/firedrake/potential_evaluation/pytential.py @@ -138,7 +138,7 @@ def to_firedrake(self, evaluated_potential, out=None): # Otherwise, we have to project back to our function space pot = \ self.target_to_meshmode_connection.from_meshmode(evaluated_potential) - pot = project(evaluated_potential, self.function_space) + pot = project(pot, self.function_space) if out is not None: out.dat.data[:] = pot.dat.data[:] pot = out diff --git a/tests/pointwiseoperator/test_layer_potentials.py b/tests/pointwiseoperator/test_layer_potentials.py index f5e29dd271..195cfad803 100644 --- a/tests/pointwiseoperator/test_layer_potentials.py +++ b/tests/pointwiseoperator/test_layer_potentials.py @@ -106,19 +106,35 @@ def test_sommerfeld_helmholtz(ctx_factory, fspace_degree, kappa): # We'll use this to test convergence eoc_recorder = EOCRecorder() - def get_true_sol_expr(spatial_coord): + def get_true_sol(fspace, kappa, cl_ctx, queue): """ - Get the ufl expression for the true solution + Get the ufl expression for the true solution (3D) + or a function with the evaluated solution (2D) """ - mesh_dim = len(spatial_coord) + mesh_dim = fspace.mesh().geometric_dimension() if mesh_dim == 3: + spatial_coord = SpatialCoordinate(fspace.mesh()) x, y, z = spatial_coord # pylint: disable=C0103 norm = sqrt(x**2 + y**2 + z**2) return Constant(1j / (4*pi)) / norm * exp(1j * kappa * norm) if mesh_dim == 2: - x, y = spatial_coord # pylint: disable=C0103 - return Constant(1j / 4) * hankel_function(kappa * sqrt(x**2 + y**2), n=80) + # Evaluate true-sol using sumpy + from sumpy.p2p import P2P + from sumpy.kernel import HelmholtzKernel + # https://github.com/inducer/sumpy/blob/900745184d2618bc27a64c847f247e01c2b90b02/examples/curve-pot.py#L87-L88 + p2p = P2P(cl_ctx, [HelmholtzKernel(dim=2)], exclude_self=False, + value_dtypes=np.complex128) + # source is just (0, 0) + sources = np.array([[0.0], [0.0]]) + strengths = np.array([[1.0], [1.0]]) + # targets are everywhere + targets = np.array([Function(fspace).interpolate(x_i).dat.data + for x_i in SpatialCoordinate(fspace.mesh())]) + evt, (true_sol_arr,) = p2p(queue, targets, sources, strengths, k=kappa) + true_sol = Function(fspace) + true_sol.dat.data[:] = true_sol_arr[:] + return true_sol raise ValueError("Only meshes of dimension 2, 3 supported") # Create mesh and build hierarchy @@ -127,21 +143,34 @@ def get_true_sol_expr(spatial_coord): element_size=0.5, levels=3, order=2, - project_refinements_to_cad=False) + project_refinements_to_cad=False, + cache=False) scatterer_bdy = 5 # (inner boundary) the circle truncated_bdy = (1, 2, 3, 4) # (outer boundary) the square # Solve for each mesh in hierarchy for h, mesh in zip([0.5 * 2**i for i in range(len(mesh_hierarchy))], mesh_hierarchy): # Build function spaces V = FunctionSpace(mesh, "CG", fspace_degree) + Vdg = FunctionSpace(mesh, "DG", fspace_degree) # Get true solution - spatial_coord = SpatialCoordinate(mesh) - true_sol_expr = get_true_sol_expr(spatial_coord) - # TODO: Fix true sol - true_sol = Function(V).interpolate(true_sol_expr) + true_sol = get_true_sol(Vdg, kappa, cl_ctx, queue) + + # {{{ Get Neumann Data + n = FacetNormal(mesh) - f = dot(grad(true_sol), n) + f_expr = dot(grad(true_sol), n) + v = TestFunction(Vdg) + u = TrialFunction(Vdg) + a = inner(u, v) * ds(scatterer_bdy) + inner(u, v) * dx + L = inner(f_expr, v) * ds(scatterer_bdy) + inner(Constant(0.0), v) * dx + from firedrake import DirichletBC + bc = DirichletBC(Vdg, Constant(0.0), truncated_bdy) + f = Function(Vdg) + solve(a == L, f, bcs=[bc]) + f = project(f, V) + + # }}} # FIXME: Handle normal signs # places has source of inner boundary and target of outer boundary @@ -155,7 +184,7 @@ def get_true_sol_expr(spatial_coord): target_region_id=truncated_bdy) # TODO: Fix RHS from sumpy.kernel import HelmholtzKernel - Sf = SingleLayerPotential(dot(grad(true_sol), n), + Sf = SingleLayerPotential(f, HelmholtzKernel(dim=2), places, actx=actx, @@ -167,6 +196,7 @@ def get_true_sol_expr(spatial_coord): inner(dot(grad(Sf), n), v) * ds(truncated_bdy) assemble(rhs) + 1/0 # TODO: Continue fixing test From 0de18393b4c167f9a304ed2cd40fd777d2e73fc2 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Sat, 20 Feb 2021 15:04:29 -0600 Subject: [PATCH 10/20] pass local_operands to ufl reconstruct --- firedrake/pointwise_operators.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/firedrake/pointwise_operators.py b/firedrake/pointwise_operators.py index fcda901577..99019ce728 100644 --- a/firedrake/pointwise_operators.py +++ b/firedrake/pointwise_operators.py @@ -196,7 +196,7 @@ def evaluate_adj_component_state(self, x, idx): def _split(self): return tuple(Function(V, val) for (V, val) in zip(self.function_space(), self.topological.split())) - def _ufl_expr_reconstruct_(self, *operands, function_space=None, derivatives=None, name=None, operator_data=None, val=None, coefficient=None, arguments=None, add_kwargs={}): + def _ufl_expr_reconstruct_(self, *operands, function_space=None, derivatives=None, name=None, operator_data=None, val=None, coefficient=None, arguments=None, local_operands=None, add_kwargs={}): "Return a new object of the same type with new operands." deriv_multiindex = derivatives or self.derivatives @@ -211,6 +211,7 @@ def _ufl_expr_reconstruct_(self, *operands, function_space=None, derivatives=Non name=name, coefficient=coefficient, arguments=arguments, + local_operands=local_operands, operator_data=operator_data, add_kwargs=add_kwargs) else: @@ -221,6 +222,7 @@ def _ufl_expr_reconstruct_(self, *operands, function_space=None, derivatives=Non name=name or self.name(), coefficient=corresponding_coefficient, arguments=arguments or (self.arguments()+self.action_coefficients()), + local_operands=self.local_operands, operator_data=operator_data or self.operator_data, **add_kwargs) From 8a40b9fcb692c6f6b78d18dab5fc7a90b61f9739 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Sat, 20 Feb 2021 15:04:49 -0600 Subject: [PATCH 11/20] Pass kwargs to Potential --- firedrake/potential_evaluation/__init__.py | 8 +++---- firedrake/potential_evaluation/potentials.py | 24 +++++++------------- firedrake/potential_evaluation/pytential.py | 2 ++ 3 files changed, 13 insertions(+), 21 deletions(-) diff --git a/firedrake/potential_evaluation/__init__.py b/firedrake/potential_evaluation/__init__.py index 4851d4a059..d7194e359f 100644 --- a/firedrake/potential_evaluation/__init__.py +++ b/firedrake/potential_evaluation/__init__.py @@ -47,8 +47,6 @@ class Potential(AbstractExternalOperator): 0 & x \notin \Sigma \end{cases} """ - _external_operator_type = 'GLOBAL' - def __init__(self, density, **kwargs): """ :arg density: A :mod:`firedrake` :class:`firedrake.function.Function` @@ -107,10 +105,10 @@ def _evaluator(self): self.function_space()) def _evaluate(self): - raise NotImplementedError + return self._evaluator._evaluate() - def _compute_derivatives(self, continuity_tolerance=None): - raise NotImplementedError + def _compute_derivatives(self): + return self._evaluator._compute_derivatives() def _evaluate_action(self, *args): return self._evaluator._evaluate_action(*args) diff --git a/firedrake/potential_evaluation/potentials.py b/firedrake/potential_evaluation/potentials.py index fabe683a02..bdec2203ec 100644 --- a/firedrake/potential_evaluation/potentials.py +++ b/firedrake/potential_evaluation/potentials.py @@ -12,6 +12,10 @@ def SingleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): :arg potential_src_and_tgt: A :class:`firedrake.potential_evaluation.PotentialSourceAndTarget` + Any kwargs not listed below are passed on to the + :class:`firedrake.potential_evaluation.Potential` class + (along with some of the kwargs which are listed below). + :kwarg function_space: If density is a firedrake function, this is optional. Otherwise it must be supplied. If supplied and density is a function, they must @@ -87,7 +91,6 @@ def _layer_potential(layer_potential_sym, actx = kwargs['actx'] # make sure we got a function_space during validation assert 'function_space' in kwargs - function_space = kwargs['function_space'] # extract Pytential operation kwargs pyt_op_kwargs = _extract_pytential_operation_kwargs(**kwargs) # now return the pytential operation as an external operator @@ -151,13 +154,6 @@ def _validate_layer_potential_args_and_set_defaults(density, raise ValueError("source of a layer potential must have co-dimension 1," f" not {mesh_gdim - src_tdim}.") - # Make sure all kwargs are recognized - allowed_kwargs = ("cl_ctx", "queue", "actx", "op_kwargs", "qbx_kwargs", - "meshmode_connection_kwargs", "function_space") - for key in kwargs: - if key not in allowed_kwargs: - raise ValueError(f"Unrecognized kwarg {key}") - # Now handle pyopencl computing contexts and build # a PyOpenCLArrayContext from meshmode.array_context import PyOpenCLArrayContext @@ -226,14 +222,10 @@ def _extract_pytential_operation_kwargs(**kwargs): """ Extract kwargs to be passed to :func:`PytentialOperation` """ - pyt_op_kwargs = {} - pyt_op_possible_keywords = ("warn_if_cg", - "meshmode_connection_kwargs", - "qbx_kwargs", - "function_space") - for key in pyt_op_possible_keywords: - if key in kwargs: - pyt_op_kwargs[key] = kwargs[key] + keywords_to_remove = ("cl_ctx", "queue", "actx", "op_kwargs") + pyt_op_kwargs = dict(kwargs) + for key in keywords_to_remove: + pyt_op_kwargs.pop(key, None) return pyt_op_kwargs diff --git a/firedrake/potential_evaluation/pytential.py b/firedrake/potential_evaluation/pytential.py index 9711b904e4..ad96eace6e 100644 --- a/firedrake/potential_evaluation/pytential.py +++ b/firedrake/potential_evaluation/pytential.py @@ -223,5 +223,7 @@ def bound_op_with_kwargs(density_arg): operator_data = {'connection': meshmode_connection, 'potential_operator': bound_op_with_kwargs, } + if 'operator_data' in kwargs: + raise ValueError("User may not pass 'operator_data' as kwarg") kwargs['operator_data'] = operator_data return Potential(density, **kwargs) From 7c3d31d3908fdae354f2a16d656812a6cc06517e Mon Sep 17 00:00:00 2001 From: benSepanski Date: Fri, 26 Feb 2021 17:01:51 -0600 Subject: [PATCH 12/20] Fixed copy-paste error --- firedrake/potential_evaluation/__init__.py | 46 ++++++++++++++++++---- 1 file changed, 39 insertions(+), 7 deletions(-) diff --git a/firedrake/potential_evaluation/__init__.py b/firedrake/potential_evaluation/__init__.py index d7194e359f..04ddf4dadf 100644 --- a/firedrake/potential_evaluation/__init__.py +++ b/firedrake/potential_evaluation/__init__.py @@ -88,7 +88,7 @@ def __init__(self, density, **kwargs): # Make sure density is a member of our function space, if it is if isinstance(density, Function): - if density.function_space() is not function_space: + if density.function_space() != function_space: raise ValueError("density.function_space() must be the " "same as function_space") # save attributes @@ -145,6 +145,17 @@ def __init__(self, potential, density, connection, potential_operator, function_ self.potential_operator = potential_operator self.function_space = function_space + import matplotlib.pyplot as plt + from meshmode.mesh.visualization import draw_2d_mesh, draw_curve + src_discr = connection.get_source_discretization() + tgt_discr = connection.get_target_discretization() + draw_curve(src_discr.mesh) + plt.title("SOURCE") + plt.show() + draw_2d_mesh(tgt_discr.mesh, set_bounding_box=True) + plt.title("TARGET") + plt.show() + def _eval_potential_operator(self, density, out=None): """ Evaluate the potential on the action coefficients and return. @@ -317,8 +328,8 @@ def __init__(self, mesh, integer or tuple) or *None*. *None* indicates either the entire mesh, or the entire exterior boundary (as determined by the value of region). - The string "everywhere" is also equivalent to a - value of *None* + The string "on_boundary"/"everywhere" is also equivalent to a + value of *None* (for the appropriate subdomain dimension) NOTE: For implementation reasons, if the target is a boundary, instead of just evaluating the @@ -363,6 +374,21 @@ def __init__(self, mesh, [source_region_dim, target_region_dim]): if id_ is None or id_ == "everywhere": continue + if id_ == "everywhere": + if dim != mesh.topological_dimension(): + raise ValueError("\"everywhere\" only valid for subdomain" + " of same dimension as mesh. Try " + "\"on_boundary\" instead") + continue + if id_ == "on_boundary": + if dim != mesh.topological_dimension()-1: + raise ValueError("\"on_boundary\" only valid for subdomain" + " of dimension one less than the mesh. Try " + "\"everywhere\" instead") + continue + if isinstance(id_, str): + raise TypeError("Only valid strings are \"on_boundary\" and " + f"\"everywhere\", not '{id_}'.") # standardize id_ to a tuple if isinstance(id_, int): id_ = tuple([id_]) @@ -394,16 +420,22 @@ def __init__(self, mesh, # handle None case if source_region_id is None: - source_region_id = "everywhere" + if source_region_dim == mesh.topological_dimension(): + source_region_id = "everywhere" + else: + source_region_id = "on_boundary" if target_region_id is None: - target_region_id = "everywhere" + if target_region_dim == mesh.topological_dimension(): + target_region_id = "everywhere" + else: + target_region_id = "on_boundary" # store mesh, region dims, and region ids self.mesh = mesh self._source_region_dim = source_region_dim - self._target_region_dim = source_region_dim + self._target_region_dim = target_region_dim self._source_region_id = source_region_id - self._target_region_id = source_region_id + self._target_region_id = target_region_id def get_source_dimension(self): """ From 45c64e1d72aed6ae55f238d5f49d5854ac777a05 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Tue, 2 Mar 2021 10:37:56 -0600 Subject: [PATCH 13/20] remove evaluator --- firedrake/potential_evaluation/__init__.py | 114 +++++++++------------ 1 file changed, 46 insertions(+), 68 deletions(-) diff --git a/firedrake/potential_evaluation/__init__.py b/firedrake/potential_evaluation/__init__.py index 04ddf4dadf..3943dffaa8 100644 --- a/firedrake/potential_evaluation/__init__.py +++ b/firedrake/potential_evaluation/__init__.py @@ -96,66 +96,6 @@ def __init__(self, density, **kwargs): self.connection = connection self.potential_operator = operator_data['potential_operator'] - @cached_property - def _evaluator(self): - return PotentialEvaluator(self, - self.density, - self.connection, - self.potential_operator, - self.function_space()) - - def _evaluate(self): - return self._evaluator._evaluate() - - def _compute_derivatives(self): - return self._evaluator._compute_derivatives() - - def _evaluate_action(self, *args): - return self._evaluator._evaluate_action(*args) - - def _evaluate_adjoint_action(self, *args): - return self._evaluator._evaluate_action(*args) - - def evaluate_adj_component_control(self, x, idx): - derivatives = tuple(dj + int(idx == j) for j, dj in enumerate(self.derivatives)) - dN = self._ufl_expr_reconstruct_(*self.ufl_operands, derivatives=derivatives) - return dN._evaluate_adjoint_action((x.function,)).vector() - - def evaluate_adj_component(self, x, idx): - print(x, type(x)) - raise NotImplementedError - - -class PotentialEvaluator: - """ - Evaluates a potential - """ - def __init__(self, potential, density, connection, potential_operator, function_space): - """ - :arg density: The UFL/firedrake density function - :arg connection: A :class:`PotentialEvaluationLibraryConnection` - :arg potential_operator: The external potential evaluation library - bound operator. - :arg function_space: the :mod:`firedrake` function space of this - operator - """ - self.potential = potential - self.density = density - self.connection = connection - self.potential_operator = potential_operator - self.function_space = function_space - - import matplotlib.pyplot as plt - from meshmode.mesh.visualization import draw_2d_mesh, draw_curve - src_discr = connection.get_source_discretization() - tgt_discr = connection.get_target_discretization() - draw_curve(src_discr.mesh) - plt.title("SOURCE") - plt.show() - draw_2d_mesh(tgt_discr.mesh, set_bounding_box=True) - plt.title("TARGET") - plt.show() - def _eval_potential_operator(self, density, out=None): """ Evaluate the potential on the action coefficients and return. @@ -168,7 +108,7 @@ def _eval_potential_operator(self, density, out=None): :class:`firedrake.function.Function` storing the evaluated potential """ - density_discrete = interpolate(density, self.function_space) + density_discrete = interpolate(density, self.function_space()) density = self.connection.from_firedrake(density_discrete) potential = self.potential_operator(density) return self.connection.to_firedrake(potential, out=out) @@ -177,29 +117,66 @@ def _evaluate(self): """ Evaluate P(self.density) into self """ - return self._eval_potential_operator(self.density, out=self.potential) + return self._eval_potential_operator(self.density, out=self) def _compute_derivatives(self): """ Return a function Derivative(P, self.derivatives, self.density) """ - # FIXME : Assumes derivatives are Jacobian + # FIXME : Assumes only one ufl operand, is that okay? + assert len(self.ufl_operands) == 1 + assert self.ufl_operands[0] is self.density + assert len(self.derivatives) == 1 + # Derivative(P, (0,), self.density)(density_star) + # = P(self.density) + if self.derivatives == (0,): + return lambda density_star, out=None: \ + self._eval_potential_operator(self.density, out=out) + # P is linear, so the nth Gateaux derivative + # Derivative(P, (n,), u)(u*) = P(u*) + # + # d/dx (P \circ u)(v) + # = lim_{t->0} (P(u(v+tx)) - P(u)) / t + # \approx lim_{t->0} (P(u(v) + t du/dx - P(u(v))) / t + # = P(du/dx) + # + # So d^n/dx^n( P \circ u ) = P(d^nu/dx^n) return self._eval_potential_operator - def _evaluate_action(self, args): + def _evaluate_action(self, *args): """ Evaluate derivatives of layer potential at action coefficients. i.e. Derivative(P, self.derivatives, self.density) evaluated at the action coefficients and store into self """ - if len(args) == 0: + assert len(args) == 1 # sanity check + # Either () or (density_star,) + density_star = args[0] + assert len(density_star) in [0, 1] # sanity check + # Evaluate Derivative(P, self.derivatives, self.density) at density_star + if len(density_star) == 0: return self._evaluate() # FIXME: Assumes just taking Jacobian - assert len(args) == 1 + density_star, = density_star operator = self._compute_derivatives() - return operator(*args, out=self.potential) + return operator(density_star, out=self) + + def _evaluate_adjoint_action(self, *args): + """ + Operator is self-adjoint, so just evaluate action + """ + return self._evaluate_action(*args) + + def evaluate_adj_component_control(self, x, idx): + derivatives = tuple(dj + int(idx == j) for j, dj in enumerate(self.derivatives)) + dN = self._ufl_expr_reconstruct_(*self.ufl_operands, derivatives=derivatives) + return dN._evaluate_adjoint_action((x.function,)).vector() + + def evaluate_adj_component(self, x, idx): + print(x, type(x)) + raise NotImplementedError class PotentialEvaluationLibraryConnection: @@ -260,8 +237,9 @@ def __init__(self, function_space, potential_source_and_target, if shape is None or len(shape) == 0: dg_function_space = FunctionSpace(mesh, "DG", degree) elif len(shape) == 1: + dim = shape[0] dg_function_space = VectorFunctionSpace(mesh, "DG", degree, - dim=shape) + dim=dim) else: dg_function_space = TensorFunctionSpace(mesh, "DG", degree, shape=shape) From 8b7ce906f57bc13c0a1c833f5b6a3546a0559170 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Tue, 2 Mar 2021 10:38:26 -0600 Subject: [PATCH 14/20] Make sure to use on_boundary not everywhere --- firedrake/potential_evaluation/pytential.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/firedrake/potential_evaluation/pytential.py b/firedrake/potential_evaluation/pytential.py index ad96eace6e..8e2d1d7323 100644 --- a/firedrake/potential_evaluation/pytential.py +++ b/firedrake/potential_evaluation/pytential.py @@ -69,11 +69,17 @@ def __init__(self, function_space, potential_source_and_target, actx, # boundary restrict_conn = None if src_bdy is not None: - # convert "everywhere" to meshmode BTAG_ALL + # convert "on_boundary" to meshmode BTAG_ALL meshmode_src_bdy = src_bdy - if meshmode_src_bdy == "everywhere": + if meshmode_src_bdy == "on_boundary": from meshmode.mesh import BTAG_ALL meshmode_src_bdy = BTAG_ALL + else: + # sanity check + assert isinstance(src_bdy, int) or isinstance(src_bdy, tuple) + if isinstance(src_bdy, tuple): + for id_ in src_bdy: + assert isinstance(id_, int) # get group factory order = self.dg_function_space.ufl_element().degree() from meshmode.discretization.poly_element import \ From 1d8a8a34cc40d801db3df84b70efd7fd86594016 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Tue, 2 Mar 2021 11:07:48 -0600 Subject: [PATCH 15/20] Broadcast pytential operation over arrays --- firedrake/potential_evaluation/pytential.py | 25 ++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/firedrake/potential_evaluation/pytential.py b/firedrake/potential_evaluation/pytential.py index 8e2d1d7323..2cfa03ef8b 100644 --- a/firedrake/potential_evaluation/pytential.py +++ b/firedrake/potential_evaluation/pytential.py @@ -1,3 +1,5 @@ +import numpy as np + from firedrake.potential_evaluation import \ Potential, PotentialEvaluationLibraryConnection from firedrake import project @@ -221,9 +223,30 @@ def PytentialOperation(actx, from pytential import bind pyt_op = bind((qbx, tgt_discr), unbound_op) + # Import for type-checks + from meshmode.dof_array import DOFArray + # build bound operator that just takes density as argument def bound_op_with_kwargs(density_arg): - return pyt_op(**{density_name: density_arg}) + """ + Apply operator to density_arg, or each DOFArray inside of + the density_arg + """ + # If given a DOFArray, compute on it + if isinstance(density_arg, DOFArray): + return pyt_op(**{density_name: density_arg}) + # TODO: Do we need these? may make things slower + # type checks + assert isinstance(density_arg, np.ndarray) + assert density_arg.dtype == np.object + # Otherwise, compute for each multi-index + result = np.empty(shape=density_arg.shape, + dtype=np.object) + for multi_index in np.ndindex(density_arg.shape): + # TODO: Probably should remove this + assert isinstance(density_arg[multi_index], DOFArray) + result[multi_index] = pyt_op(**{density_name: density_arg[multi_index]}) + return result # Now build and return Potential object operator_data = {'connection': meshmode_connection, From 921d9e341c19e60880de8df567cb37d41648cfc5 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Tue, 2 Mar 2021 11:07:59 -0600 Subject: [PATCH 16/20] Add module description --- firedrake/potential_evaluation/potentials.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/firedrake/potential_evaluation/potentials.py b/firedrake/potential_evaluation/potentials.py index bdec2203ec..c1dbb4e950 100644 --- a/firedrake/potential_evaluation/potentials.py +++ b/firedrake/potential_evaluation/potentials.py @@ -1,3 +1,8 @@ +""" +User-facing functions to create layer potentials. + +Produces a Potential either through volumential or pytential +""" from firedrake.functionspaceimpl import WithGeometry @@ -86,6 +91,7 @@ def _layer_potential(layer_potential_sym, unbound_op = layer_potential_sym(kernel, sym.var("density"), **op_kwargs) + # make sure we got an actx during validation assert 'actx' in kwargs actx = kwargs['actx'] From 5dc8f438e54f35206b3af5f3be6dd522e6384df2 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Tue, 2 Mar 2021 11:35:33 -0600 Subject: [PATCH 17/20] cache connection on mesh --- firedrake/potential_evaluation/pytential.py | 62 +++++++++++++++------ 1 file changed, 45 insertions(+), 17 deletions(-) diff --git a/firedrake/potential_evaluation/pytential.py b/firedrake/potential_evaluation/pytential.py index 2cfa03ef8b..84d4108a3e 100644 --- a/firedrake/potential_evaluation/pytential.py +++ b/firedrake/potential_evaluation/pytential.py @@ -3,6 +3,30 @@ from firedrake.potential_evaluation import \ Potential, PotentialEvaluationLibraryConnection from firedrake import project +from firedrake.functionspacedata import cached + + +@cached +def _build_connection_from_firedrake(mesh, key, + actx, + function_space, + restrict_to_boundary, + **meshmode_connection_kwargs): + """ + Build connection from firedrake to meshmode and cache on mesh based on key. + + :arg mesh: the mesh to cache on + :arg key: (function space ufl element, + restrict_to_boundary, + frozenset(meshmode_connection_kwargs.items())) + Other args are just passed to + :func:`meshmode.interop.firedrake.build_connection_from_firedrake` + """ + from meshmode.interop.firedrake import build_connection_from_firedrake + return build_connection_from_firedrake(actx, + function_space, + restrict_to_boundary=restrict_to_boundary, + **meshmode_connection_kwargs) class MeshmodeConnection(PotentialEvaluationLibraryConnection): @@ -62,11 +86,17 @@ def __init__(self, function_space, potential_source_and_target, actx, assert potential_source_and_target.get_source_dimension() == tdim - 1 src_bdy = potential_source_and_target.get_source_id() - from meshmode.interop.firedrake import build_connection_from_firedrake - src_conn = build_connection_from_firedrake(actx, - self.dg_function_space, - restrict_to_boundary=src_bdy, - **meshmode_connection_kwargs) + # Build source connection to (or retrieve cached connection from mesh) + mesh = self.dg_function_space.mesh() + key = (self.dg_function_space.ufl_element(), + src_bdy, + frozenset(meshmode_connection_kwargs.items())) + src_conn = _build_connection_from_firedrake(mesh, key, + actx, + self.dg_function_space, + restrict_to_boundary=src_bdy, + **meshmode_connection_kwargs) + # If source is a boundary, build a connection to restrict it to the # boundary restrict_conn = None @@ -103,18 +133,16 @@ def __init__(self, function_space, potential_source_and_target, actx, assert potential_source_and_target.get_target_dimension() == tdim - 1 tgt_bdy = potential_source_and_target.get_target_id() - # Can we re-use the source connection? - src_dim = potential_source_and_target.get_source_dimension() - tgt_dim = potential_source_and_target.get_target_dimension() - if tgt_bdy == src_bdy and src_dim == tgt_dim: - tgt_conn = src_conn - else: - # If not, build a new one - tgt_conn = build_connection_from_firedrake( - actx, - self.dg_function_space, - restrict_to_boundary=tgt_bdy, - **meshmode_connection_kwargs) + # Build target connection to (or retrieve cached connection from mesh) + mesh = self.dg_function_space.mesh() + key = (self.dg_function_space.ufl_element(), + tgt_bdy, + frozenset(meshmode_connection_kwargs.items())) + tgt_conn = _build_connection_from_firedrake(mesh, key, + actx, + self.dg_function_space, + restrict_to_boundary=tgt_bdy, + **meshmode_connection_kwargs) # store computing context self.actx = actx From b311247731706cddeb40a6bc045121931f9fc20d Mon Sep 17 00:00:00 2001 From: benSepanski Date: Tue, 2 Mar 2021 14:00:41 -0600 Subject: [PATCH 18/20] Fix caching to allow reuse even with vector versions of an fspace --- firedrake/potential_evaluation/pytential.py | 44 ++++++++++++--------- 1 file changed, 26 insertions(+), 18 deletions(-) diff --git a/firedrake/potential_evaluation/pytential.py b/firedrake/potential_evaluation/pytential.py index 84d4108a3e..0b06916a1f 100644 --- a/firedrake/potential_evaluation/pytential.py +++ b/firedrake/potential_evaluation/pytential.py @@ -7,16 +7,17 @@ @cached -def _build_connection_from_firedrake(mesh, key, - actx, - function_space, - restrict_to_boundary, - **meshmode_connection_kwargs): +def _cached_build_connection_from_firedrake(mesh, key, + actx, + function_space, + restrict_to_boundary, + **meshmode_connection_kwargs): """ Build connection from firedrake to meshmode and cache on mesh based on key. :arg mesh: the mesh to cache on - :arg key: (function space ufl element, + :arg key: the key this connection is looked up by + (function space ufl element, restrict_to_boundary, frozenset(meshmode_connection_kwargs.items())) Other args are just passed to @@ -29,6 +30,23 @@ def _build_connection_from_firedrake(mesh, key, **meshmode_connection_kwargs) +def _build_connection_from_firedrake(actx, function_space, restrict_to_boundary, + **meshmode_connection_kwargs): + """ + Constructs key and invokes _cached_build_connection_from_firedrake + """ + mesh = function_space.mesh() + key = (function_space.ufl_element().family(), + function_space.ufl_element().degree(), + restrict_to_boundary, + frozenset(meshmode_connection_kwargs.items())) + return _cached_build_connection_from_firedrake(mesh, key, + actx, + function_space, + restrict_to_boundary, + **meshmode_connection_kwargs) + + class MeshmodeConnection(PotentialEvaluationLibraryConnection): """ Build a connection into :mod:`meshmode` @@ -87,12 +105,7 @@ def __init__(self, function_space, potential_source_and_target, actx, src_bdy = potential_source_and_target.get_source_id() # Build source connection to (or retrieve cached connection from mesh) - mesh = self.dg_function_space.mesh() - key = (self.dg_function_space.ufl_element(), - src_bdy, - frozenset(meshmode_connection_kwargs.items())) - src_conn = _build_connection_from_firedrake(mesh, key, - actx, + src_conn = _build_connection_from_firedrake(actx, self.dg_function_space, restrict_to_boundary=src_bdy, **meshmode_connection_kwargs) @@ -134,12 +147,7 @@ def __init__(self, function_space, potential_source_and_target, actx, tgt_bdy = potential_source_and_target.get_target_id() # Build target connection to (or retrieve cached connection from mesh) - mesh = self.dg_function_space.mesh() - key = (self.dg_function_space.ufl_element(), - tgt_bdy, - frozenset(meshmode_connection_kwargs.items())) - tgt_conn = _build_connection_from_firedrake(mesh, key, - actx, + tgt_conn = _build_connection_from_firedrake(actx, self.dg_function_space, restrict_to_boundary=tgt_bdy, **meshmode_connection_kwargs) From 8361abefedf8ca5d5a56aab03b8290cb9fbd7370 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Fri, 16 Apr 2021 16:38:13 -0500 Subject: [PATCH 19/20] rearrange files for a clearer organization --- firedrake/potential_evaluation/__init__.py | 489 +++++++++---------- firedrake/potential_evaluation/potentials.py | 472 +++++++++--------- 2 files changed, 475 insertions(+), 486 deletions(-) diff --git a/firedrake/potential_evaluation/__init__.py b/firedrake/potential_evaluation/__init__.py index 3943dffaa8..a6da88f3ec 100644 --- a/firedrake/potential_evaluation/__init__.py +++ b/firedrake/potential_evaluation/__init__.py @@ -1,285 +1,248 @@ -import numpy as np -import firedrake.function +""" +User-facing functions to create layer potentials. +Produces a Potential either through volumential or pytential +""" from firedrake.functionspaceimpl import WithGeometry from firedrake.mesh import MeshGeometry -from firedrake.pointwise_operators import AbstractExternalOperator -from firedrake.utils import cached_property - -from firedrake import Function, FunctionSpace, interpolate, Interpolator, \ - SpatialCoordinate, VectorFunctionSpace, TensorFunctionSpace, project - -from ufl.algorithms.apply_derivatives import VariableRuleset -from ufl.algorithms import extract_coefficients -from ufl.constantvalue import as_ufl -from ufl.core.multiindex import indices -from ufl.tensors import as_tensor - -from pyop2.datatypes import ScalarType -from warnings import warn - - -from firedrake.potential_evaluation.potentials import \ - SingleLayerPotential, DoubleLayerPotential, VolumePotential __all__ = ("SingleLayerPotential", "DoubleLayerPotential", "VolumePotential", "PotentialSourceAndTarget") -class Potential(AbstractExternalOperator): - r""" - This is a function which represents a potential - computed on a firedrake function. - - For a firedrake function :math:`u:\Omega\to\mathbb C^n` - with :math:`\Omega \subseteq \mathbb R^m`, - the potential :math:`P(u):\Omega \to\mathbb C^n` is a - convolution of :math:`u` against some kernel function. - More concretely, given a source region :math:`\Gamma \subseteq \Omega` - a target region :math:`\Sigma \subseteq \Omega`, and - a kernel function :math:`K`, we define - - .. math:: - - P(u)(x) = \begin{cases} - \int_\Gamma K(x, y) u(y) \,dy & x \in \Sigma - 0 & x \notin \Sigma - \end{cases} +def SingleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): """ - def __init__(self, density, **kwargs): - """ - :arg density: A :mod:`firedrake` :class:`firedrake.function.Function` - or UFL expression which represents the density :math:`u` - - :kwarg operator_data: Should have keys - - :kwarg connection: A :class:`PotentialEvaluationLibraryConnection` - :kwarg potential_operator: The external potential evaluation library - bound operator. - :kwarg function_space: the function space - """ - # super - AbstractExternalOperator.__init__(self, density, **kwargs) - - if 'operator_data' not in kwargs: - raise ValueError("Missing kwarg 'operator_data'") - operator_data = kwargs['operator_data'] - - # Get connection & bound op and validate - if 'connection' not in operator_data: - raise ValueError("Missing operator_data 'connection'") - connection = operator_data['connection'] - if not isinstance(connection, PotentialEvaluationLibraryConnection): - raise TypeError("connection must be of type " - "PotentialEvaluationLibraryConnection, not %s." - % type(connection)) - if 'potential_operator' not in operator_data: - raise ValueError("Missing operator_data 'potential_operator'") - - # Get function space and validate aginst bound op - if 'function_space' not in kwargs: - raise ValueError("Missing kwarg 'function_space'") - - function_space = kwargs['function_space'] - if not isinstance(function_space, WithGeometry): - raise TypeError("function_space must be of type WithGeometry, " - f"not {type(function_space)}") - - # Make sure density is a member of our function space, if it is - if isinstance(density, Function): - if density.function_space() != function_space: - raise ValueError("density.function_space() must be the " - "same as function_space") - # save attributes - self.density = density - self.connection = connection - self.potential_operator = operator_data['potential_operator'] - - def _eval_potential_operator(self, density, out=None): - """ - Evaluate the potential on the action coefficients and return. - If *out* is not *None*, stores the result in *out* - - :arg density: the density - :arg out: If not *None*, then a :class:`~firedrake.function.Function` - in which to store the evaluated potential - :return: *out* if it is not *None*, otherwise a new - :class:`firedrake.function.Function` storing the evaluated - potential - """ - density_discrete = interpolate(density, self.function_space()) - density = self.connection.from_firedrake(density_discrete) - potential = self.potential_operator(density) - return self.connection.to_firedrake(potential, out=out) - - def _evaluate(self): - """ - Evaluate P(self.density) into self - """ - return self._eval_potential_operator(self.density, out=self) - - def _compute_derivatives(self): - """ - Return a function - Derivative(P, self.derivatives, self.density) - """ - # FIXME : Assumes only one ufl operand, is that okay? - assert len(self.ufl_operands) == 1 - assert self.ufl_operands[0] is self.density - assert len(self.derivatives) == 1 - # Derivative(P, (0,), self.density)(density_star) - # = P(self.density) - if self.derivatives == (0,): - return lambda density_star, out=None: \ - self._eval_potential_operator(self.density, out=out) - # P is linear, so the nth Gateaux derivative - # Derivative(P, (n,), u)(u*) = P(u*) - # - # d/dx (P \circ u)(v) - # = lim_{t->0} (P(u(v+tx)) - P(u)) / t - # \approx lim_{t->0} (P(u(v) + t du/dx - P(u(v))) / t - # = P(du/dx) - # - # So d^n/dx^n( P \circ u ) = P(d^nu/dx^n) - return self._eval_potential_operator - - def _evaluate_action(self, *args): - """ - Evaluate derivatives of layer potential at action coefficients. - i.e. Derivative(P, self.derivatives, self.density) evaluated at - the action coefficients and store into self - """ - assert len(args) == 1 # sanity check - # Either () or (density_star,) - density_star = args[0] - assert len(density_star) in [0, 1] # sanity check - # Evaluate Derivative(P, self.derivatives, self.density) at density_star - if len(density_star) == 0: - return self._evaluate() - - # FIXME: Assumes just taking Jacobian - density_star, = density_star - operator = self._compute_derivatives() - return operator(density_star, out=self) - - def _evaluate_adjoint_action(self, *args): - """ - Operator is self-adjoint, so just evaluate action - """ - return self._evaluate_action(*args) + Evaluate the single layer potential of the density + function convoluted against the kernel on the source + at each target point. + + :arg density: The :mod:`firedrake` density function + :arg kernel: the :class:`sumpy.kernel.Kernel` + :arg potential_src_and_tgt: A + :class:`firedrake.potential_evaluation.PotentialSourceAndTarget` + + Any kwargs not listed below are passed on to the + :class:`firedrake.potential_evaluation.Potential` class + (along with some of the kwargs which are listed below). + + :kwarg function_space: If density is a firedrake function, + this is optional. Otherwise it must be supplied. + If supplied and density is a function, they must + be the same object. + :kwarg cl_ctx: (Optional) a :class:`pyopencl.Context` used + to create a command queue. + At most one of *cl_ctx*, *queue*, *actx* should be included. + If none is included, a *cl_ctx* will be created. + :kwarg queue: (Optional) a :class:`pyopencl.CommandQueue` used + to create an array context. + At most one of *cl_ctx*, *queue*, *actx* should be included. + If none is included, a *cl_ctx* will be created. + :kwarg actx: (Optional) a :class:`meshmode.array_context.PyOpenCLArrayContext` + used for :mod:`pytential` and :mod:`meshmode` computations. + At most one of *cl_ctx*, *queue*, *actx* should be included. + If none is included, a *cl_ctx* will be created. + :kwarg op_kwargs: (Optional) kwargs passed to the construction + of the unbound pytential operator (e.g. {'k': 0.5} for + a :class:`sumpy.kernel.HelmholtzKernel` must be passed + to :func:`pytential.symbolic.primitives.S`), + DEFAULT: + * {'qbx_forced_limit': None} + :kwarg qbx_kwargs: (Optional) *None*, or a dict with arguments to pass to + :class:`pytential.qbx.QBXLayerPotentialSource`. + DEFAULTS: + - fine_order: 4 8 function space degree + - qbx_order: function space degree + - fmm_order: *None* + - fmm_level_to_order: A + :class:`sumpy.level_to_order.FMMLibExpansionOrderFinder` set + to provide tolerance of double machine epsilon + 2**-53 with an extra order of 1. + This allows for taking 1 derivative, but should be + increased for taking more derivatives. + - fmm_backend: "fmmlib" + :kwarg meshmode_connection_kwargs: (Optional) + *None*, or a dict with arguments to pass to + :func:`meshmode.interop.firedrake.build_connection_from_firedrake` + :kwarg warn_if_cg: (Optional) Pass *False* to suppress the warning if + a "CG" space is used. *True* by default + """ + # TODO: Make fmm_level_to_order depend on number of derivatives somehow? + from pytential.symbolic.primitives import S + return _layer_potential(S, density, kernel, potential_src_and_tgt, **kwargs) - def evaluate_adj_component_control(self, x, idx): - derivatives = tuple(dj + int(idx == j) for j, dj in enumerate(self.derivatives)) - dN = self._ufl_expr_reconstruct_(*self.ufl_operands, derivatives=derivatives) - return dN._evaluate_adjoint_action((x.function,)).vector() - def evaluate_adj_component(self, x, idx): - print(x, type(x)) - raise NotImplementedError +def DoubleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): + """ + Same as SingleLayerPotential, but evaluates the double layer potential + instead of the single + """ + from pytential.symbolic.primitives import D + return _layer_potential(D, density, kernel, potential_src_and_tgt, **kwargs) -class PotentialEvaluationLibraryConnection: +def _layer_potential(layer_potential_sym, + density, kernel, potential_src_and_tgt, **kwargs): """ - A connection to an external library for potential evaluation + Build a layer potential. For usage, see + :func:`SingleLayerPotential` or :func:`DoubleLayerPotential`. """ - def __init__(self, function_space, potential_source_and_target, - warn_if_cg=True): - """ - Initialize self to work on the function space - - :arg function_space: The :mod:`firedrake` function space - (of type :class:`~firedrake.functionspaceimpl.WithGeometry`) on - which to convert to/from. Must be a 'Lagrange' or - 'Discontinuous Lagrange' space. - :arg potential_source_and_target: A :class:`PotentialSourceAndTarget`. - mesh must match that of function_space. - :arg warn_if_cg: If *True*, warn if the space is "CG" - """ - # validate function space - if not isinstance(function_space, WithGeometry): - raise TypeError("function_space must be of type WithGeometry, not " - "%s" % type(function_space)) - - family = function_space.ufl_element().family() - acceptable_families = ['Discontinuous Lagrange', 'Lagrange'] - if family not in acceptable_families: - raise ValueError("function_space.ufl_element().family() must be " - "one of %s, not '%s'" % - (acceptable_families, family)) - if family == 'Lagrange' and warn_if_cg: - warn("Functions in continuous function space will be projected " - "to/from a 'Discontinuous Lagrange' space. Make sure " - "any operators evaluated are continuous. " - "Pass warn_if_cg=False to suppress this warning.") - - # validate potential_source_and_targets - if not isinstance(potential_source_and_target, - PotentialSourceAndTarget): - raise TypeError("potential_source_and_targets must be of type " - "PotentialSourceAndTarget, not '%s'." - % type(potential_source_and_target)) - - # make sure meshes match - if potential_source_and_target.mesh is not function_space.mesh(): - raise ValueError("function_space.mesh() and " - "potential_source_and_target.mesh must be the same" - " obejct.") - - # build DG space if necessary - family = function_space.ufl_element().family() - if family == 'Discontinuous Lagrange': - dg_function_space = function_space - elif family == 'Lagrange': - mesh = function_space.mesh() - degree = function_space.ufl_element().degree() - shape = function_space.shape - if shape is None or len(shape) == 0: - dg_function_space = FunctionSpace(mesh, "DG", degree) - elif len(shape) == 1: - dim = shape[0] - dg_function_space = VectorFunctionSpace(mesh, "DG", degree, - dim=dim) - else: - dg_function_space = TensorFunctionSpace(mesh, "DG", degree, - shape=shape) + kwargs = _validate_layer_potential_args_and_set_defaults( + density, kernel, potential_src_and_tgt, **kwargs) + # build our unbound operator + assert 'op_kwargs' in kwargs # make sure we got op kwargs + op_kwargs = kwargs['op_kwargs'] + from pytential import sym + unbound_op = layer_potential_sym(kernel, + sym.var("density"), + **op_kwargs) + + # make sure we got an actx during validation + assert 'actx' in kwargs + actx = kwargs['actx'] + # make sure we got a function_space during validation + assert 'function_space' in kwargs + # extract Pytential operation kwargs + pyt_op_kwargs = _extract_pytential_operation_kwargs(**kwargs) + # now return the pytential operation as an external operator + from firedrake.potential_evaluation.pytential import PytentialOperation + return PytentialOperation(actx, density, unbound_op, "density", + potential_src_and_tgt, + **pyt_op_kwargs) + + +def _validate_layer_potential_args_and_set_defaults(density, + kernel, + places, + **kwargs): + """ + Validate the arguments for single/double layer potential. + + Returns dictionary updated with default arguments + """ + # validate density function space + if hasattr(density, 'function_space'): + if not isinstance(density.function_space(), WithGeometry): + raise TypeError("density.function_space() must be of type " + f"WithGeometry, not {type(density.function_space())}") + function_space = kwargs.get('function_space', None) + if function_space is not None: + if function_space is not density.function_space(): + raise ValueError("density.function_space() and function_space" + " must be the same object") else: - acceptable_families = ['Discontinuous Lagrange', 'Lagrange'] - raise ValueError("function_space.ufl_element().family() must be " - "one of %s, not '%s'" % - (acceptable_families, family)) - - # store function space and dg function space - self.function_space = function_space - self.dg_function_space = dg_function_space - self.is_dg = function_space == dg_function_space - # store source and targets - self.source_and_target = potential_source_and_target - - def from_firedrake(self, density): - """ - Convert the density into a form acceptable by an bound operation - in an external library + kwargs['function_space'] = density.function_space() + else: + function_space = kwargs.get('function_space', None) + if function_space is None: + raise ValueError("density has no function_space method, so" + "function_space kwarg must be supplied") + if not isinstance(function_space, WithGeometry): + raise TypeError("function_space must be of type " + f"WithGeometry, not {type(function_space)}") + function_space = kwargs['function_space'] + + # validate kernel + from sumpy.kernel import Kernel + if not isinstance(kernel, Kernel): + raise TypeError("kernel must be a sumpy.kernel.Kernel, not of " + f"type '{type(kernel)}'.") + # validate src and tgts + from firedrake.potential_evaluation import PotentialSourceAndTarget + if not isinstance(places, PotentialSourceAndTarget): + raise TypeError("potential_src_and_tgt must be a sumpy.kernel.Kernel, " + f"not of type '{type(places)}'.") + # Make sure src is of right dimension + mesh_gdim = places.mesh.geometric_dimension() + mesh_tdim = places.mesh.topological_dimension() + src_tdim = places.get_source_dimension() + # sanity check + assert mesh_gdim - mesh_tdim in [0, 1] + assert mesh_gdim - src_tdim in [0, 1] + assert mesh_tdim in (mesh_gdim, src_tdim) + # now do the real user-input check + if mesh_gdim - src_tdim != 1: + raise ValueError("source of a layer potential must have co-dimension 1," + f" not {mesh_gdim - src_tdim}.") + + # Now handle pyopencl computing contexts and build + # a PyOpenCLArrayContext + from meshmode.array_context import PyOpenCLArrayContext + cl_ctx = None + queue = None + actx = None + if 'cl_ctx' in kwargs: + if 'actx' in kwargs or 'queue' in kwargs: + raise ValueError("At most one of 'actx', 'queue', 'cl_ctx' should " + "be supplied") + cl_ctx = kwargs['cl_ctx'] + from pyopencl import Context + if not isinstance(cl_ctx, Context): + raise TypeError("cl_ctx must be of type Context, not " + f"{type(cl_ctx)}") + queue = None + actx = None + elif 'queue' in kwargs: + if 'actx' in kwargs: + raise ValueError("At most one of 'actx', 'queue' should " + "be supplied") + queue = kwargs['queue'] + from pyopencl import CommandQueue + if not isinstance(queue, CommandQueue): + raise TypeError("queue must be of type CommandQueue, not " + f"{type(queue)}") + actx = None + elif 'actx' in kwargs: + actx = kwargs['actx'] + if not isinstance(actx, PyOpenCLArrayContext): + raise TypeError("actx must be of type PyOpenCLArrayContext, not " + f"{type(actx)}") + + # now make sure we actually get an actx in kwargs + if actx is None: + if queue is None: + if cl_ctx is None: + from pyopencl import create_some_context + cl_ctx = create_some_context() + from pyopencl import CommandQueue + queue = CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + kwargs['actx'] = actx + + # Set qbx_kwargs defaults + if 'qbx_kwargs' not in kwargs: + degree = function_space.ufl_element().degree() + from sumpy.expansion.level_to_order import FMMLibExpansionOrderFinder + qbx_kwargs = {'fine_order': 4 * degree, + 'qbx_order': degree, + 'fmm_order': None, + 'fmm_level_to_order': + FMMLibExpansionOrderFinder(2**-53, extra_order=1), + 'fmm_backend': "fmmlib", + } + kwargs['qbx_kwargs'] = qbx_kwargs + + # Set op_kwargs defaults + if 'op_kwargs' not in kwargs: + kwargs['op_kwargs'] = {'qbx_forced_limit': None} + + return kwargs + + +def _extract_pytential_operation_kwargs(**kwargs): + """ + Extract kwargs to be passed to :func:`PytentialOperation` + """ + keywords_to_remove = ("cl_ctx", "queue", "actx", "op_kwargs") + pyt_op_kwargs = dict(kwargs) + for key in keywords_to_remove: + pyt_op_kwargs.pop(key, None) - :arg density: A :class:`~firedrake.function.Function` holding the - density. - :returns: The converted density - """ - raise NotImplementedError + return pyt_op_kwargs - def to_firedrake(self, evaluated_potential, out=None): - """ - Convert the evaluated potential from an external library - into a firedrake function - - :arg evaluated_potential: the evaluated potential - :arg out: If not *None*, store the converted potential into this - :class:`firedrake.function.Function`. - :return: *out* if it is not *None*, otherwise a new - :class:`firedrake.function.Function` storing the evaluated - potential - """ - raise NotImplementedError + +def VolumePotential(density, kernel, potential_src_and_tgt, **kwargs): + raise NotImplementedError class PotentialSourceAndTarget: diff --git a/firedrake/potential_evaluation/potentials.py b/firedrake/potential_evaluation/potentials.py index c1dbb4e950..68ba119a71 100644 --- a/firedrake/potential_evaluation/potentials.py +++ b/firedrake/potential_evaluation/potentials.py @@ -1,240 +1,266 @@ """ -User-facing functions to create layer potentials. - -Produces a Potential either through volumential or pytential +Fundamental classes for potential evaluation, including +Potential and a connection to potential-evaluaiton libraries """ +from warnings import warn + from firedrake.functionspaceimpl import WithGeometry +from firedrake.pointwise_operators import AbstractExternalOperator +from firedrake.potential_evaluation import PotentialSourceAndTarget +from firedrake import Function, FunctionSpace, interpolate, \ + VectorFunctionSpace, TensorFunctionSpace -def SingleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): - """ - Evaluate the single layer potential of the density - function convoluted against the kernel on the source - at each target point. - - :arg density: The :mod:`firedrake` density function - :arg kernel: the :class:`sumpy.kernel.Kernel` - :arg potential_src_and_tgt: A - :class:`firedrake.potential_evaluation.PotentialSourceAndTarget` - - Any kwargs not listed below are passed on to the - :class:`firedrake.potential_evaluation.Potential` class - (along with some of the kwargs which are listed below). - - :kwarg function_space: If density is a firedrake function, - this is optional. Otherwise it must be supplied. - If supplied and density is a function, they must - be the same object. - :kwarg cl_ctx: (Optional) a :class:`pyopencl.Context` used - to create a command queue. - At most one of *cl_ctx*, *queue*, *actx* should be included. - If none is included, a *cl_ctx* will be created. - :kwarg queue: (Optional) a :class:`pyopencl.CommandQueue` used - to create an array context. - At most one of *cl_ctx*, *queue*, *actx* should be included. - If none is included, a *cl_ctx* will be created. - :kwarg actx: (Optional) a :class:`meshmode.array_context.PyOpenCLArrayContext` - used for :mod:`pytential` and :mod:`meshmode` computations. - At most one of *cl_ctx*, *queue*, *actx* should be included. - If none is included, a *cl_ctx* will be created. - :kwarg op_kwargs: (Optional) kwargs passed to the construction - of the unbound pytential operator (e.g. {'k': 0.5} for - a :class:`sumpy.kernel.HelmholtzKernel` must be passed - to :func:`pytential.symbolic.primitives.S`), - DEFAULT: - * {'qbx_forced_limit': None} - :kwarg qbx_kwargs: (Optional) *None*, or a dict with arguments to pass to - :class:`pytential.qbx.QBXLayerPotentialSource`. - DEFAULTS: - - fine_order: 4 8 function space degree - - qbx_order: function space degree - - fmm_order: *None* - - fmm_level_to_order: A - :class:`sumpy.level_to_order.FMMLibExpansionOrderFinder` set - to provide tolerance of double machine epsilon - 2**-53 with an extra order of 1. - This allows for taking 1 derivative, but should be - increased for taking more derivatives. - - fmm_backend: "fmmlib" - :kwarg meshmode_connection_kwargs: (Optional) - *None*, or a dict with arguments to pass to - :func:`meshmode.interop.firedrake.build_connection_from_firedrake` - :kwarg warn_if_cg: (Optional) Pass *False* to suppress the warning if - a "CG" space is used. *True* by default - """ - # TODO: Make fmm_level_to_order depend on number of derivatives somehow? - from pytential.symbolic.primitives import S - return _layer_potential(S, density, kernel, potential_src_and_tgt, **kwargs) +class Potential(AbstractExternalOperator): + r""" + This is a function which represents a potential + computed on a firedrake function. -def DoubleLayerPotential(density, kernel, potential_src_and_tgt, **kwargs): - """ - Same as SingleLayerPotential, but evaluates the double layer potential - instead of the single - """ - from pytential.symbolic.primitives import D - return _layer_potential(D, density, kernel, potential_src_and_tgt, **kwargs) + For a firedrake function :math:`u:\Omega\to\mathbb C^n` + with :math:`\Omega \subseteq \mathbb R^m`, + the potential :math:`P(u):\Omega \to\mathbb C^n` is a + convolution of :math:`u` against some kernel function. + More concretely, given a source region :math:`\Gamma \subseteq \Omega` + a target region :math:`\Sigma \subseteq \Omega`, and + a kernel function :math:`K`, we define + .. math:: -def _layer_potential(layer_potential_sym, - density, kernel, potential_src_and_tgt, **kwargs): - """ - Build a layer potential. For usage, see - :func:`SingleLayerPotential` or :func:`DoubleLayerPotential`. + P(u)(x) = \begin{cases} + \int_\Gamma K(x, y) u(y) \,dy & x \in \Sigma + 0 & x \notin \Sigma + \end{cases} """ - kwargs = _validate_layer_potential_args_and_set_defaults( - density, kernel, potential_src_and_tgt, **kwargs) - # build our unbound operator - assert 'op_kwargs' in kwargs # make sure we got op kwargs - op_kwargs = kwargs['op_kwargs'] - from pytential import sym - unbound_op = layer_potential_sym(kernel, - sym.var("density"), - **op_kwargs) - - # make sure we got an actx during validation - assert 'actx' in kwargs - actx = kwargs['actx'] - # make sure we got a function_space during validation - assert 'function_space' in kwargs - # extract Pytential operation kwargs - pyt_op_kwargs = _extract_pytential_operation_kwargs(**kwargs) - # now return the pytential operation as an external operator - from firedrake.potential_evaluation.pytential import PytentialOperation - return PytentialOperation(actx, density, unbound_op, "density", - potential_src_and_tgt, - **pyt_op_kwargs) - - -def _validate_layer_potential_args_and_set_defaults(density, - kernel, - places, - **kwargs): - """ - Validate the arguments for single/double layer potential. + def __init__(self, density, **kwargs): + """ + :arg density: A :mod:`firedrake` :class:`firedrake.function.Function` + or UFL expression which represents the density :math:`u` - Returns dictionary updated with default arguments - """ - # validate density function space - if hasattr(density, 'function_space'): - if not isinstance(density.function_space(), WithGeometry): - raise TypeError("density.function_space() must be of type " - f"WithGeometry, not {type(density.function_space())}") - function_space = kwargs.get('function_space', None) - if function_space is not None: - if function_space is not density.function_space(): - raise ValueError("density.function_space() and function_space" - " must be the same object") - else: - kwargs['function_space'] = density.function_space() - else: - function_space = kwargs.get('function_space', None) - if function_space is None: - raise ValueError("density has no function_space method, so" - "function_space kwarg must be supplied") + :kwarg operator_data: Should have keys + + :kwarg connection: A :class:`PotentialEvaluationLibraryConnection` + :kwarg potential_operator: The external potential evaluation library + bound operator. + :kwarg function_space: the function space + """ + # super + AbstractExternalOperator.__init__(self, density, **kwargs) + + if 'operator_data' not in kwargs: + raise ValueError("Missing kwarg 'operator_data'") + operator_data = kwargs['operator_data'] + + # Get connection & bound op and validate + if 'connection' not in operator_data: + raise ValueError("Missing operator_data 'connection'") + connection = operator_data['connection'] + if not isinstance(connection, PotentialEvaluationLibraryConnection): + raise TypeError("connection must be of type " + "PotentialEvaluationLibraryConnection, not %s." + % type(connection)) + if 'potential_operator' not in operator_data: + raise ValueError("Missing operator_data 'potential_operator'") + + # Get function space and validate aginst bound op + if 'function_space' not in kwargs: + raise ValueError("Missing kwarg 'function_space'") + + function_space = kwargs['function_space'] if not isinstance(function_space, WithGeometry): - raise TypeError("function_space must be of type " - f"WithGeometry, not {type(function_space)}") - function_space = kwargs['function_space'] - - # validate kernel - from sumpy.kernel import Kernel - if not isinstance(kernel, Kernel): - raise TypeError("kernel must be a sumpy.kernel.Kernel, not of " - f"type '{type(kernel)}'.") - # validate src and tgts - from firedrake.potential_evaluation import PotentialSourceAndTarget - if not isinstance(places, PotentialSourceAndTarget): - raise TypeError("potential_src_and_tgt must be a sumpy.kernel.Kernel, " - f"not of type '{type(places)}'.") - # Make sure src is of right dimension - mesh_gdim = places.mesh.geometric_dimension() - mesh_tdim = places.mesh.topological_dimension() - src_tdim = places.get_source_dimension() - # sanity check - assert mesh_gdim - mesh_tdim in [0, 1] - assert mesh_gdim - src_tdim in [0, 1] - assert mesh_tdim in (mesh_gdim, src_tdim) - # now do the real user-input check - if mesh_gdim - src_tdim != 1: - raise ValueError("source of a layer potential must have co-dimension 1," - f" not {mesh_gdim - src_tdim}.") - - # Now handle pyopencl computing contexts and build - # a PyOpenCLArrayContext - from meshmode.array_context import PyOpenCLArrayContext - cl_ctx = None - queue = None - actx = None - if 'cl_ctx' in kwargs: - if 'actx' in kwargs or 'queue' in kwargs: - raise ValueError("At most one of 'actx', 'queue', 'cl_ctx' should " - "be supplied") - cl_ctx = kwargs['cl_ctx'] - from pyopencl import Context - if not isinstance(cl_ctx, Context): - raise TypeError("cl_ctx must be of type Context, not " - f"{type(cl_ctx)}") - queue = None - actx = None - elif 'queue' in kwargs: - if 'actx' in kwargs: - raise ValueError("At most one of 'actx', 'queue' should " - "be supplied") - queue = kwargs['queue'] - from pyopencl import CommandQueue - if not isinstance(queue, CommandQueue): - raise TypeError("queue must be of type CommandQueue, not " - f"{type(queue)}") - actx = None - elif 'actx' in kwargs: - actx = kwargs['actx'] - if not isinstance(actx, PyOpenCLArrayContext): - raise TypeError("actx must be of type PyOpenCLArrayContext, not " - f"{type(actx)}") - - # now make sure we actually get an actx in kwargs - if actx is None: - if queue is None: - if cl_ctx is None: - from pyopencl import create_some_context - cl_ctx = create_some_context() - from pyopencl import CommandQueue - queue = CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) - kwargs['actx'] = actx - - # Set qbx_kwargs defaults - if 'qbx_kwargs' not in kwargs: - degree = function_space.ufl_element().degree() - from sumpy.expansion.level_to_order import FMMLibExpansionOrderFinder - qbx_kwargs = {'fine_order': 4 * degree, - 'qbx_order': degree, - 'fmm_order': None, - 'fmm_level_to_order': - FMMLibExpansionOrderFinder(2**-53, extra_order=1), - 'fmm_backend': "fmmlib", - } - kwargs['qbx_kwargs'] = qbx_kwargs - - # Set op_kwargs defaults - if 'op_kwargs' not in kwargs: - kwargs['op_kwargs'] = {'qbx_forced_limit': None} - - return kwargs - - -def _extract_pytential_operation_kwargs(**kwargs): + raise TypeError("function_space must be of type WithGeometry, " + f"not {type(function_space)}") + + # Make sure density is a member of our function space, if it is + if isinstance(density, Function): + if density.function_space() != function_space: + raise ValueError("density.function_space() must be the " + "same as function_space") + # save attributes + self.density = density + self.connection = connection + self.potential_operator = operator_data['potential_operator'] + + def _eval_potential_operator(self, density, out=None): + """ + Evaluate the potential on the action coefficients and return. + If *out* is not *None*, stores the result in *out* + + :arg density: the density + :arg out: If not *None*, then a :class:`~firedrake.function.Function` + in which to store the evaluated potential + :return: *out* if it is not *None*, otherwise a new + :class:`firedrake.function.Function` storing the evaluated + potential + """ + density_discrete = interpolate(density, self.function_space()) + density = self.connection.from_firedrake(density_discrete) + potential = self.potential_operator(density) + return self.connection.to_firedrake(potential, out=out) + + def _evaluate(self): + """ + Evaluate P(self.density) into self + """ + return self._eval_potential_operator(self.density, out=self) + + def _compute_derivatives(self): + """ + Return a function + Derivative(P, self.derivatives, self.density) + """ + # FIXME : Assumes only one ufl operand, is that okay? + assert len(self.ufl_operands) == 1 + assert self.ufl_operands[0] is self.density + assert len(self.derivatives) == 1 + # Derivative(P, (0,), self.density)(density_star) + # = P(self.density) + if self.derivatives == (0,): + return lambda density_star, out=None: \ + self._eval_potential_operator(self.density, out=out) + # P is linear, so the nth Gateaux derivative + # Derivative(P, (n,), u)(u*) = P(u*) + # + # d/dx (P \circ u)(v) + # = lim_{t->0} (P(u(v+tx)) - P(u)) / t + # \approx lim_{t->0} (P(u(v) + t du/dx - P(u(v))) / t + # = P(du/dx) + # + # So d^n/dx^n( P \circ u ) = P(d^nu/dx^n) + return self._eval_potential_operator + + def _evaluate_action(self, *args): + """ + Evaluate derivatives of layer potential at action coefficients. + i.e. Derivative(P, self.derivatives, self.density) evaluated at + the action coefficients and store into self + """ + assert len(args) == 1 # sanity check + # Either () or (density_star,) + density_star = args[0] + assert len(density_star) in [0, 1] # sanity check + # Evaluate Derivative(P, self.derivatives, self.density) at density_star + if len(density_star) == 0: + return self._evaluate() + + # FIXME: Assumes just taking Jacobian + density_star, = density_star + operator = self._compute_derivatives() + return operator(density_star, out=self) + + def _evaluate_adjoint_action(self, *args): + """ + Operator is self-adjoint, so just evaluate action + """ + return self._evaluate_action(*args) + + def evaluate_adj_component_control(self, x, idx): + derivatives = tuple(dj + int(idx == j) for j, dj in enumerate(self.derivatives)) + dN = self._ufl_expr_reconstruct_(*self.ufl_operands, derivatives=derivatives) + return dN._evaluate_adjoint_action((x.function,)).vector() + + def evaluate_adj_component(self, x, idx): + print(x, type(x)) + raise NotImplementedError + + +class PotentialEvaluationLibraryConnection: """ - Extract kwargs to be passed to :func:`PytentialOperation` + A connection to an external library for potential evaluation """ - keywords_to_remove = ("cl_ctx", "queue", "actx", "op_kwargs") - pyt_op_kwargs = dict(kwargs) - for key in keywords_to_remove: - pyt_op_kwargs.pop(key, None) + def __init__(self, function_space, potential_source_and_target, + warn_if_cg=True): + """ + Initialize self to work on the function space + + :arg function_space: The :mod:`firedrake` function space + (of type :class:`~firedrake.functionspaceimpl.WithGeometry`) on + which to convert to/from. Must be a 'Lagrange' or + 'Discontinuous Lagrange' space. + :arg potential_source_and_target: A :class:`PotentialSourceAndTarget`. + mesh must match that of function_space. + :arg warn_if_cg: If *True*, warn if the space is "CG" + """ + # validate function space + if not isinstance(function_space, WithGeometry): + raise TypeError("function_space must be of type WithGeometry, not " + "%s" % type(function_space)) + + family = function_space.ufl_element().family() + acceptable_families = ['Discontinuous Lagrange', 'Lagrange'] + if family not in acceptable_families: + raise ValueError("function_space.ufl_element().family() must be " + "one of %s, not '%s'" % + (acceptable_families, family)) + if family == 'Lagrange' and warn_if_cg: + warn("Functions in continuous function space will be projected " + "to/from a 'Discontinuous Lagrange' space. Make sure " + "any operators evaluated are continuous. " + "Pass warn_if_cg=False to suppress this warning.") + + # validate potential_source_and_targets + if not isinstance(potential_source_and_target, PotentialSourceAndTarget): + raise TypeError("potential_source_and_targets must be of type " + "PotentialSourceAndTarget, not '%s'." + % type(potential_source_and_target)) + + # make sure meshes match + if potential_source_and_target.mesh is not function_space.mesh(): + raise ValueError("function_space.mesh() and " + "potential_source_and_target.mesh must be the same" + " obejct.") + + # build DG space if necessary + family = function_space.ufl_element().family() + if family == 'Discontinuous Lagrange': + dg_function_space = function_space + elif family == 'Lagrange': + mesh = function_space.mesh() + degree = function_space.ufl_element().degree() + shape = function_space.shape + if shape is None or len(shape) == 0: + dg_function_space = FunctionSpace(mesh, "DG", degree) + elif len(shape) == 1: + dim = shape[0] + dg_function_space = VectorFunctionSpace(mesh, "DG", degree, + dim=dim) + else: + dg_function_space = TensorFunctionSpace(mesh, "DG", degree, + shape=shape) + else: + acceptable_families = ['Discontinuous Lagrange', 'Lagrange'] + raise ValueError("function_space.ufl_element().family() must be " + "one of %s, not '%s'" % + (acceptable_families, family)) + + # store function space and dg function space + self.function_space = function_space + self.dg_function_space = dg_function_space + self.is_dg = function_space == dg_function_space + # store source and targets + self.source_and_target = potential_source_and_target + + def from_firedrake(self, density): + """ + Convert the density into a form acceptable by an bound operation + in an external library - return pyt_op_kwargs + :arg density: A :class:`~firedrake.function.Function` holding the + density. + :returns: The converted density + """ + raise NotImplementedError + def to_firedrake(self, evaluated_potential, out=None): + """ + Convert the evaluated potential from an external library + into a firedrake function -def VolumePotential(density, kernel, potential_src_and_tgt, **kwargs): - raise NotImplementedError + :arg evaluated_potential: the evaluated potential + :arg out: If not *None*, store the converted potential into this + :class:`firedrake.function.Function`. + :return: *out* if it is not *None*, otherwise a new + :class:`firedrake.function.Function` storing the evaluated + potential + """ + raise NotImplementedError From 319a94ba31eb4443804981ce33f7776206c648ee Mon Sep 17 00:00:00 2001 From: benSepanski Date: Sun, 18 Apr 2021 10:58:13 -0500 Subject: [PATCH 20/20] Update volumential.py so it caches the volumential lookups on recreation --- firedrake/potential_evaluation/pytential.py | 50 ++-- firedrake/potential_evaluation/volumential.py | 216 ++++++++++++------ 2 files changed, 171 insertions(+), 95 deletions(-) diff --git a/firedrake/potential_evaluation/pytential.py b/firedrake/potential_evaluation/pytential.py index 0b06916a1f..65482b0592 100644 --- a/firedrake/potential_evaluation/pytential.py +++ b/firedrake/potential_evaluation/pytential.py @@ -1,17 +1,17 @@ import numpy as np -from firedrake.potential_evaluation import \ +from firedrake.potential_evaluation.potentials import \ Potential, PotentialEvaluationLibraryConnection from firedrake import project from firedrake.functionspacedata import cached @cached -def _cached_build_connection_from_firedrake(mesh, key, - actx, - function_space, - restrict_to_boundary, - **meshmode_connection_kwargs): +def _cached_build_meshmode_connection_from_firedrake(mesh, key, + actx, + function_space, + restrict_to_boundary, + **meshmode_connection_kwargs): """ Build connection from firedrake to meshmode and cache on mesh based on key. @@ -30,8 +30,8 @@ def _cached_build_connection_from_firedrake(mesh, key, **meshmode_connection_kwargs) -def _build_connection_from_firedrake(actx, function_space, restrict_to_boundary, - **meshmode_connection_kwargs): +def _build_meshmode_connection_from_firedrake(actx, function_space, restrict_to_boundary, + **meshmode_connection_kwargs): """ Constructs key and invokes _cached_build_connection_from_firedrake """ @@ -40,21 +40,26 @@ def _build_connection_from_firedrake(actx, function_space, restrict_to_boundary, function_space.ufl_element().degree(), restrict_to_boundary, frozenset(meshmode_connection_kwargs.items())) - return _cached_build_connection_from_firedrake(mesh, key, - actx, - function_space, - restrict_to_boundary, - **meshmode_connection_kwargs) + return _cached_build_meshmode_connection_from_firedrake( + mesh, + key, + actx, + function_space, + restrict_to_boundary, + **meshmode_connection_kwargs) class MeshmodeConnection(PotentialEvaluationLibraryConnection): """ - Build a connection into :mod:`meshmode` + A connection into :mod:`meshmode` """ def __init__(self, function_space, potential_source_and_target, actx, warn_if_cg=True, meshmode_connection_kwargs=None): """ + Construct a connection into :mod:`meshmode`, or retrieve + a cached one if the desired connection already exists + For other args see :class:`PotentialEvaluationLibraryConnection` :arg actx: a :class:`meshmode.array_context.PyOpenCLArrayContext` @@ -105,10 +110,10 @@ def __init__(self, function_space, potential_source_and_target, actx, src_bdy = potential_source_and_target.get_source_id() # Build source connection to (or retrieve cached connection from mesh) - src_conn = _build_connection_from_firedrake(actx, - self.dg_function_space, - restrict_to_boundary=src_bdy, - **meshmode_connection_kwargs) + src_conn = _build_meshmode_connection_from_firedrake(actx, + self.dg_function_space, + restrict_to_boundary=src_bdy, + **meshmode_connection_kwargs) # If source is a boundary, build a connection to restrict it to the # boundary @@ -147,10 +152,11 @@ def __init__(self, function_space, potential_source_and_target, actx, tgt_bdy = potential_source_and_target.get_target_id() # Build target connection to (or retrieve cached connection from mesh) - tgt_conn = _build_connection_from_firedrake(actx, - self.dg_function_space, - restrict_to_boundary=tgt_bdy, - **meshmode_connection_kwargs) + tgt_conn = _build_meshmode_connection_from_firedrake( + actx, + self.dg_function_space, + restrict_to_boundary=tgt_bdy, + **meshmode_connection_kwargs) # store computing context self.actx = actx diff --git a/firedrake/potential_evaluation/volumential.py b/firedrake/potential_evaluation/volumential.py index 829d375f80..1ed4f451ce 100644 --- a/firedrake/potential_evaluation/volumential.py +++ b/firedrake/potential_evaluation/volumential.py @@ -1,6 +1,86 @@ -from firedrake.potential_evaluation import PotentialEvaluationLibraryConnection -from firedrake.potential_evaluation.pytential import \ - MeshmodeConnection +from firedrake.functionspacedata import cached +from firedrake.potential_evaluation.potentials import PotentialEvaluationLibraryConnection +from firedrake.potential_evaluation.pytential import MeshmodeConnection + + +@cached +def _cached_build_volumential_connection(mesh, key, actx, + meshmode_discr, + order, + nlevels, + **box_fmm_geometry_factory_kwargs): + # TODO: Fix documentation + """ + Build connections (meshmode -> volumential, volumential -> meshmode) + + :arg mesh: the mesh to cache on + :arg key: the key this connection is looked up by + (function space ufl element, + restrict_to_boundary, + order, + nlevels, + frozenset(box_fmm_geometry_factory_kwargs.items())) + :arg meshmode_discr: the :class:`meshmode.discretization.Discretization` + a connection is being built to + :arg order: Passed to :class:`volumential.geometry.BoxFMMGeometryFactory`. + :arg nlevels: Passed to :class:`volumential.geometry.BoxFMMGeometryFactory`. + :kwarg box_fmm_geometry_factory_kwargs: Extra kwargs passed to + :class:`volumential.geometry.BoxFMMGeometryFactory`. Must + not include :code:`"expand_to_hold_mesh"`:. + + :return: A tuple *(meshmode -> volumential connection, + volumential -> meshmode connection)*. + """ + # Build connection from meshmode into volumential + # (following https://gitlab.tiker.net/xywei/volumential/-/blob/fe2c3e7af355d5c527060e783237c124c95397b5/test/test_interpolation.py#L72 ) # noqa: E501 + from volumential.geometry import ( + BoundingBoxFactory, BoxFMMGeometryFactory) + from volumential.interpolation import ElementsToSourcesLookupBuilder + bbox_fac = BoundingBoxFactory(dim=meshmode_discr.ambient_dim) + # make sure desired kwargs are present/not present + assert "expand_to_hold_mesh" not in box_fmm_geometry_factory_kwargs + # build our BoxFMMGeometryFactory + bboxfmm_fac = BoxFMMGeometryFactory( + cl_ctx=actx.context, + dim=meshmode_discr.ambient_dim, + bbox_getter=bbox_fac, + expand_to_hold_mesh=meshmode_discr.mesh, + **box_fmm_geometry_factory_kwargs) + boxgeo = bboxfmm_fac(actx.queue) + elt_to_src_lookup_fac = ElementsToSourcesLookupBuilder( + actx.context, tree=boxgeo.tree, discr=meshmode_discr) + elt_to_src_lookup, _ = elt_to_src_lookup_fac(actx.queue) + + # Build connection from volumential into meshmode + from volumential.interpolation import LeavesToNodesLookupBuilder + leaves_to_node_lookup_fac = LeavesToNodesLookupBuilder( + actx.context, trav=boxgeo.trav, discr=meshmode_discr) + leaves_to_node_lookup, _ = leaves_to_node_lookup_fac(actx.queue) + return (leaves_to_node_lookup, elt_to_src_lookup) + + +def _build_volumential_connection(actx, + function_space, + order, + nlevels, + **box_fmm_geometry_factory_kwargs): + """ + Constructs key and invokes _cached_build_volumential_connection + """ + mesh = function_space.mesh() + key = (function_space.ufl_element().family(), + function_space.ufl_element().degree(), + order, + nlevels, + frozenset(box_fmm_geometry_factory_kwargs.items())) + return _cached_build_volumential_connection( + mesh, + key, + actx, + function_space, + order, + nlevels, + **box_fmm_geometry_factory_kwargs) class VolumentialConnection(PotentialEvaluationLibraryConnection): @@ -10,9 +90,8 @@ class VolumentialConnection(PotentialEvaluationLibraryConnection): def __init__(self, potential_source_and_target, function_space, - cl_ctx, - queue, - box_fmm_geometry_kwargs, + actx, + box_fmm_geometry_factory_kwargs, reorder_potentials, meshmode_connection_kwargs, warn_if_cg=True): @@ -21,17 +100,21 @@ def __init__(self, The *potential_source_and_target* must have a 3D source, 3D target, and 3D mesh of co-dimension 0. - :arg cl_ctx: A :class:`pyopencl.Context` - :arg queue: A :class:`pyopencl.CommandQueue` - :arg box_fmm_geometry_kwargs: kwargs passed to + Currently, potential_source_and_target must have the + same source and target, and both target "everywhere" + + :arg actx: A :class:`meshmode.array_context.PyOpenCLArrayContext` + :arg box_fmm_geometry_factory_kwargs: arguments passed to :class:`volumential.geometry.BoxFMMGeometryFactory`. - Must include the key *'nlevels'* + Must not contain :code:`"cl_ctx"`, + :code:`"dim"`, :code:`"bbox_getter"`, or :code:`"expand_to_hold_mesh"` + Defaults: + * order: the degree of the firedrake functions pace + * nlevels: 50 :arg reorder_potentials: *True* iff the volumential FMM is reordering potentials. *False* otherwise. - :arg meshmode_connection_kwargs: A dict passes as kwargs - to :func:`meshmode.interop.build_connection_from_firedrake`. - :mod:`meshmode` is used an intermediate level between - :mod:`firedrake` and :mod:`volumential`. + :arg meshmode_connection_kwargs: As used by + :class:`firedrake.potential_evaluation.pytential.MeshmodeConnection` """ # super PotentialEvaluationLibraryConnection.__init__( @@ -50,75 +133,62 @@ def __init__(self, if potential_source_and_target.mesh().topological_dimension() != 3: raise ValueError("mesh must have topological dimension 3") - # FIXME : allow subdomain regions + # FIXME : allow subdomain regions. + # Note that this will require creating + # volumential connections separately for to vs. from if potential_source_and_target.get_source_id() != "everywhere": raise NotImplementedError("subdomain sources not implemented") if potential_source_and_target.get_target_id() != "everywhere": raise NotImplementedError("subdomain targets not implemented") - # validate cl_ctx argument - from pyopencl import Context - if not isinstance(cl_ctx, Context): - raise TypeError("cl_ctx must be of type Context, not " - "%s." % type(cl_ctx)) - - # validate queue argument - from pyopencl import CommandQueue - if not isinstance(queue, CommandQueue): - raise TypeError("queue must be of type CommandQueue, not " - "%s." % type(queue)) + # validate actx argument + from meshmode.array_context import PyOpenCLArrayContext + if not isinstance(actx, PyOpenCLArrayContext): + raise TypeError("actx must be of type PyOpenCLArrayContext, not " + "%s." % type(actx)) # validate reorder_potentials if not isinstance(reorder_potentials, bool): raise TypeError("reorder_potentials must be of type bool, not " "%s." % type(reorder_potentials)) - # validate box_fmm_geometry_kwargs - if not isinstance(box_fmm_geometry_kwargs, dict): - raise TypeError("box_fmm_geometry_kwargs must be of " + # validate box_fmm_geometry_factory_kwargs + if not isinstance(box_fmm_geometry_factory_kwargs, dict): + raise TypeError("box_fmm_geometry_factory_kwargs must be of " "type dict, not '%s'." % - type(box_fmm_geometry_kwargs)) - if 'nlevels' not in box_fmm_geometry_kwargs: - raise ValueError("box_fmm_geometry_kwargs missing required keyword" - " 'nlevels'.") - - # validate meshmode_connection_kwargs - if meshmode_connection_kwargs is None: - meshmode_connection_kwargs = {} - if not isinstance(meshmode_connection_kwargs, dict): - raise TypeError("meshmode_connection_kwargs must be *None* or of " - "type dict, not '%s'." % - type(meshmode_connection_kwargs)) + type(box_fmm_geometry_factory_kwargs)) + for illegal_key in ["cl_ctx", "dim", "bbox_getter", "expand_to_hold_mesh"]: + if illegal_key in box_fmm_geometry_factory_kwargs: + raise ValueError("box_fmm_geometry_factory_kwargs " + + f"must not have key {illegal_key}.") + # set default values in box_fmm_geometry_factory_kwargs and remove them + # from the dict + order = box_fmm_geometry_factory_kwargs.get("order", function_space.degree()) + del box_fmm_geometry_factory_kwargs["order"] + nlevels = box_fmm_geometry_factory_kwargs.get("nlevels", 50) + del box_fmm_geometry_factory_kwargs["nlevels"] # Build intermediate connection into meshmode - from meshmode.interop.firedrake import build_connection_from_firedrake from meshmode.array_context import PyOpenCLArrayContext - actx = PyOpenCLArrayContext(queue) meshmode_connection = MeshmodeConnection( - function_space, - potential_source_and_target, - actx, - warn_if_cg=warn_if_cg, - meshmode_connection_kwargs=meshmode_connection_kwargs) - - # Build connection from meshmode into volumential - # (following https://gitlab.tiker.net/xywei/volumential/-/blob/fe2c3e7af355d5c527060e783237c124c95397b5/test/test_interpolation.py#L72 ) # noqa : E501 - from volumential.geometry import ( - BoundingBoxFactory, BoxFMMGeometryFactory) - from volumential.interpolation import ElementsToSourcesLookupBuilder - dim = function_space.mesh().geometric_dimension() - bboxfmm_fac = BoxFMMGeometryFactory(cl_ctx, dim=dim, - **box_fmm_geometry_kwargs) - boxgeo = bboxfmm_fac(queue) - elt_to_src_lookup_fac = ElementsToSourcesLookupBuilder( - cl_ctx, tree=boxgeo.tree, discr=meshmode_connection.discr) - elt_to_src_lookup, evt = elt_to_src_lookup_fac(queue) - - # Build connection from volumential into meshmode - from volumential.interpolation import LeavesToNodesLookupBuilder - leaves_to_node_lookup_fac = LeavesToNodesLookupBuilder( - cl_ctx, trav=boxgeo.trav, discr=meshmode_connection.discr) - leaves_to_node_lookup, evt = leaves_to_node_lookup_fac(queue) + function_space, + potential_source_and_target, + actx, + warn_if_cg=warn_if_cg, + meshmode_connection_kwargs=meshmode_connection_kwargs) + + # set to hold appropriate mesh + # FIXME: Once we allow subdomains, will may have to supply + # both source and target meshes separately? + box_fmm_geometry_factory_kwargs["expand_to_hold_mesh"] = \ + meshmode_connection.source_to_meshmode_connection.discr.mesh + # build lookup tables + to_volumential_lookup, from_volumential_lookup = \ + _build_volumential_connection(actx, + function_space, + order, + nlevels, + **box_fmm_geometry_factory_kwargs) # Figure out whether or not conversion from volumential -> meshmode # should user order or tree order @@ -126,11 +196,10 @@ def __init__(self, # Store maps for conversion to/from volumential self.meshmode_connection = meshmode_connection - self.to_volumential_lookup = elt_to_src_lookup - self.from_volumential_lookup = leaves_to_node_lookup + self.to_volumential_lookup = to_volumential_lookup + self.from_volumential_lookup = from_volumential_lookup self.from_volumential_order = order # store necessary computing contextx - self.queue = queue self.actx = actx def from_firedrake(self, density): @@ -142,7 +211,7 @@ def from_firedrake(self, density): # pass flattened operand into volumential interpolator from volumential.interpolation import interpolate_from_meshmode volumential_density = \ - interpolate_from_meshmode(self.queue, + interpolate_from_meshmode(self.actx.queue, meshmode_density, self.to_volumential_lookup, order="user") # user order is more intuitive @@ -153,12 +222,13 @@ def to_firedrake(self, evaluated_potential, out=None): from volumential.interpolation import interpolate_to_meshmode from meshmode.dof_array import unflatten meshmode_pot_vals = \ - interpolate_to_meshmode(self.queue, + interpolate_to_meshmode(self.actx.queue, evaluated_potential, self.from_volumential_lookup, order=self.from_volumential_order) + mm_tgt_conn = self.meshmode_connection.target_to_meshmode_connection meshmode_pot_vals = unflatten(self.actx, - self.meshmode_connection.discr, + mm_tgt_conn.discr, meshmode_pot_vals) # get meshmode data back as firedrake function return self.meshmode_connection.to_firedrake(meshmode_pot_vals, out=self)