Broken FEEC projections on polar domains by alisa-kirkinskaia · Pull Request #576 · pyccel/psydac

alisa-kirkinskaia · 2026-02-18T11:47:02Z

This PR implements broken FEEC projections for the $C^0$ and $C^1$ sequences in 2D (conga_projections.py). The formulas for the projections can be found on pages 19-23 of the arXiv preprint https://arxiv.org/pdf/2505.15996.

Main additions

Broken FEEC polar projections

Create new module psydac.feec.polar.conga_projections.
Add LinearOperator subclasses C0PolarProjection_V0/1/2 acting on the coefficients of scalar- or vector-valued tensor-product splines defined on the logical domain. The tensor-product splines are projected onto the subspaces $\{V_h^0, V_h^1, V_h^2\}$ of pre-polar 0/1/2-form splines with $C^0$ smoothness after the push-forward.u
Add LinearOperator subclasses C1PolarProjection_U0/1/2 acting on the coefficients of scalar- or vector-valued tensor-product splines defined on the logical domain. The tensor-product splines are projected onto the subspaces $\{U_h^0, U_h^1, U_h^2\}$ of pre-polar 0/1/2-form splines with $C^1$ smoothness after the push-forward.
Make sure that the dot methods of the new projection operators correctly execute in parallel for arbitrary domain decompositions in both the angular and radial dimensions.

Tests

Add unit tests for the projections checking:
- the projection property P(P(x)) = P(x)
- tosparse methods against reference matrix operators
- consistency between dot and tosparse methods

Poisson and Maxwell examples

Add 2D Poisson example in script psydac/feec/polar/examples/poisson_2d.py for solving manufactured Poisson problems on polar mapped domains. The script supports disk, target, and Czarny domains, analytical or spline mappings, and different treatments of the polar singularity (polar-spec, polar-std, C0conga, and C1conga).
Add time-dependent 2D Maxwell example for TE (transverse electric) wave in script psydac/feec/polar/examples/maxwell_2d.py. The script provides two field configurations defined in psydac/feec/polar/examples/analytical_solutions.py:
- CircularCavitySolution: Time-harmonic solution of Maxwell's equations in a disk-like domain with
  perfectly conducting walls
- GaussianInitialCondition: localized rotational Gaussian initial condition for the electric field, with the magnetic field initialized from $B_z = \mathrm{curl} E$

Further changes

Use latest version of Igakit (commit dalcinl/igakit@92ee097 of 2026/07/24) which supports NumPy >= 2.4
Use SymPDE version 0.19.3 which fixes a bug in the linearity checks
Fix a bug in eval_field and eval_field_gradient of TensorFemSpace related to floating point round-off at MPI subdomain boundaries. Add a unit test in psydac/fem/tests/test_eval_fields_parallel.py

Additional info

Poisson example
Example of run: mpirun -n 2 python poisson_2d.py -S -n 16 24 -d 2 2 -t disk -D 0.2 -m 'C0conga'
Exact solution, approximate solution and error plot:

TE Maxwell example
Example of run: mpirun -n 2 python maxwell_2d.py -S -n 16 20 -d 2 2 -T 1 -D 0.2 -s 1
Plot of exact solution and approximate solution at final time T = 1:

TODO

…into polar_splines_alisa

- Remove polar_mapping and use_logical_sol options - Change the analytical solution - Avoid Laplacian() for right hand side - Add disk radius as parse argument

The folder 'subprojects/igakit' was accidentally deleted in commit c51b81b

yguclu

Good job addressing my previous review comments, @alisa-kirkinskaia!

I have gone through the whole PR once again

yguclu · 2026-08-24T13:41:32Z

+    @property
+    def T(self):
+        return self.transpose()


The base class LinearOperator already provides an identical implementation for this property. Hence, there is no reason to repeat it here

yguclu · 2026-08-24T13:42:40Z

+    Parameters:
+    -----------
+
+    W1 : VectorFemSpace (former ProductFemSpace)
+         Full tensor product spline space of the 1-forms S^{p1-1, p2} x S^{p1, p2-1}
+
+    transposed : Boolean
+         Switch between P1 and P1 transposed (default is False)


As in a previous comment

Suggested change

Parameters:

-----------

W1 : VectorFemSpace (former ProductFemSpace)

Full tensor product spline space of the 1-forms S^{p1-1, p2} x S^{p1, p2-1}

transposed : Boolean

Switch between P1 and P1 transposed (default is False)

Parameters:

-----------

