Symbolic operator representation#

Based on pymbolic.

DOF Description#

class pytential.symbolic.dof_desc.DEFAULT_SOURCE[source]#: Symbolic identifier for a default source.

class pytential.symbolic.dof_desc.DEFAULT_TARGET[source]#: Symbolic identifier for a default target.

class pytential.symbolic.dof_desc.QBX_SOURCE_STAGE1[source]#: Symbolic identifier for the Stage 1 discretization of a pytential.qbx.QBXLayerPotentialSource.

class pytential.symbolic.dof_desc.QBX_SOURCE_STAGE2[source]#: Symbolic identifier for the Stage 2 discretization of a pytential.qbx.QBXLayerPotentialSource.

class pytential.symbolic.dof_desc.QBX_SOURCE_QUAD_STAGE2[source]#: Symbolic identifier for the upsampled Stage 2 discretization of a pytential.qbx.QBXLayerPotentialSource.

class pytential.symbolic.dof_desc.GRANULARITY_NODE[source]#: DOFs are per node.

class pytential.symbolic.dof_desc.GRANULARITY_CENTER[source]#: DOFs interleaved per expansion center (two per node, one on each side).

class pytential.symbolic.dof_desc.GRANULARITY_ELEMENT[source]#: DOFs per discretization element.

class pytential.symbolic.dof_desc.DOFDescriptor(geometry: Optional[Hashable] = None, discr_stage: Optional[pytential.symbolic.dof_desc.DiscretizationStages] = None, granularity: Optional[pytential.symbolic.dof_desc.DOFGranularities] = None)[source]#

A data structure specifying the meaning of a vector of degrees of freedom that is handled by pytential (a “DOF vector”). In particular, using geometry, this data structure describes the geometric object on which the (scalar) function described by the DOF vector exists. Using granularity, the data structure describes how the geometric object is discretized (e.g. conventional nodal data, per-element scalars, etc.)

geometry#: An identifier for the geometry on which the DOFs exist. This can be a simple string or any other hashable identifier for the geometric object. The geometric objects are generally subclasses of PotentialSource, TargetBase or Discretization.

discr_stage#: Specific to a pytential.source.LayerPotentialSourceBase, this describes on which of the discretizations the DOFs are defined. Can be one of QBX_SOURCE_STAGE1, QBX_SOURCE_STAGE2 or QBX_SOURCE_QUAD_STAGE2.

granularity#: Describes the level of granularity of the DOF vector. Can be one of GRANULARITY_NODE (one DOF per node), GRANULARITY_CENTER (two DOFs per node, one per side) or GRANULARITY_ELEMENT (one DOF per element).

copy(geometry: ~typing.Optional[~typing.Hashable] = None, discr_stage: ~typing.Optional[pytential.symbolic.dof_desc.DiscretizationStages] = <class 'pytential.symbolic.dof_desc._NoArgSentinel'>, granularity: ~typing.Optional[pytential.symbolic.dof_desc.DOFGranularities] = None) → DOFDescriptor[source]#

to_stage1() → DOFDescriptor[source]#

to_stage2() → DOFDescriptor[source]#

to_quad_stage2() → DOFDescriptor[source]#

pytential.symbolic.dof_desc.as_dofdesc(desc: pytential.symbolic.dof_desc.DOFDescriptorLike) → DOFDescriptor[source]#

class pytential.symbolic.dof_desc.DiscretizationStages#: A Union of all the allowed discretization stages.

class pytential.symbolic.dof_desc.DOFGranularities#: A Union of all the allowed DOF granularity types.

class pytential.symbolic.dof_desc.DOFDescriptorLike#: Types convertible to a DOFDescriptor by as_dofdesc().

Basic objects#

Object types#

Based on the mathematical quantity being represented, the following types of objects occur as part of a symbolic operator representation:

If a quantity is a scalar, it is just a symbolic expression–that is, an element of the set of formal expressions recursively generated by the placeholders (see Placeholders), constants, and arithmetic operations on them (see pymbolic.primitives). Objects of this type are created simply by doing arithmetic on placeholders and scalar constants.
If the quantity is “just a bunch of scalars” (like, say, rows in a system of integral equations), the symbolic representation is an object array. Each element of the object array contains a symbolic expression.

pytools.obj_array.make_obj_array() and pytools.obj_array.flat_obj_array() can help create those.
If it is a geometric quantity (that makes sense without explicit reference to coordinates), it is a pymbolic.geometric_algebra.MultiVector. This can be converted to an object array by calling: pymbolic.geometric_algebra.MultiVector.as_vector().

pyopencl.array.Array and meshmode.dof_array.DOFArray instances hold per-node degrees of freedom (and only those). Such instances do not occur on the symbolic side of pytential at all. They’re only visible either as bound inputs (see pytential.bind()) or outputs of evaluation. Which one is used depends on the meaning of the data being represented. If the data is associated with a Discretization, then DOFArray is used and otherwise Array is used.

