volumential.nearfield_potential_table

volumential.nearfield_potential_table.DUFFY_NO_FALLBACK_ENV_VAR = 'VOLUMENTIAL_DUFFY_NO_FALLBACK'

Environment variable that turns the batched-to-scalar DuffyRadial fallback into a hard error. Campaign runs set it so that a table which quietly dropped to the (orders-of-magnitude slower, differently converged) scalar per-entry builder cannot be mistaken for a batched build. It is an environment switch rather than a DuffyBuildConfig field on purpose: the build config is hashed into the table-cache fingerprint, so adding a field there would invalidate every cached table and would make an operational strictness policy part of the numerical cache key.

volumential.nearfield_potential_table.DUFFY_BUILD_ROUTINGS = ('batched', 'scalar', 'scalar-adaptive', 'scalar-fallback')

Recognized values of table.build_routing.

volumential.nearfield_potential_table.kernel_global_scaling_const_to_pymbolic(kernel)

kernel’s global scaling constant as a pymbolic expression.

For every kernel Volumential tabulates the constant is a pure number – \(-1/(2\pi)\), \(1/(4\pi)\), \(i/4\) – so it is folded to a literal here. That keeps the generated knl_scaling assignment a constant regardless of how the upstream mappers happen to shape the expression tree, which is what NearFieldInteractionTable._rewrite_complex_exponentials needs in order to see (and split) an imaginary constant instead of an opaque node. Kernels whose scaling still carries a symbol go through the backend-aware mapper unchanged.

class volumential.nearfield_potential_table.ComplexExponentialRewriter(unproven_arg_names=frozenset({}))

Bases: CSECachingMapperMixin, IdentityMapper

Rewrite exp(re + 1j*im) as exp(re) * (cos(im) + 1j*sin(im)).

pyopencl’s pyopencl-complex.h implements cdouble_exp (and the complex trigonometric functions) with the OpenCL sincos(x, &cosx) out-parameter builtin. On the PoCL 7.0 / LLVM 19.1.7 CPU driver that builtin costs roughly 200 ns per call against roughly 1.6 ns for a separate sin(x) plus cos(x) – a ~130x pathology, measured in isolation and not reproduced on a second CPU OpenCL runtime on the same processor. Because the Duffy table quadrature evaluates the kernel at every node, that single builtin made the 3D Helmholtz direct-table build about ten times slower than the otherwise identical real-valued Yukawa build at the same quadrature orders, and accounted for essentially the whole gap.

Splitting the exponential keeps complex-valued kernels on real exp/cos/sin calls, which the driver compiles normally. The rewrite is exp(a + b) = exp(a) exp(b) composed with Euler’s formula, both of which hold for complex a and b, applied to an exact structural split of the exponent (see _split_complex_expression), so it is valid for genuinely complex exponents too – the damped complex-frequency form exp((-a + 1j b) r) included – and not only for the purely imaginary exp(1j k r) of the Helmholtz kernel. Exponents with no complex constant (Yukawa, Laplace) are left untouched and keep their plain real exp.

The rewrite is value-exact but not conditioning-exact once the phase itself may be complex. For \(z = x + \mathrm{i}y\), both \(\cos z\) and \(\sin z\) grow like \(e^{|y|}/2\) while \(e^{\mathrm{i}z}\) decays like \(e^{-y}\), so Euler’s formula turns a decaying exponential into a cancelling difference of two large terms: severe relative error for moderate \(|y|\) and overflow to infinity or NaN beyond that. A complex exponent is not hypothetical – HelmholtzKernel(dim, allow_evanescent=True) declares its wave number k as complex128 – so a phase that is not provably real keeps its cdouble_exp, which evaluates the decaying result directly and stably. The test is _is_known_real, a positive node-by-node proof rather than a search for complex dependencies: the fused Duffy expressions are post-CSE, and _split_complex_expression hands an opaque CommonSubexpression to the real part wholesale, so a phase can be complex without naming a complex argument anywhere the split can see. unproven_arg_names names the kernel arguments that are not provably real; see _kernel_arg_names_not_known_real.