W1 : VectorFemSpace

Full tensor product spline space of 1-forms S^{p1-1, p2} x S^{p1, p2-1}.

transposed : bool, default=False

Switch between P1 and P1 transposed.

yguclu · 2026-08-24T13:43:15Z

+    @property
+    def T(self):
+        return self.transpose()


The base class LinearOperator already provides an identical implementation for this property. Hence, there is no reason to repeat it here

yguclu · 2026-08-24T13:56:46Z

+    Parameters:
+    -----------
+
+    W1 : VectorFemSpace (ProductFemSpace)
+         Full tensor product spline space of the 1-forms S^{p1-1, p2} x S^{p1, p2-1}
+
+    transposed : Boolean
+         Switch between P1 and P1 transposed (default is False)
+
+    hbc : Boolean
+         Switch on and off the imposition of homogeneous Dirichlet boundary
+         conditions on the tangential (angular) direction (default is False)
+    """
+


See a previous comment about ProductFemSpace. Further, I believe that the homogeneous Dirichlet BCs are imposed along the radial direction for the tangential component of the field. Is that right?

Suggested change

Parameters:

-----------

W1 : VectorFemSpace (ProductFemSpace)

Full tensor product spline space of the 1-forms S^{p1-1, p2} x S^{p1, p2-1}

transposed : Boolean

Switch between P1 and P1 transposed (default is False)

hbc : Boolean

Switch on and off the imposition of homogeneous Dirichlet boundary

conditions on the tangential (angular) direction (default is False)

"""

Parameters:

-----------

W1 : VectorFemSpace

Full tensor product spline space of 1-forms S^{p1-1, p2} x S^{p1, p2-1}.

transposed : bool, default=False

Switch between P1 and P1 transposed.

hbc : bool, default=False

If True, impose homogeneous Dirichlet boundary conditions on the

tangential (angular) component of the field.

"""

yguclu · 2026-08-24T13:57:44Z

+    Parameters:
+    -----------
+
+    W2 : TensorFemSpace
+         Full tensor product spline space of the 2-forms S^{p1-1, p2-1}
+
+    transposed : Boolean
+         Switch between P2 and P2 transposed (default is False)
+    """
+


Suggested change

Parameters:

-----------

W2 : TensorFemSpace

Full tensor product spline space of the 2-forms S^{p1-1, p2-1}

transposed : Boolean

Switch between P2 and P2 transposed (default is False)

"""

Parameters:

-----------

W2 : TensorFemSpace

Full tensor product spline space of 2-forms S^{p1-1, p2-1}.

transposed : bool, default=False

Switch between P2 and P2 transposed.

"""

Co-authored-by: Yaman Güçlü <yaman.guclu@gmail.com>

yguclu

Thanks for the latest changes @alisa-kirkinskaia!

yguclu · 2026-09-14T09:17:52Z

+
+    Attributes
+    ----------
+    logical_bounds : tuple


If these are tuples of Python floats:

Suggested change

logical_bounds : tuple

logical_bounds : tuple[tuple[float, ...], ...]

yguclu · 2026-09-14T09:22:58Z

+    mapping : sympde.topology.mapping.Mapping
+        Mapping used by the solver. It is either the analytical mapping or
+        its spline approximation. Initialized by calling ``build_geometry``.


Does the PSYDAC class SplineMapping currently subclass the SymPDE class Mapping? @campospinto

yguclu · 2026-09-14T09:26:22Z

+    mapping : sympde.topology.mapping.Mapping
+        Mapping used by the solver. It is either the analytical mapping or
+        its spline approximation. Initialized by calling ``build_geometry``.
+
+    analytical_mapping : sympde.topology.mapping.Mapping
+        Original analytical mapping from the logical to the physical domain.
+
+    domain : sympde.topology.Domain
+        Physical domain associated with ``mapping``. Initialized by
+        calling ``build_geometry``.


Are the double tick marks around build_geometry and mapping really needed? I would expect single ticks to be sufficient

yguclu · 2026-09-14T10:09:20Z

