Stokes/Elasticity using biharmonic/Laplace

[isuruf]

• Point sumpy back to main

[alexfikl]

@isuruf I've been using #29 for my stuff. Should I switch to this? Does it contain any fixes or is it just a clean squash on top of main?

[isuruf]

@alexfikl, this PR has only parts of #29 so that it is easier to review. I'm not sure if this PR has all the features that you need. Probably not.

[inducer]

I think it'd be good to have CalculusPatch tests to ensure that PDE rediual for Stokes and elasticity is zero.

[isuruf]

Added a test to check the stokes PDE in a9931e5 and elasticity in bb31471

[inducer]

That's technically not testing the Stokes potential, it's testing pytential-computed derivatives thereof.

[alexfikl]

Added a bunch of nitpicks to pytential.symbolic.stokes (mostly based on how I've been using the new stuff). Many are ignorable, so feel free!

The main thing that bothered me while going through it was the nu_sym everywhere. It feels like that should be handled in elasticity so the user never gets a chance to pass nu_sym=3 from the Stokes wrappers.

[isuruf]

@inducer, @alexfikl, this is ready for another round of reviews.

[alexfikl]

Left a bunch more nitpicks in pde.system_utils this time around :D

This is looking very cool!

[inducer]

That's technically not testing the Stokes potential, it's testing pytential-computed derivatives thereof.

[inducer]

Pull the docs of this into like an "internals" chapter? internals.rst?

[isuruf]

This is ready for another round of reviews

[alexfikl]

Left a bunch of nitpicks while reading this again. I don't think any of them are stoppers of any sort, so go for it when you and @inducer are done! 🚀

Also, thanks for all this Stokes performance work! It's been immensely helpful in the past two-ish years!

[inducer]

Unsubscribing... @-mention or request review once it's ready for a look or needs attention.

[isuruf]

@inducer, ready for a review

[alexfikl]

@inducer Went ahead and squashed this into a single commit to make it easier to rebase/merge from main.

There seem to be some new and exciting test failures in there too (possibly because of a bad rebase). I'll look into those next week or so.

[Copilot]

Pull request overview

This PR expands pytential’s symbolic Stokes/elasticity support by introducing alternative kernel representations (naive/Laplace/biharmonic) and adding infrastructure to rewrite IntG expressions in terms of a chosen base kernel, along with tests and docs updates.

Changes:

• Add new symbolic elasticity wrappers and refactor symbolic Stokes wrappers/operators to support multiple representation methods.
• Introduce pytential.symbolic.pde.system_utils to rewrite layer-potential expressions (including target-to-source transformation handling) and add unit tests for these transformations.
• Add/adjust tests for new utilities (LU solve helper) and update convergence/robustness checks.

Reviewed changes

Copilot reviewed 14 out of 14 changed files in this pull request and generated 6 comments.

Show a summary per file

File	Description
pytential/symbolic/elasticity.py	New elasticity wrapper framework (naive/Laplace/biharmonic) shared with Stokes.
pytential/symbolic/stokes.py	Major refactor of Stokes wrappers/operators to use the new wrapper framework and support method dispatch.
pytential/symbolic/pde/system_utils.py	New machinery for rewriting IntG expressions using a chosen base kernel and derivative relations.
pytential/utils.py	Adds symbolic utilities (chop, forward/back substitution, LU solve helper).
pytential/symbolic/mappers.py	Extends flattening behavior to rewrite IntG densities/kernel args when flattening.
pytential/symbolic/primitives.py	Adjusts Kernel imports/usage to avoid import-time issues and keep ambient_dim inference working.
pytential/qbx/__init__.py	Relaxes density type check to allow scalar densities (but current implementation has a runtime issue).
test/test_tools.py	Adds a unit test for solve_from_lu.
test/test_stokes.py	Expands Stokes tests to cover multiple methods/nu values and adds PDE-based verification tests.
test/test_pde_system_utils.py	New tests for target->source transformation and base-kernel rewriting utilities.
test/test_linalg_skeletonization.py	Makes a zero-comparison more robust via tolerance.
doc/symbolic.rst	Documents new elasticity module and system-utils rewriting utilities.
doc/conf.py	Adds Sphinx intersphinx mapping for ExpressionKernel.
pyproject.toml	Adds mis to typos dictionary.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

The LUdecomposition() call is wrapped in except RewriteFailedError, but LUdecomposition will not raise RewriteFailedError here (and rank-deficient errors from symengine/sympy are different exception types). As written, rank-deficiency will escape the retry logic and crash. Catch the actual exceptions raised by LUdecomposition (or a broader exception type) and mark failed=True so the retry path is exercised.

Suggested change

	try:
	L, U, perm = sym_mat.LUdecomposition()
	except RewriteFailedError:
	# symengine throws an error when rank deficient
	# and sympy returns U with last row zero
	failed = True
	if not sp.USE_SYMENGINE and all(expr == 0 for expr in U[-1, :]):
	L = U = perm = None
	try:
	L, U, perm = sym_mat.LUdecomposition()
	except Exception:
	# symengine throws an error when rank deficient
	# and sympy returns U with last row zero
	failed = True
	if (not failed
	and not sp.USE_SYMENGINE
	and all(expr == 0 for expr in U[-1, :])):

[alexfikl]

This seems fair. The logic in _get_base_kernel_matrix_lu_factorization seems incomplete, but I'm not quite sure what to do there..

[alexfikl]

Spent a bit more time debugging those errors, mainly in test_elasticity_pde. For future reference:

• pytest.param(3, "Laplace", 0.5, False) is failing, but pytest.param(3, "Laplace", 0.5, True) is passing just fine (first uses the Stokeslet/single-layer and the second uses the Stresslet/double-layer).
• All the logic for these two cases is in StressletWrapperTornberg.apply_single_and_double_layer and it does not use anything from system_utils.
• The "pytential" and CalculusPatch versions agree, so all the derivatives are taken correctly.
• I couldn't fine any differences in stokes.py from some older version (that I think worked just fine): https://github.com/inducer/pytential/blob/88cb6ca94d610612518bd577ed0f1e736223bea4/pytential/symbolic/stokes.py.
• The test checks both $\Delta u - p$ and $\nabla \cdot u$. The first one is zero order, but the divergence is "converging" to first order, so likely something wrong in the velocity part (?). Also the $x$ component seems to converge to second-order-ish when setting $\mu = 3$..

EDIT:

• pytest.param(3, "Naive", 0.5, False, marks=pytest.mark.slowtest) seems to be failing pretty much the same as the Laplace version with $\mu = 3$.