Broken FEEC projections on polar domains

[alisa-kirkinskaia]

Implement broken FEEC projections for C0 and C1 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

• 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.
• Add LinearOperator subclasses C1PolarProjection_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 $\{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.
• 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
• Add 2D Poisson example in script psydac/feec/polar/examples/poisson_2d.py
• Add time-dependent 2D Maxwell example for TE (transverse electric) wave in script psydac/feec/polar/examples/maxwell_2d.py

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
• Check if we should 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
• Replace sympy.lambdify with pyccel.lambdify
• Update PR description
• Update CHANGELOG.md

[yguclu]

@alisa-kirkinskaia @campospinto Can you remember why we have no class C1PolarProjection_V2?

Even if it should be identical to C0PolarProjection_V2, for consistency I would at least provide an alias for it, e.g.

C1PolarProjection_V2 = C0PolarProjection_V2

A better solution could be to subclass C0PolarProjection_V2, because it allows us to write a dedicated docstring for it:

class C1PolarProjection_V2(C0PolarProjection_V2):
    """
    CONGA Projector P2 from the full spline space S^{p1-1, p2-1} on logical
    domain to U2, the pre-polar 2-forms splines. The associate matrix
    is square, as in the CONGA approach we keep using the tensor B-spline basis,
    instead of the polar basis of Toshniwal. P2 enforces coefficient relations
    to be in U2.

    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)
    """

    def transpose(self, conjugate=False):
        return C1PolarProjection_V2(self.W2, transposed=not self.transposed)

There is no need to override the other methods, given that they are identical to those of the superclass.

[campospinto]

@alisa-kirkinskaia @campospinto Can you remember why we have no class C1PolarProjection_V2?

Even if it should be identical to C0PolarProjection_V2, for consistency I would at least provide an alias for it, e.g.

C1PolarProjection_V2 = C0PolarProjection_V2

A better solution could be to subclass C0PolarProjection_V2, because it allows us to write a dedicated docstring for it:

class C1PolarProjection_V2(C0PolarProjection_V2):
    """
    CONGA Projector P2 from the full spline space S^{p1-1, p2-1} on logical
    domain to U2, the pre-polar 2-forms splines. The associate matrix
    is square, as in the CONGA approach we keep using the tensor B-spline basis,
    instead of the polar basis of Toshniwal. P2 enforces coefficient relations
    to be in U2.

    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)
    """

    def transpose(self, conjugate=False):
        return C1PolarProjection_V2(self.W2, transposed=not self.transposed)

There is no need to override the other methods, given that they are identical to those of the superclass.

No, if I remember correctly the C0 projection was just used in the simulations but this is much cleaner.

[alisa-kirkinskaia]

@alisa-kirkinskaia @campospinto Can you remember why we have no class C1PolarProjection_V2?

Even if it should be identical to C0PolarProjection_V2, for consistency I would at least provide an alias for it, e.g.

C1PolarProjection_V2 = C0PolarProjection_V2

A better solution could be to subclass C0PolarProjection_V2, because it allows us to write a dedicated docstring for it:

class C1PolarProjection_V2(C0PolarProjection_V2):
    """
    CONGA Projector P2 from the full spline space S^{p1-1, p2-1} on logical
    domain to U2, the pre-polar 2-forms splines. The associate matrix
    is square, as in the CONGA approach we keep using the tensor B-spline basis,
    instead of the polar basis of Toshniwal. P2 enforces coefficient relations
    to be in U2.

    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)
    """

    def transpose(self, conjugate=False):
        return C1PolarProjection_V2(self.W2, transposed=not self.transposed)

There is no need to override the other methods, given that they are identical to those of the superclass.

Yes, it was because the two projections are the same. I added your suggestion in 5ecf431. Do we also need to add C1PolarProjection_V2 to the unit tests?

[yguclu]

@alisa-kirkinskaia @campospinto Can you remember why we have no class C1PolarProjection_V2?
Even if it should be identical to C0PolarProjection_V2, for consistency I would at least provide an alias for it, e.g.

C1PolarProjection_V2 = C0PolarProjection_V2

A better solution could be to subclass C0PolarProjection_V2, because it allows us to write a dedicated docstring for it:

class C1PolarProjection_V2(C0PolarProjection_V2):
    """
    CONGA Projector P2 from the full spline space S^{p1-1, p2-1} on logical
    domain to U2, the pre-polar 2-forms splines. The associate matrix
    is square, as in the CONGA approach we keep using the tensor B-spline basis,
    instead of the polar basis of Toshniwal. P2 enforces coefficient relations
    to be in U2.

    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)
    """

    def transpose(self, conjugate=False):
        return C1PolarProjection_V2(self.W2, transposed=not self.transposed)

There is no need to override the other methods, given that they are identical to those of the superclass.

Yes, it was because the two projections are the same. I added your suggestion in 5ecf431. Do we also need to add C1PolarProjection_V2 to the unit tests?

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

[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