Diagnostics#

class pytential.symbolic.primitives.ErrorExpression(message)[source]#

An expression that, if evaluated, causes an error with the supplied message.

__init__(message)[source]#

Placeholders#

pytential.symbolic.primitives.var[source]#: alias of Variable

pytential.symbolic.primitives.make_sym_mv(name, num_components)[source]#

pytential.symbolic.primitives.make_sym_surface_mv(name, ambient_dim, dim, dofdesc=None)[source]#

Functions#

pytential.symbolic.primitives.real#

pytential.symbolic.primitives.imag#

pytential.symbolic.primitives.conj#

pytential.symbolic.primitives.abs#

pytential.symbolic.primitives.sqrt#

pytential.symbolic.primitives.sin#

pytential.symbolic.primitives.cos#

pytential.symbolic.primitives.tan#

pytential.symbolic.primitives.asin#

pytential.symbolic.primitives.acos#

pytential.symbolic.primitives.atan#

pytential.symbolic.primitives.atan2#

pytential.symbolic.primitives.sinh#

pytential.symbolic.primitives.cosh#

pytential.symbolic.primitives.tanh#

pytential.symbolic.primitives.asinh#

pytential.symbolic.primitives.acosh#

pytential.symbolic.primitives.atanh#

pytential.symbolic.primitives.exp#

pytential.symbolic.primitives.log#

Discretization properties#

class pytential.symbolic.primitives.IsShapeClass(shape, dofdesc)[source]#

A predicate that is True if the elements of the discretization have a unique type that matches shape.

shape#: A modepy.Shape subclass.

class pytential.symbolic.primitives.QWeight(dofdesc=None)[source]#: Bare quadrature weights (without Jacobians).

pytential.symbolic.primitives.nodes(ambient_dim, dofdesc=None)[source]#: Return a pymbolic.geometric_algebra.MultiVector of node locations.

pytential.symbolic.primitives.parametrization_derivative(ambient_dim, dim, dofdesc=None)[source]#: Return a pymbolic.geometric_algebra.MultiVector representing the derivative of the reference-to-global parametrization.

pytential.symbolic.primitives.parametrization_derivative_matrix(ambient_dim, dim, dofdesc=None)[source]#: Return a numpy.ndarray representing the derivative of the reference-to-global parametrization.

pytential.symbolic.primitives.pseudoscalar(ambient_dim, dim=None, dofdesc=None)[source]#: Same as the outer product of all parametrization derivative columns.

pytential.symbolic.primitives.area_element(ambient_dim, dim=None, dofdesc=None)[source]#

pytential.symbolic.primitives.sqrt_jac_q_weight(ambient_dim, dim=None, dofdesc=None)[source]#

pytential.symbolic.primitives.normal(ambient_dim, dim=None, dofdesc=None)[source]#: Exterior unit normals.

pytential.symbolic.primitives.mean_curvature(ambient_dim, dim=None, dofdesc=None)[source]#: (Numerical) mean curvature.

pytential.symbolic.primitives.first_fundamental_form(ambient_dim, dim=None, dofdesc=None)[source]#

pytential.symbolic.primitives.second_fundamental_form(ambient_dim, dim=None, dofdesc=None)[source]#: Compute the second fundamental form of a surface. This is in reference to the reference-to-global mapping in use for each element.

Note

Some references assume that the second fundamental form is computed with respect to an orthonormal basis, which this is not.

pytential.symbolic.primitives.shape_operator(ambient_dim, dim=None, dofdesc=None)[source]#

pytential.symbolic.primitives.expansion_radii(ambient_dim, dim=None, granularity=None, dofdesc=None)[source]#

pytential.symbolic.primitives.expansion_centers(ambient_dim, side, dim=None, dofdesc=None)[source]#

pytential.symbolic.primitives.h_max(ambient_dim, dim=None, dofdesc=None)[source]#: Defines a maximum element size in the discretization.

pytential.symbolic.primitives.weights_and_area_elements(ambient_dim, dim=None, dofdesc=None)[source]#: Combines area_element() and QWeight.

Elementary numerics#

class pytential.symbolic.primitives.NumReferenceDerivative(ref_axes=None, operand=None, dofdesc=None)[source]#: An operator that takes a derivative of operand with respect to the the element reference coordinates.

class pytential.symbolic.primitives.NodeSum(operand=None)[source]#: Implements a global sum over all discretization nodes.

class pytential.symbolic.primitives.NodeMax(operand=None)[source]#: Implements a global maximum over all discretization nodes.

class pytential.symbolic.primitives.NodeMin(operand=None)[source]#: Implements a global minimum over all discretization nodes.

class pytential.symbolic.primitives.ElementwiseSum(operand=None, dofdesc=None)[source]#: Returns a vector of DOFs with all entries on each element set to the sum of DOFs on that element.