Mixes in the common-subexpression cache so a shared CSE node in the post-CSE expression DAG is visited once rather than once per reference. CSECachingMapperMixin has to come first in the bases: both it and IdentityMapper define map_common_subexpression, so with the other order the MRO picks the uncached one and the cache – and therefore map_common_subexpression_uncached() – is never reached.

map_common_subexpression_uncached(expr, /, *args, **kwargs)

Rewrite the body of a common subexpression that is not yet cached.

CSECachingMapperMixin calls this once per distinct CommonSubexpression node and caches the result; descending with IdentityMapper is what makes the rewrite reach inside the node instead of stopping at it.

map_call(expr, /, *args, **kwargs)

Split an exp of a complex argument into real exp/cos/sin.

Any other call, an exp whose imaginary part is structurally zero, and an exp whose halves are not provably double-precision real are all returned unchanged; see the class docstring for why, and _is_double_precision_real for the last of those tests.

class volumential.nearfield_potential_table.DuffyBuildConfig(radial_rule: str = 'tanh-sinh-fast', regular_quad_order: object = 20, radial_quad_order: object = 61, mp_dps: int = 50, auto_tune_orders: bool = False, auto_tune_samples: int = 5, auto_tune_floor_factor: float = 8.0, auto_tune_candidates: object = None)

Bases: object

Quadrature settings for one Duffy near-field table build.

The single object that carries a build’s accuracy knobs from the caller to NearFieldInteractionTable.build_table_via_duffy_radial(): the radial rule and its order, the regular (angular/tensor) order, the mpmath working precision, and the optional auto-tuning of the two orders. It is frozen, so it can be recorded verbatim in a table’s provenance.

radial_rule: str = 'tanh-sinh-fast'

regular_quad_order: object = 20

radial_quad_order: object = 61

mp_dps: int = 50

auto_tune_orders: bool = False

auto_tune_samples: int = 5

auto_tune_floor_factor: float = 8.0

auto_tune_candidates: object = None

class volumential.nearfield_potential_table.SymmetryReductionDiagnostics(full_entry_count: int, representative_count: int, compression_ratio: float, orbit_size_histogram: tuple, sign_metadata_count: int, negative_scale_count: int, unreduced_payload_bytes: int, reconstructed_payload_bytes: int, metadata_payload_bytes: int, max_reconstruction_error: object = None, l2_reconstruction_error: object = None)

Bases: object

What a table’s symmetry reduction bought, and what it cost.

Reported by NearFieldInteractionTable.get_symmetry_reduction_diagnostics(): the entry counts before and after reduction and their ratio, the orbit-size histogram and sign bookkeeping behind it, the payload sizes of the three representations, and the error of reconstructing the full table from the representatives.

The two error fields are None unless there is finite reference data to reconstruct against. A table that has not been reduced yet supplies its own data for that, so the usual dense call fills them in; an already reduced table has no full array to compare with, and there the caller has to pass reference_data.

full_entry_count: int

representative_count: int

compression_ratio: float

orbit_size_histogram: tuple

sign_metadata_count: int

negative_scale_count: int

unreduced_payload_bytes: int

reconstructed_payload_bytes: int

metadata_payload_bytes: int

max_reconstruction_error: object = None

l2_reconstruction_error: object = None

volumential.nearfield_potential_table.constant_one(x, y=None, z=None)

The constant function one, broadcast to the shape of x.

volumential.nearfield_potential_table.get_laplace(dim)

The free-space Laplace kernel as a plain callable of its components.

2D only: \(-\log r / (2\pi)\), evaluated from the Cartesian components of the displacement so that it can be handed to the table builder as an ordinary Python function.

volumential.nearfield_potential_table.get_cahn_hilliard(dim, b=0, c=0, approx_at_origin=False)

