Upgrade SymPy

[FrederikSchnack]

Summary

Update SymPDE to support the maintained SymPy release range:

SymPy >= 1.12, < 1.15

To this end, address compatibility problems caused by stricter symbolic construction, traversal, reconstruction, hashing, and sympification behavior in recent SymPy versions.

Closes #106, #197, #198.

Compatibility issues and fixes

1. Incorrect construction of custom Symbol subclasses

ScalarFunction, VectorFunction, and DifferentialForm were constructed through Expr.__new__ or Basic.__new__ instead of SymPy’s Symbol constructor.

Recent SymPy versions expect symbol subclasses to initialize Symbol-specific state, including assumptions, commutativity, and free-symbol behavior. Without this initialization:

• free_symbols could be incorrect;
• assumptions were not initialized consistently;
• scalar and vector functions could have incorrect commutativity;
• equality and hashing did not follow SymPy’s symbol semantics.

The affected objects are now constructed through Symbol.__xnew__.

SymPDE-specific metadata, such as the function space or differential-form degree and dimension, is included through _hashable_content rather than being passed as invalid raw symbol arguments.

Regression tests verify:

• correct scalar and vector commutativity;
• correct free_symbols;
• stable equality and hashing;
• distinct differential forms for different degree or dimension values.

2. Non-SymPy values stored in Basic.args

Several SymPDE expression classes stored ordinary Python objects directly in their structural arguments, including:

• strings;
• integers;
• None;
• lists and tuples;
• dictionaries.

Modern SymPy traverses, hashes, prints, sympifies, and reconstructs every object exposed through Basic.args. Raw Python metadata could therefore cause failures during operations such as:

expr.func(*expr.args)
expr.xreplace(...)
preorder_traversal(expr)

A common structural-argument normalization layer was introduced:

• strings are stored as SymPy Str objects;
• integers are sympified as Integer;
• tuples and lists are stored as SymPy Tuple objects;
• dictionaries are stored as SymPy Dict objects;
• optional None values use a private SymPy-compatible sentinel.

Matching restoration helpers convert these values back to their historical Python representation when accessed through public properties. This preserves the existing SymPDE API while keeping expression trees valid for SymPy.

The normalization is applied to objects including domains, boundaries, interfaces, equations, matrix elements, traces, one-dimensional forms, differential forms, and multipatch mappings.

3. Expressions could not be reconstructed from func and args

SymPy assumes that symbolic expressions can generally be rebuilt using:

expr.func(*expr.args)

This mechanism is used internally by substitution, xreplace, simplification, and other tree transformations.

A number of SymPDE constructors did not satisfy this contract because:

• structural metadata was omitted from args;
• constructors expected keyword-only arguments;
• reconstructed arguments were SymPy containers instead of Python containers;
• subclass constructors interpreted structural arguments as user-facing arguments;
• optional arguments were not converted back from their symbolic representation.

Constructors were updated to accept and restore normalized structural arguments. Reconstruction now retains information such as:

• domain interiors and boundaries;
• mappings and logical domains;
• interface orientations;
• domain connectivity;
• line, square, and cube bounds;
• equation boundary conditions and constraints;
• tensor-form names, axes, and weights;
• multipatch mapping dictionaries.

Regression tests now exercise both reconstruction and substitution:

rebuilt = expr.func(*expr.args)
replaced = expr.xreplace({...})

4. Multipatch mappings used raw dictionaries as symbolic arguments

MultiPatchMapping previously placed a Python dictionary directly in Basic.args and implemented a custom __hash__.

Recent SymPy operations expect structural arguments themselves to be symbolic and reconstructible. This caused problems during hashing, replacement, and rebuilding.

The mapping dictionary is now stored structurally as a SymPy Dict, while the public mappings property continues to return a regular Python dictionary. Hashing is derived from the normalized structural arguments instead of a separate custom implementation.

Tests verify that multipatch mappings survive both reconstruction and xreplace.

5. Function symbols with the same name collapsed across spaces

Function symbols were primarily identified by their printed name. As a result, functions with the same name but belonging to different function spaces could compare as equal and collapse during SymPy simplification.

For example, an expression conceptually equivalent to:

u_from_V + u_from_W

could be simplified as though both functions were the same symbol.

Function spaces now define structural equality from their name, domain, shape, and kind. Scalar and vector function symbols include their function space in _hashable_content.

This ensures that:

• equivalent spaces compare and hash equally;
• different spaces remain distinct;
• functions with the same name in different spaces remain distinct symbols;
• additions do not incorrectly combine functions from different spaces;
• mapped and logical test functions retain the correct domain identity.

6. Constructor postprocessing returned a Python list

The integral multiplication postprocessor wrapped a transformed IntAdd expression in a Python list before passing it back to Mul.

Recent SymPy versions strictly sympify constructor-postprocessor results and reject this list during form-linearity checks.

The postprocessor now returns the transformed IntAdd expression directly:

IntAdd(...)

instead of:

[IntAdd(...)]

7. Legacy SymPy compatibility helpers were outdated

SymPDE carried local copies of compatibility utilities for:

• is_sequence;
• Python 2-style metaclass construction.

These helpers were replaced with maintained functionality:

• use sympy.utilities.iterables.is_sequence;
• use native Python 3 metaclass declarations.

VectorFunction handling remains explicit where its indexed behavior differs from ordinary Python sequences.

The compatibility module itself remains temporarily available because it is still exported by the legacy package namespace. It is removed by the subsequent public-API PR.

8. Domain subclasses were not printed correctly in LaTeX

The custom LaTeX printer only implemented _print_Domain. Some reconstructed domain-related objects are subclasses of BasicDomain rather than Domain, causing modern SymPy’s printer dispatch to miss the custom representation.

A _print_BasicDomain implementation was added so domain interiors and interfaces consistently render using their symbolic names.

Regression tests cover multipatch interiors and interfaces.

Dependency and CI changes

• Raise the minimum supported SymPy version from 1.5 to 1.12.
• Support SymPy releases below 1.15.
• Keep the normal CI matrix on the newest version allowed by the dependency range.
• Add a dedicated Python 3.11 job pinned to SymPy 1.12.
• Run the complete test suite in both minimum- and maximum-supported dependency configurations.

The upper bound remains intentionally below SymPy 1.15 until that release range receives a separate compatibility review.

Validation

The test suite passes with the upgraded symbolic implementation:

170 passed, 2 skipped

New regression coverage includes:

• symbol assumptions and commutativity;
• free_symbols behavior;
• structural equality and hashing;
• functions with identical names in different spaces;
• SymPy-compatible expression arguments;
• func(*args) reconstruction;
• xreplace substitutions;
• multipatch mapping reconstruction;
• domain and interface LaTeX output;
• integral multiplication postprocessing.

[yguclu]

Good job!

Since we have set default_role = 'py:obj' in doc/conf.py, I believe that we can use single backtick marks in all the (new) docstrings.

I have left a few additional comments below

[yguclu]

Why is _restore_optional_argument not used here, as it is done in the functions just below?

Maybe the asymmetry is harmless today, but it is inconsistent and fragile. If a class ever stored a None name, restoring it would silently return the _none object instead of None

[yguclu]

Every usage of this function outside of this module appears to be redundant. Every call site outside of this module is equivalent to one of the restore helpers (restore_optional_argument in domain.py and equation.py, restore_int_argument for axis/ext, restore_container_argument for ornt in topology/basic.py). Please use those.

Additionally, the inline value is _none can be used in restore_optional_argument in this module.

[yguclu]

You have added several helper functions to this file, all with names starting with an underscore. This is usually reserved to private functions, but these are used outside of this module:

• _new_basic
• _restore_optional_argument
• _restore_string_argument
• _restore_int_argument
• _restore_container_argument

Hence, I recommend that you remove the leading _ from their name.

Additionally, I think we can remove function _is_none_argument as I explain in a separate comment below.

[yguclu]

Now that _hashable_content is defined, I believe the __hash__ method is redundant: Basic.__hash__ already computes hash((type(self).__name__,) + self._hashable_content()) and caches it in _mhash. Keeping both means the identity tuple is duplicated and could drift apart in the future. Please remove __hash__ here:

Suggested change

	def __hash__(self):
	return hash((self.name, self.domain, self.shape, self.kind))
	def _hashable_content(self):
	return (self.name, self.domain, self.shape, self.kind)
	def _hashable_content(self):
	return (self.name, self.domain, self.shape, self.kind)

The same applies to Integral.__hash__ in sympde/expr/expr.py (pre-existing, but it would be good to clean it up together).

[yguclu]

Why was an empty line added here?

[yguclu]

Since we are already modifying these lines, let's avoid the useless parentheses:

Suggested change

	assert logical_u.space.domain == domain.logical_domain
	assert(LogicalExpr(dx(u), domain) == dx1(logical_u))
	assert(LogicalExpr(dy(u), domain) == dx2(logical_u))
	assert logical_u.space.domain == domain.logical_domain
	assert LogicalExpr(dx(u), domain) == dx1(logical_u)
	assert LogicalExpr(dy(u), domain) == dx2(logical_u)

[yguclu]

Instead of collecting the reconstruction arguments in *structural_args, could we make the signature follow the order of self.args explicitly?

def __new__(cls, name, interiors=None, boundaries=None, mapping=None,
            logical_domain=None, connectivity=None, *, dim=None):

Here connectivity would accept either a Connectivity object or an iterable of Interface objects. The latter is what is stored in args[5], and Connectivity({str(i.name): i for i in interfaces}) restores it without loss.

With this signature:

• func(*args) works without a separate code path;
• the len(structural_args) != 5 check and the "cannot be mixed" check are no longer needed;
• the signature documents itself, and the order of the arguments is no longer coupled to the _new_basic(...) call at the end of the constructor.

The _restore_*_argument calls stay as they are. The subclasses NCube, Line, Square and Cube already forward to Domain.__new__ positionally, so they are not affected. The only behavior change is that positional calls such as Domain('D', I) become legal, which I think is acceptable.

Please also update the docstrings of the class and of the constructor:

• Both claim that interiors or connectivity must be given. However, Domain('D', connectivity=Connectivity()) raises TypeError: object of type 'NoneType' has no len() (this is already the case on master), while dim alone is accepted and creates a single interior domain.
• mapping and logical_domain must be given together (or not at all).
• interiors may also be a Union, or a nested iterable which is flattened.

[yguclu]

In NCube.__new__ and its subclasses Line, Square and Cube, the if structural_args: branch (taken when SymPy rebuilds the object via func(*args)) forwards directly to Domain.__new__. In this branch the user parameters (dim, min_coords, max_coords, or bounds, bounds1, ...) actually receive the stored Domain arguments (interiors, boundaries, ...), which only works because the positions happen to line up.

More importantly, this branch skips the initialization of obj._coordinates, so the rebuilt object falls back to BasicDomain.coordinates, which is limited to 3 dimensions and returns a list:

Domain	D.coordinates	D.func(*D.args).coordinates
Line('L')	x1	x1
Square('S')	(x1, x2)	[x1, x2]
NCube('N', 4, (0,)*4, (1,)*4)	(x1, x2, x3, x4)	[x1, x2, x3]

The 4D case is wrong. The others only differ in type, but they still show that the rebuilt object is not identical to the original. The existing tests don't catch this because the coordinates are not part of args (so == still returns True), and xreplace returns the original object when nothing is replaced.

Could you:

1. set _coordinates in this "if branch" as well;
2. add a test that compares coordinates before and after func(*args), including a 4D NCube;
3. document in the constructors that, when reconstructing, the positional parameters receive the Domain arguments?

[yguclu]

In the subclasses of Basic1dForm, *structural_args is only needed because name is stored in args, although it is fully determined by the class (Mass is always 'M', Stiffness is always 'S', etc.). When SymPy reconstructs an object, Mass(*args) receives the name in axis, the axis in weight, and the weight in structural_args, which is quite confusing.

Could we store the name as a class attribute instead, and keep only (axis, weight) in args? For example:

class Basic1dForm(AtomicExpr):
    _name = None

    def __new__(cls, axis, weight=S.One):
        return _new_basic(cls, axis, weight)

    @property
    def name(self):
        return self._name

    @property
    def axis(self):
        return int(self.args[0])

    @property
    def weight(self):
        return self.args[1]

    ...

class Mass(Basic1dForm):
    _name = 'M'

class Stiffness(Basic1dForm):
    _name = 'S'

...

Then:

• the __new__ methods of the subclasses and their structural_args branches can be removed entirely;
• func(*args) and xreplace work directly, since the constructor signature matches args;
• the Str conversion of the name is no longer needed;
• equality and hashing still distinguish the subclasses, because Basic includes the type in both (e.g. Mass(0) != Stiffness(0)).

[yguclu]

These old # ... separator are quite confusing.. I suggest replacing them with two blank lines, since we are modifying those lines already

Suggested change

	assert replaced.mappings == {patch: replacement}
	# ...
	def test_mapping_1d():
	assert replaced.mappings == {patch: replacement}
	def test_mapping_1d():

Maybe we could apply the same change to other places in this module?