class pytential.symbolic.primitives.ElementwiseMax(operand=None, dofdesc=None)[source]#: Returns a vector of DOFs with all entries on each element set to the maximum of DOFs on that element.

pytential.symbolic.primitives.integral(ambient_dim, dim, operand, dofdesc=None)[source]#: A volume integral of operand.

class pytential.symbolic.primitives.Ones(dofdesc=None)[source]#: A DOF-vector that is constant one on the whole discretization.

pytential.symbolic.primitives.ones_vec(dim, dofdesc=None)[source]#

pytential.symbolic.primitives.area(ambient_dim, dim, dofdesc=None)[source]#

pytential.symbolic.primitives.mean(ambient_dim, dim, operand, dofdesc=None)[source]#

class pytential.symbolic.primitives.IterativeInverse(expression, rhs, variable_name, extra_vars=None, dofdesc=None)[source]#

Operators#

class pytential.symbolic.primitives.Interpolation(from_dd, to_dd, operand)[source]#: Interpolate quantity from a DOF described by from_dd to a DOF described by to_dd.

pytential.symbolic.primitives.interp(from_dd, to_dd, operand)[source]#

Geometric Calculus (based on Geometric/Clifford Algebra)#

class pytential.symbolic.primitives.Derivative[source]#

Conventional Calculus#

pytential.symbolic.primitives.dd_axis(axis, ambient_dim, operand)[source]#: Return the derivative along (XYZ) axis axis (in ambient_dim-dimensional space) of operand.

pytential.symbolic.primitives.d_dx()#

pytential.symbolic.primitives.d_dy()#

pytential.symbolic.primitives.d_dz()#

pytential.symbolic.primitives.grad_mv(ambient_dim, operand)[source]#: Return the ambient_dim-dimensional gradient of operand as a pymbolic.geometric_algebra.MultiVector.

pytential.symbolic.primitives.grad(ambient_dim, operand)[source]#: Return the ambient_dim-dimensional gradient of operand as a numpy.ndarray.

pytential.symbolic.primitives.laplace(ambient_dim, operand)[source]#

Layer potentials#

class pytential.symbolic.primitives.IntG(target_kernel, source_kernels, densities, qbx_forced_limit, source=None, target=None, kernel_arguments=None, **kwargs)[source]#

\[\int_\Gamma T (\sum S_k G(x-y) \sigma_k(y)) dS_y\]

where \(\sigma_k\) is the k-th density, \(G\) is a Green’s function, \(S_k\) are source derivative operators and \(T\) is a target derivative operator.

target_kernel#

source_kernels#

densities#

qbx_forced_limit#

kernel_arguments#

pytential.symbolic.primitives.int_g_dsource(ambient_dim, dsource, kernel, density, qbx_forced_limit, source=None, target=None, kernel_arguments=None, **kwargs)[source]#

\[\int_\Gamma \operatorname{dsource} \dot \nabla_y \dot g(x-y) \sigma(y) dS_y\]

where \(\sigma\) is density, and dsource, a multivector. Note that the first product in the integrand is a geometric product.

pytential.symbolic.primitives.dsource#: A pymbolic.geometric_algebra.MultiVector.

pytential.symbolic.primitives.int_g_vec(kernel, density, qbx_forced_limit, source=None, target=None, kernel_arguments=None, **kwargs)[source]#

Creates a vector of IntG objects from one kernel with source and target derivatives and maps a vector of densities into a vector of IntG objects.

Historically IntG objects supported only one source kernel and allowed multiple densities to get a vector of objects as a convenience function. Now that IntG objects supports multiple source kernels with one density associated to each source kernel, the previous interface would lead to ambiguity. This function is intended to preserve the “vectorizing” behavior of of the constructor of IntG for use cases where that is preferred.

pytential.symbolic.primitives.S(kernel, density, qbx_forced_limit=<class 'pytential.symbolic.primitives._unspecified'>, source=None, target=None, kernel_arguments=None, **kwargs)[source]#

pytential.symbolic.primitives.Sp(kernel, *args, **kwargs)[source]#

pytential.symbolic.primitives.Spp(kernel, *args, **kwargs)[source]#

pytential.symbolic.primitives.D(kernel, *args, **kwargs)[source]#

pytential.symbolic.primitives.Dp(kernel, *args, **kwargs)[source]#

pytential.symbolic.primitives.normal_derivative(ambient_dim, operand, dim=None, dofdesc=None)[source]#

pytential.symbolic.primitives.tangential_derivative(ambient_dim, operand, dim=None, dofdesc=None)[source]#

“Conventional” Vector Calculus#

pytential.symbolic.primitives.tangential_onb(ambient_dim, dim=None, dofdesc=None)[source]#: Return a matrix of shape (ambient_dim, dim) with orthogonal columns spanning the tangential space of the surface of dofdesc.