The 2D Cahn-Hilliard kernel as a plain callable of its components.

With \(\lambda_1^2\) and \(\lambda_2^2\) the two roots of \(\lambda^2 - b\lambda + c\), the returned function evaluates

\[-\frac{1}{2\pi(\lambda_1^2 - \lambda_2^2)} \bigl( K_0(\lambda_1 r) - K_0(\lambda_2 r) \bigr),\]

the same expression that volumential.table_manager.CahnHilliardKernel carries symbolically. With approx_at_origin, each \(K_0\) is replaced by a small-argument series with the leading logarithmic term removed analytically.

Distinct positive roots only, and all three inequalities are strict:

\[b > 0, \qquad c > 0, \qquad b^2 > 4c.\]

The first two make both roots of \(\lambda^2 - b\lambda + c\) positive (they are its sum and its product) and the third makes them distinct, which is what the three steps here need: the real numpy.sqrt of each root, a nonzero \(\lambda_1^2 - \lambda_2^2\) in the prefactor, and a nonzero denominator in the stabilized quadratic formula. Each boundary case breaks a different one – \(b^2 = 4c\) collapses the prefactor, \(c = 0\) reaches a \(0/0\), a negative root reaches \(\sqrt{\text{negative}}\) – and all of them produce a callable that returns nan rather than raising, so check the coefficients before tabulating with it. volumential.table_manager.CahnHilliardKernel is not restricted this way: it takes the complex square root and rejects coincident roots outright.

volumential.nearfield_potential_table.get_cahn_hilliard_laplacian(dim, b=0, c=0)

Not implemented; always raises NotImplementedError.

volumential.nearfield_potential_table.sumpy_kernel_to_lambda(sknl, fallback_dim=None, parameter_values=None)

Turn a sumpy kernel into a plain callable of its components.

The table builders evaluate the kernel at quadrature nodes as an ordinary Python function, so the kernel’s symbolic expression – source and target post-processing applied, multiplied by its global scaling constant – is lambdified, with the Hankel and modified-Bessel calls routed to scipy.special. parameter_values substitutes the kernel’s free symbols, such as a wave number, before lambdification; fallback_dim supplies the dimension for a kernel object that carries neither dim nor ambient_dim.

class volumential.nearfield_potential_table.NearFieldInteractionTable(quad_order, method='gauss-legendre', dim=2, kernel_func=None, kernel_type=None, sumpy_kernel=None, build_method='DuffyRadial', source_box_extent=1, dtype=<class 'numpy.float64'>, derive_kernel_func=True, progress_bar=True, **kwargs)

Bases: object

Class for a near-field interaction table.

A near-field interaction table stores precomputed singular integrals on template boxes and supports transforms to actual boxes on lookup. The query process is done through scaling the entries based on actual box sized.

Orientations are ordered counter-clockwise.

A template box is one of [0,1]^dim

property data

The table entries, allocated full of NaN on first access.

Assigning to it drops the orbit representative ids and marks the table as no longer symmetry-reduced, so the value is taken as the full entry array from then on. The setter only coerces the dtype; giving it anything but a one-dimensional array of _full_entry_count() values is a caller error that later entry lookups will discover, not something it rejects.

get_entry_index(source_mode_index, target_point_index, case_id)

Full table entry id of one (source mode, target, case) interaction.

On a symmetry-reduced table this is the id of the entry’s orbit representative, so an interaction and its images under the kernel’s symmetries resolve to one id. decode_index() inverts the unreduced form of the same encoding.

The id addresses the full entry space, which is an index into data only while the table is stored densely. Compact reduced storage (set_reduced_table_data()) keeps the values in a shorter array addressed through reduced_entry_ids, so read a value with get_entry_data() or get_entry_data_for_full_indices() rather than by indexing data with what this returns.

decode_index(entry_id)

This is the inverse function of get_entry_index()

unwrap_mode_index(mode_index)