+    :code
+    $(\partial^2_{xx} + \partial^2_{yy}) \phi(x,y) = -\rho(x,y)$
+
+    """


Excellent, thanks! Could you please mention the new module psydac.utilities.operators in the file CHANGELOG.md?

yguclu · 2026-09-14T10:18:15Z

Have you checked that this script still runs after removing the local class Laplacian? I think you are missing an import statement

yguclu · 2026-09-14T10:23:45Z

+    ----------
+    solution_name : str
+        Name of the analytical solution to visualize.
+        Must be either ``cavity`` or ``gaussian``.


Are double tick marks needed for generating the documentation?

Alisa Kirkinskaia and others added 30 commits September 16, 2025 15:05

poisson can run in parallel

c3c801d

reorganize

ad8adad

only one geometry file is created

364d64c

minor changes

04b6e8c

update names in examples/poisson_2d_mapping.py

5012856

error in plot_2d_decomposition

51d0592

merge 'devel' into polar_splines_alisa

63f9b65

fix an error after merging

7354e97

error in plot_2d_decomposition

2c28462

Use tight layout

a86faf1

Avoid string warnings

6d1304f

some changes to visualization

39d14f4

Merge branch 'polar_splines_alisa' of https://github.com/pyccel/psydac …

0ca428e

…into polar_splines_alisa

computing average in parallel

9d2a618

fixed errors in parallel average

77cd4e4

Merge branch 'devel' into polar_splines_alisa

e6c8a1a

update plot_2d_decomposition

dbdf3a1

Changes in disk domain

8ce233f

- Remove polar_mapping and use_logical_sol options - Change the analytical solution - Avoid Laplacian() for right hand side - Add disk radius as parse argument

Callable exact solution in logical coordinates for plotting

c4eab7c

Fix right hand side

1026b3d

Fix smoothing methods 'None', 'polar-std' and 'polar-spec'

ad0a0db

Print info

fb2213b

remove trick

529897b

fixed issue with polar-spec and polar-std (cause unknown)

6feacff

fix crash with analytical mapping

6e74dd0

it now works with analytical mapping

3682565

typo

0e7f56e

Merge branch 'devel' into polar_splines_alisa

573501c

initial commit for maxwell

04e9c4b

fixed some errors

bd3e1d5

alisa-kirkinskaia added 3 commits August 6, 2026 14:48

command line argument for analytical solutions, deleted waveTE

17517ae

fix eval_field

c6c1a32

test for eval_field that fails

0a959da

alisa-kirkinskaia marked this pull request as draft August 18, 2026 14:40

alisa-kirkinskaia marked this pull request as ready for review August 18, 2026 14:40

alisa-kirkinskaia marked this pull request as draft August 18, 2026 15:03

alisa-kirkinskaia marked this pull request as ready for review August 18, 2026 15:03

Put back igakit subproject:

b24a65b

The folder 'subprojects/igakit' was accidentally deleted in commit c51b81b

yguclu marked this pull request as draft August 19, 2026 15:08

yguclu marked this pull request as ready for review August 19, 2026 15:08

fix unit test

d330b0e

alisa-kirkinskaia marked this pull request as draft August 20, 2026 10:16

alisa-kirkinskaia marked this pull request as ready for review August 20, 2026 10:16

updated changelog

ac5e398

yguclu requested changes Aug 24, 2026

View reviewed changes

Merge branch 'devel' into polar_splines_alisa

855c5d1

yguclu added the Next Release Must be in next release label Sep 9, 2026

alisa-kirkinskaia and others added 12 commits September 13, 2026 18:11

docstrings for polar_model

1d155d5

docstrings in utils_congapol.py

3a3d2de

docstring for parallel_utils.py

9b294ea

changelog

cbd3e74

spaces

3ff25fd

Update psydac/feec/polar/examples/analytical_solutions.py

f48bc10

Co-authored-by: Yaman Güçlü <yaman.guclu@gmail.com>

Update psydac/feec/polar/examples/analytical_solutions.py

b474cc8

Co-authored-by: Yaman Güçlü <yaman.guclu@gmail.com>

change argument in analytical_solutions.py

605cd21

Update psydac/feec/polar/examples/analytical_solutions.py

ffed36c

Co-authored-by: Yaman Güçlü <yaman.guclu@gmail.com>

move Laplacian class

54d3856

docstring and remove duplicates for poisson

b49ad06

improve update_plot

2796dff

yguclu requested changes Sep 14, 2026

View reviewed changes

	logical_bounds : tuple
	logical_bounds : tuple[tuple[float, ...], ...]

Conversation

alisa-kirkinskaia commented Feb 18, 2026 • edited by yguclu Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Main additions

Broken FEEC polar projections

Tests

Poisson and Maxwell examples

Further changes

Additional info

TODO

Uh oh!

yguclu left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

yguclu left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

alisa-kirkinskaia commented Feb 18, 2026 •

edited by yguclu

Loading