pytential.symbolic.primitives.xyz_to_tangential(xyz_vec, dofdesc=None)[source]#

pytential.symbolic.primitives.tangential_to_xyz(tangential_vec, dofdesc=None)[source]#

pytential.symbolic.primitives.project_to_tangential(xyz_vec, dofdesc=None)[source]#

pytential.symbolic.primitives.cross(vec_a, vec_b)[source]#

pytential.symbolic.primitives.n_dot(vec, dofdesc=None)[source]#

pytential.symbolic.primitives.n_cross(vec, dofdesc=None)[source]#

pytential.symbolic.primitives.curl(vec)[source]#

Pretty-printing expressions#

pytential.symbolic.primitives.pretty(expr)[source]#

Binding an operator to a discretization#

pytential.bind(places, expr, auto_where=None)[source]#

Parameters:

places – a pytential.collection.GeometryCollection. Alternatively, any list or mapping that is a valid argument for its constructor can also be used.
auto_where – for simple source-to-self or source-to-target evaluations, find ‘where’ attributes automatically.
expr – one or multiple expressions consisting of primitives form pytential.symbolic.primitives (aka pytential.sym). Multiple expressions can be combined into one object to pass here in the form of a numpy object array

Returns:

a pytential.symbolic.execution.BoundExpression

PDE operators#

Scalar PDEs#

class pytential.symbolic.pde.scalar.L2WeightedPDEOperator(kernel: Kernel, use_l2_weighting: bool)[source]#

\(L^2\)-weighting for scalar IEs.

\(L^2\)-weighting is performed to help with the solution of IEs that yield locally singular densities. It does this by matching the \(\ell^2\)-norm used by the iterative method (e.g. GMRES) with the (approximate) \(L^2\)-norm. Clearly, singular densities might not have finite \(\ell^2\)-norm, hampering convergence of the iterative method. See Bremer, J. On the Nyström discretization of integral equations on planar curves with corners. ACHA, 2011..

kernel#

use_l2_weighting#

get_weight(dofdesc=None)[source]#

Returns:: a symbolic expression for the weights (quadrature weights and area elements) on dofdesc if use_l2_weighting is True and 1 otherwise.

get_sqrt_weight(dofdesc=None)[source]#

Returns:: the square root of get_weight().

get_density_var(name: str)[source]#

Parameters:: name – a string name for the density.
Returns:: a symbolic variable or array (problem dependent) corresponding to the density with the given name.

prepare_rhs(b)[source]#: Modify the right-hand side (e.g. boundary conditions) to match the operator constructed in operator().

representation(u, qbx_forced_limit=None, **kwargs)[source]#

Returns:: a representation for the unknowns based on an integral equation with density u. If qbx_forced_limit denotes an on-surface evaluation, the corresponding jump relations are not added to the representation.

operator(u, **kwargs)[source]#

Returns:: a boundary integral operator with corresponding jump relations that can be used to solve for the density u.

__init__(kernel: Kernel, use_l2_weighting: bool)[source]#

class pytential.symbolic.pde.scalar.DirichletOperator(kernel: Kernel, loc_sign: int, *, alpha: Optional[Number] = None, use_l2_weighting: bool = False, kernel_arguments: Optional[Dict[str, Any]] = None)[source]#

IE operator and field representation for solving Dirichlet boundary value problems with scalar kernels (e.g. LaplaceKernel, HelmholtzKernel, YukawaKernel)

Note

When testing this as a potential matcher, note that it can only access potentials that come from charge distributions having no net charge. (This is true at least in 2D)

Inherits from L2WeightedPDEOperator.

is_unique_only_up_to_constant()[source]#

representation(u, map_potentials=None, qbx_forced_limit=None, **kwargs)[source]#

Parameters:

u – symbolic variable for the density.
map_potentials – a callable that can be used to apply additional transformations on all the layer potentials in the representation, e.g. to take a target derivative.
kwargs – additional keyword arguments passed on to the layer potential constructor.

operator(u, **kwargs)[source]#

Parameters:

u – symbolic variable for the density.
kwargs – additional keyword arguments passed on to the layer potential constructor.

__init__(kernel: Kernel, loc_sign: int, *, alpha: Optional[Number] = None, use_l2_weighting: bool = False, kernel_arguments: Optional[Dict[str, Any]] = None)[source]#

Parameters:

loc_sign – \(+1\) for exterior or \(-1\) for interior problems.
alpha – a complex coefficient with positive imaginary part for the combined-field integral representation (CFIE) for the Helmholtz equation (based on Brakhage and Werner). For other kernels, it does does not offer any benefits.

class pytential.symbolic.pde.scalar.NeumannOperator(kernel: Kernel, loc_sign: int, *, alpha: Optional[Number] = None, use_improved_operator: bool = True, use_l2_weighting: bool = False, kernel_arguments: Optional[Dict[str, Any]] = None)[source]#

