Broken FEEC projections on polar domains

[alisa-kirkinskaia]

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
• Use newer GitHub Actions for repository checkout and Python setup

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

• Add missing class C1PolarProjection_V2 and test it
• Rename C1PolarProjection_V0/1/2 to C1PolarProjection_U0/1/2 (i.e. replace V with U)
• Put analytical solutions in a single module (instead of two) given the amount of duplicated code
• Mention analytical solutions in PR description
• Speed up help messages in example scripts
• Reduce size of output figures so that they fit on a laptop screen
• Provide utility function(s) to handle parallel scripts
• Collect MPI utilities in a single module
• Address TODOs in Maxwell script
• Replace sympy.lambdify with pyccel.lambdify
• Add license header to all new files
• Investigate crash when running maxwell_2d.py in parallel for a certain combination of degree and ncells
• Savage main function in waveTE.py
• Update PR description
• Update CHANGELOG.md

[alisa-kirkinskaia]

Yes, in principle we should test the new class, too. This would be an insurance against future changes to the class. I think we could simply add parametrization to test_PolarProjection_V2

Done in 0c0eae0

[yguclu]

Hi there, I would like to reopen this conversation 😉

This is an example script where we choose the manufactured solution, not a situation where the user can choose an arbitrary function. In any case, even if the solution were more oscillatory, a standard approach would be to increase the number of cells and the spline degree.

Please consider that increasing nquads slows down the matrix assembly a lot. Luckily this is only a 2D example which we do not run with a high number of cells.. but in 3D we would notice a big difference.

In conclusion, I feel that the default number of quadrature points per cell (in each direction) of degree + 1 is appropriate in most cases. Here I would not increase nquads unless strictly necessary.

[alisa-kirkinskaia]

Hi there, I would like to reopen this conversation 😉

This is an example script where we choose the manufactured solution, not a situation where the user can choose an arbitrary function. In any case, even if the solution were more oscillatory, a standard approach would be to increase the number of cells and the spline degree.

Please consider that increasing nquads slows down the matrix assembly a lot. Luckily this is only a 2D example which we do not run with a high number of cells.. but in 3D we would notice a big difference.

In conclusion, I feel that the default number of quadrature points per cell (in each direction) of degree + 1 is appropriate in most cases. Here I would not increase nquads unless strictly necessary.

I see, thanks! f99f5e4