Split a flat basis-mode index into one index per axis.

The flattening is row-major over quad_order nodes per axis and has to agree with the mesh generator’s node ordering.

get_template_mode(mode_index)

The mode_index-th basis mode of the template box, as a callable.

Template modes are defined on an l_infty circle.

get_mode(mode_index)

normal modes are defined on the source box

get_mode_cheb_coeffs(mode_index, cheb_order)

Cheb coeffs of a mode. The projection process is performed on [0,1]^dim.

get_symmetry_transform(source_mode_index)

Apply proper transforms to map source mode to a reduced region

Returns: - a transform that can be applied on the interaction case vectors connection box centers. - a transform that can be applied to the mode/point indices.

find_target_point(target_point_index, case_index)

Apply proper transforms to find the target point’s coordinate.

Only translations and scalings are allowed in this step, avoiding the indices of quad points to be messed up.

lookup_by_symmetry(entry_id)

Loop up table entry that is mapped to a region where: - k_i <= q/2 in all direction i - k_i’s are sorted in ascending order

Returns the mapped entry_id

get_entry_data(entry_id)

Return the requested entry value from full or orbit-reduced storage.

get_reduced_entry_ids()

Return full entry IDs stored by the current reduced representation.

get_reduced_table_data()

Return (full_entry_ids, values) for stored table entries.

set_reduced_table_data(entry_ids, values)

Store reduced table values compactly by full entry ID.

get_entry_data_for_full_indices(entry_ids)

Return stored values addressed by full table entry IDs.

has_entry_data_for_full_indices(entry_ids)

Return whether all full entry IDs have stored finite values.

compute_table_entry_duffy_radial(entry_id, radial_rule='tanh-sinh-fast', deg_theta=20, radial_quad_order=61, mp_dps=50)

reconstruct_full_table_from_symmetry(canonical_data=None)

Reconstruct full table data from orbit-canonical representatives.

get_symmetry_reduction_diagnostics(reference_data=None)

Return paper-facing symmetry-reduction counts and reconstruction errors.

build_table_via_duffy_radial_batched(queue, radial_rule='tanh-sinh-fast', deg_theta=20, radial_quad_order=61, mp_dps=50, kernel_kwargs=None)

build_table_via_duffy_radial_batched_2d(queue, radial_rule='tanh-sinh-fast', deg_theta=20, radial_quad_order=61, mp_dps=50, kernel_kwargs=None)

compute_nmlz(mode_id)

build_normalizer_table(pool=None, pb=None)

Build normalizers, used for log-scaled kernels, currently only supported in 2D.

build_table_via_duffy_radial(radial_rule='tanh-sinh-fast', regular_quad_order=20, radial_quad_order=61, mp_dps=50, queue=None, cl_ctx=None, build_config=None, auto_tune_orders=False, auto_tune_samples=5, auto_tune_floor_factor=8.0, auto_tune_candidates=None, **kwargs)

build_table(cl_ctx=None, queue=None, build_config=None, **kwargs)

Fill in this table, by way of the Duffy radial builder.

The entry point the table manager calls. Which Duffy path actually runs – and, when the batched one gives way to the scalar one, why – is decided a level down in build_table_via_duffy_radial() and recorded on the table as build_routing and build_fallback_reason.

build_kernel_exterior_normalizer_table(cl_ctx, queue, pool=None, ncpus=None, mesh_order=5, quad_order=10, mesh_size=0.03, remove_tmp_files=True, **kwargs)

Build the kernel exterior normalizer table for fractional Laplacians.

An exterior normalizer for kernel \(G(r)\) and target \(x\) is defined as

\[\int_{B^c} G(\lVert x - y \rVert) dy\]

where \(B\) is the source box \([0, source_box_extent]^dim\).

get_potential_scaler(entry_id, source_box_size=1, kernel_type=None, kernel_power=None)

Returns a helper function to rescale the table entry based on source_box’s actual size (edge length).