IE operator and field representation for solving Neumann boundary value problems with scalar kernels (e.g. LaplaceKernel, HelmholtzKernel, YukawaKernel)

Note

When testing this as a potential matcher, note that it can only access potentials that come from charge distributions having no net charge. (This is true at least in 2D)

Inherits from L2WeightedPDEOperator.

is_unique_only_up_to_constant()[source]#

representation(u, map_potentials=None, qbx_forced_limit=None, **kwargs)[source]#

Parameters:

u – symbolic variable for the density.
map_potentials – a callable that can be used to apply additional transformations on all the layer potentials in the representation, e.g. to take a target derivative.

operator(u, **kwargs)[source]#

Parameters:

u – symbolic variable for the density.
kwargs – additional keyword arguments passed on to the layer potential constructor.

__init__(kernel: Kernel, loc_sign: int, *, alpha: Optional[Number] = None, use_improved_operator: bool = True, use_l2_weighting: bool = False, kernel_arguments: Optional[Dict[str, Any]] = None)[source]#

Parameters:

loc_sign – \(+1\) for exterior or \(-1\) for interior problems.
alpha – a complex coefficient with positive imaginary part for the combined-field integral representation (CFIE) for the Helmholtz equation (based on Brakhage and Werner). For other kernels, it does does not offer any benefits.
use_improved_operator – if True use the least singular operator available. Only used when alpha is not \(0\).

class pytential.symbolic.pde.scalar.BiharmonicClampedPlateOperator(kernel: Kernel, loc_sign: int)[source]#

IE operator and field representation for solving clamped plate Biharmonic equation where,

\[\begin{split}\begin{cases} \Delta^2 u = 0, & \quad \text{ on } D, \\ u = g_1, & \quad \text{ on } \partial D, \\ \mathbf{n} \cdot \nabla u = g_2, & \quad \text{ on } \partial D. \end{cases}\end{split}\]

This operator assumes that the boundary data \((g_1, g_2)\) are represented as column vectors and vertically stacked. For details on the formulation see Problem C in [Farkas1990].

Note

This operator supports only interior problem.

[Farkas1990] (1,2)

Farkas, Peter, Mathematical foundations for fast algorithms for the biharmonic equation, Technical Report 765, Department of Computer Science, Yale University, 1990, PDF.

representation(sigma, map_potentials=None, qbx_forced_limit=None, **kwargs)[source]#

Parameters:

sigma – symbolic variable for the density.
map_potentials – a callable that can be used to apply additional transformations on all the layer potentials in the representation, e.g. to take a target derivative.
kwargs – additional keyword arguments passed on to the layer potential constructor.

operator(sigma, **kwargs)[source]#

Parameters:

u – symbolic variable for the density.
kwargs – additional keyword arguments passed on to the layer potential constructor.

Returns:

the second kind integral operator for the clamped plate problem from [Farkas1990].

__init__(kernel: Kernel, loc_sign: int)[source]#

Parameters:: loc_sign – \(+1\) for exterior or \(-1\) for interior problems.

Maxwell’s equations#

pytential.symbolic.pde.maxwell.get_sym_maxwell_point_source(kernel, jxyz, k)[source]#

Return a symbolic expression that, when bound to a pytential.source.PointPotentialSource will yield a field satisfying Maxwell’s equations.

Uses the sign convention \(\exp(-i \omega t)\) for the time dependency.

This will return an object of six entries, the first three of which represent the electric, and the second three of which represent the magnetic field. This satisfies the time-domain Maxwell’s equations as verified by sumpy.point_calculus.frequency_domain_maxwell().

pytential.symbolic.pde.maxwell.get_sym_maxwell_plane_wave(amplitude_vec, v, omega, epsilon=1, mu=1, dofdesc=None)[source]#

Return a symbolic expression that, when bound to a pytential.source.PointPotentialSource will yield a field satisfying Maxwell’s equations.

Parameters:

amplitude_vec – should be orthogonal to v. If it is not, it will be orthogonalized.
v – a three-vector representing the phase velocity of the wave (may be an object array of variables or a vector of concrete numbers) While v may mathematically be complex-valued, this function is for now only tested for real values.
omega – Accepts the “Helmholtz k” to be compatible with other parts of this module.

Uses the sign convention \(\exp(-1 \omega t)\) for the time dependency.

class pytential.symbolic.pde.maxwell.PECChargeCurrentMFIEOperator(k=None)[source]#

Magnetic Field Integral Equation operator with PEC boundary conditions, under the assumption of no surface charges.

See contrib/notes/mfie.tm in the repository for a derivation.

The arguments loc below decide between the exterior and the interior MFIE. The “exterior” (loc=1) MFIE enforces a zero field on the interior of the integration surface, whereas the “interior” MFIE (loc=-1) enforces a zero field on the exterior.

Uses the sign convention \(\exp(-1 \omega t)\) for the time dependency.

for the sinusoidal time dependency.

j_operator(loc, Jt)[source]#

j_rhs(Hinc_xyz)[source]#

rho_operator(loc, rho)[source]#

rho_rhs(Jt, Einc_xyz)[source]#

scattered_volume_field(Jt, rho, qbx_forced_limit=None)[source]#: This will return an object array of six entries, the first three of which represent the electric, and the second three of which represent the magnetic field. This satisfies the time-domain Maxwell’s equations as verified by sumpy.point_calculus.frequency_domain_maxwell().

Stokes’ equations#

class pytential.symbolic.stokes.StokesletWrapper(dim=None)[source]#

Wrapper class for the StokesletKernel kernel.

This class is meant to shield the user from the messiness of writing out every term in the expansion of the double-indexed Stokeslet kernel applied to the density vector. The object is created to do some of the set-up and bookkeeping once, rather than every time we want to create a symbolic expression based on the kernel – say, once when we solve for the density, and once when we want a symbolic representation for the solution, for example.

The apply() function returns the integral expressions needed for the vector velocity resulting from convolution with the vector density, and is meant to work similarly to calling S() (which is IntG).

Similar functions are available for other useful things related to the flow: apply_pressure(), apply_derivative() (target derivative), apply_stress() (applies symmetric viscous stress tensor in the requested direction).

kernel_dict#: The dictionary allows us to exploit symmetry – that \(S_{01}\) is identical to \(S_{10}\) – and avoid creating multiple expansions for the same kernel in a different ordering.

__init__(dim=None)[source]#

apply(density_vec_sym, mu_sym, qbx_forced_limit)[source]#

Symbolic expressions for integrating Stokeslet kernel.

Returns an object array of symbolic expressions for the vector resulting from integrating the dyadic Stokeslet kernel with variable density_vec_sym.

Parameters:

density_vec_sym – a symbolic vector variable for the density vector.
mu_sym – a symbolic variable for the viscosity.
qbx_forced_limit – the qbx_forced_limit argument to be passed on to IntG.

apply_pressure(density_vec_sym, mu_sym, qbx_forced_limit)[source]#: Symbolic expression for pressure field associated with the Stokeslet.

apply_derivative(deriv_dir, density_vec_sym, mu_sym, qbx_forced_limit)[source]#

Symbolic derivative of velocity from Stokeslet.

Returns an object array of symbolic expressions for the vector resulting from integrating the deriv_dir target derivative of the dyadic Stokeslet kernel with variable density_vec_sym.

Parameters:

deriv_dir – integer denoting the axis direction for the derivative.
density_vec_sym – a symbolic vector variable for the density vector.
mu_sym – a symbolic variable for the viscosity.
qbx_forced_limit – the qbx_forced_limit argument to be passed on to IntG.

apply_stress(density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit)[source]#

Symbolic expression for viscous stress applied to a direction.

Returns a vector of symbolic expressions for the force resulting from the viscous stress

\[-p \delta_{ij} + \mu (\nabla_i u_j + \nabla_j u_i)\]

applied in the direction of dir_vec_sym.

Note that this computation is very similar to computing a double-layer potential with the Stresslet kernel in StressletWrapper. The difference is that here the direction vector is applied at the target points, while in the Stresslet the direction is applied at the source points.

Parameters:

density_vec_sym – a symbolic vector variable for the density vector.
dir_vec_sym – a symbolic vector for the application direction.
mu_sym – a symbolic variable for the viscosity.
qbx_forced_limit – the qbx_forced_limit argument to be passed on to IntG.

class pytential.symbolic.stokes.StressletWrapper(dim=None)[source]#

Wrapper class for the StressletKernel kernel.

This class is meant to shield the user from the messiness of writing out every term in the expansion of the triple-indexed Stresslet kernel applied to both a normal vector and the density vector. The object is created to do some of the set-up and bookkeeping once, rather than every time we want to create a symbolic expression based on the kernel – say, once when we solve for the density, and once when we want a symbolic representation for the solution, for example.

The apply() function returns the integral expressions needed for convolving the kernel with a vector density, and is meant to work similarly to S() (which is IntG).

kernel_dict#: The dictionary allows us to exploit symmetry – that \(T_{012}\) is identical to \(T_{120}\) – and avoid creating multiple expansions for the same kernel in a different ordering.

__init__(dim=None)[source]#

apply(density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit)[source]#

Symbolic expressions for integrating Stresslet kernel.

Returns an object array of symbolic expressions for the vector resulting from integrating the dyadic Stresslet kernel with variable density_vec_sym and source direction vectors dir_vec_sym.

Parameters:

density_vec_sym – a symbolic vector variable for the density vector.
dir_vec_sym – a symbolic vector variable for the direction vector.
mu_sym – a symbolic variable for the viscosity.
qbx_forced_limit – the qbx_forced_limit argument to be passed on to IntG.

apply_pressure(density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit)[source]#: Symbolic expression for pressure field associated with the Stresslet.

apply_derivative(deriv_dir, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit)[source]#

Symbolic derivative of velocity from stresslet.

Returns an object array of symbolic expressions for the vector resulting from integrating the deriv_dir target derivative of the dyadic Stresslet kernel with variable density_vec_sym and source direction vectors dir_vec_sym.

Parameters:

deriv_dir – integer denoting the axis direction for the derivative.
density_vec_sym – a symbolic vector variable for the density vector.
dir_vec_sym – a symbolic vector variable for the normal direction.
mu_sym – a symbolic variable for the viscosity.
qbx_forced_limit – the qbx_forced_limit argument to be passed on to IntG.

apply_stress(density_vec_sym, normal_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit)[source]#

Symbolic expression for viscous stress applied to a direction.

Returns a vector of symbolic expressions for the force resulting from the viscous stress

\[-p \delta_{ij} + \mu (\nabla_i u_j + \nabla_j u_i)\]

applied in the direction of dir_vec_sym.

Parameters:

density_vec_sym – a symbolic vector variable for the density vector.
normal_vec_sym – a symbolic vector variable for the normal vectors (outward facing normals at source locations).
dir_vec_sym – a symbolic vector for the application direction.
mu_sym – a symbolic variable for the viscosity.
qbx_forced_limit – the qbx_forced_limit argument to be passed on to IntG.

class pytential.symbolic.stokes.StokesOperator(ambient_dim, side)[source]#

ambient_dim#

side#

__init__(ambient_dim, side)[source]#

Parameters:

ambient_dim – dimension of the ambient space.
side – \(+1\) for exterior or \(-1\) for interior.

get_density_var(name='sigma')[source]#

Returns:: a symbolic vector corresponding to the density.

prepare_rhs(b, *, mu)[source]#

Returns:: a (potentially) modified right-hand side b that matches requirements of the representation.

operator(sigma)[source]#

Returns:: the integral operator that should be solved to obtain the density sigma.

velocity(sigma, *, normal, mu, qbx_forced_limit=None)[source]#

Returns:: a representation of the velocity field in the Stokes flow.

pressure(sigma, *, normal, mu, qbx_forced_limit=None)[source]#

Returns:: a representation of the pressure in the Stokes flow.

class pytential.symbolic.stokes.HsiaoKressExteriorStokesOperator(*, omega, alpha=None, eta=None)[source]#

Representation for 2D Stokes Flow based on [HsiaoKress1985].

Inherits from StokesOperator.

[HsiaoKress1985] (1,2)

G. C. Hsiao and R. Kress, On an Integral Equation for the Two-Dimensional Exterior Stokes Problem, Applied Numerical Mathematics, Vol. 1, 1985, DOI.

__init__(*, omega, alpha=None, eta=None)[source]#

Parameters:

omega – farfield behaviour of the velocity field, as defined by \(A\) in [HsiaoKress1985] Equation 2.3.
alpha – real parameter \(\alpha > 0\).
eta – real parameter \(\eta > 0\). Choosing this parameter well can have a non-trivial effect on the conditioning.

class pytential.symbolic.stokes.HebekerExteriorStokesOperator(*, eta=None)[source]#

Representation for 3D Stokes Flow based on [Hebeker1986].

Inherits from StokesOperator.

[Hebeker1986]

F. C. Hebeker, Efficient Boundary Element Methods for Three-Dimensional Exterior Viscous Flow, Numerical Methods for Partial Differential Equations, Vol. 2, 1986, DOI.

__init__(*, eta=None)[source]#

Parameters:: eta – a parameter \(\eta > 0\). Choosing this parameter well can have a non-trivial effect on the conditioning of the operator.

Scalar Beltrami equations#

class pytential.symbolic.pde.beltrami.BeltramiOperator(kernel: Kernel, *, dim: Optional[int] = None, precond: str = 'left', kernel_arguments: Optional[Dict[str, Any]] = None)[source]#

Beltrami-type operators on closed surfaces.

The construction of the operators is based on [ONeil2018] extended to take any scalar PDE kernel. However, specific kernels may require additional work to allow for unique solutions. For example, the Laplace-Beltrami equation is only unique up to a constant, so using LaplaceBeltramiOperator is recommended to take this into account. In general, the boundary integral equation can be constructed as

sigma = beltrami.get_density_var("sigma")
bc = beltrami.get_density_var("bc")

rhs = beltrami.prepare_rhs(bc)
solution = beltrami.prepare_solution(sigma)
op = beltrami.operator(sigma)

where prepare_solution() is required to recover the solution from th density due to the different types of preconditioning.

[ONeil2018]

M. O’Neil, Second-Kind Integral Equations for the Laplace-Beltrami Problem on Surfaces in Three Dimensions, Advances in Computational Mathematics, Vol. 44, pp. 1385-1409, 2018, DOI.

dim#

__init__(kernel: Kernel, *, dim: Optional[int] = None, precond: str = 'left', kernel_arguments: Optional[Dict[str, Any]] = None) → None[source]#

get_density_var(name: str = 'sigma') → Variable[source]#

Returns:: a symbolic expression for the density.

prepare_solution(sigma: Variable) → Variable[source]#

Returns:: an expression for the solution to the original Beltrami equation based on the density sigma and the type of preconditioning used in the operator.

prepare_rhs(b: Variable) → Variable[source]#

Returns:: a modified expression for the right-hands ide b based on the preconditiong used in the operator.

operator(sigma: Variable, mean_curvature: Optional[Variable] = None, **kwargs) → Variable[source]#

Parameters:: mean_curvature – an expression for the mean curvature that can be used in the construction of the operator.
Returns:: a Fredholm integral equation of the second kind for the Beltrami PDE with the unknown density sigma.

class pytential.symbolic.pde.beltrami.LaplaceBeltramiOperator(ambient_dim, *, dim: Optional[int] = None, precond: str = 'left')[source]#

Laplace-Beltrami operator on a closed surface \(\Sigma\)

\[-\Delta_\Sigma u = b\]

Inherits from BeltramiOperator.

__init__(ambient_dim, *, dim: Optional[int] = None, precond: str = 'left') → None[source]#

class pytential.symbolic.pde.beltrami.YukawaBeltramiOperator(ambient_dim: int, *, dim: Optional[int] = None, precond: str = 'left', yukawa_k_name: str = 'k')[source]#

Yukawa-Beltrami operator on a closed surface \(\Sigma\)

\[-\Delta_\Sigma u + k^2 u = b\]

Inherits from BeltramiOperator.

__init__(ambient_dim: int, *, dim: Optional[int] = None, precond: str = 'left', yukawa_k_name: str = 'k') → None[source]#

class pytential.symbolic.pde.beltrami.HelmholtzBeltramiOperator(ambient_dim: int, *, dim: Optional[int] = None, precond: str = 'left', helmholtz_k_name: str = 'k')[source]#

Helmholtz-Beltrami operator on a closed surface \(\Sigma\)

\[-\Delta_\Sigma u - k^2 u = b\]

Inherits from BeltramiOperator.

__init__(ambient_dim: int, *, dim: Optional[int] = None, precond: str = 'left', helmholtz_k_name: str = 'k') → None[source]#

Internal affairs#

Mappers#

How a symbolic operator gets executed#

class pytential.symbolic.execution.BoundExpression(places, sym_op_expr)[source]#

An expression readied for evaluation by binding it to a GeometryCollection.

cost_per_stage(calibration_params, **kwargs)[source]#

Parameters:: calibration_params – either a dict returned by estimate_kernel_specific_calibration_params, or a str “constant_one”.
Returns:: a dict mapping from instruction to per-stage cost. Each per-stage cost is represented by a dict mapping from the stage name to the predicted time.

cost_per_box(calibration_params, **kwargs)[source]#

Parameters:: calibration_params – either a dict returned by estimate_kernel_specific_calibration_params, or a str “constant_one”.
Returns:: a dict mapping from instruction to per-box cost. Each per-box cost is represented by a numpy.ndarray or pyopencl.array.Array of shape (nboxes,), where the ith entry represents the cost of all stages for box i.

scipy_op(actx: PyOpenCLArrayContext, arg_name, dtype, domains=None, **extra_args)[source]#

Parameters:: domains – a list of discretization identifiers or None values indicating the domains on which each component of the solution vector lives. None values indicate that the component is a scalar. If domains is None, DEFAULT_TARGET is required to be a key in places.
Returns:: An object that (mostly) satisfies the scipy.sparse.linalg.LinearOperator protocol, except for accepting and returning pyopencl.array.Array arrays.

eval(context=None, timing_data=None, array_context: Optional[PyOpenCLArrayContext] = None)[source]#

Evaluate the expression in self, using the input variables given in the dictionary context.

Parameters:

timing_data – A dictionary into which timing data will be inserted during evaluation. (experimental)
array_context – only needs to be supplied if no instances of DOFArray with a PyOpenCLArrayContext are supplied as part of context.

Returns:

the value of the expression, as a scalar, array or an arraycontext.ArrayContainer of these.

__call__(*args, **kwargs)[source]#

Evaluate the expression in self, using the input variables given in the dictionary context.

Returns:: the value of the expression, as a scalar, meshmode.dof_array.DOFArray, or an object array of these.

places#

Created by calling pytential.bind().