H(curl) basis and CEED_EVAL_CURL

[sebastiangrimberg]

Preliminary support for H(curl). Changes are very similar to #878.

This handles CeedBasis updates for H(curl) and also extends some of the missing H(div) functionality in ceed-preconditioning.c, for example. Any changes to element restriction for high-order Nedelec elements are left to a future PR, as is backend support (first non-tensor cuda/hip-ref, OCCA, and MAGMA, then later on cuda/hip-shared/gen).

Note: I opened this PR from a fork since I didn't have write access to CEED/libCEED, but if it's preferred to create the branch there then I can close and reopen a new PR once I have access.

TODO

• Add tests similar to H(div) basis constructor #878.
• Backend support: So far have only touched ref-basis, ref-operator, and blocked-operator. What do I need to do for CUDA/HIP/MAGMA/OCCA support for non-tensor H(div) and H(curl)? Maybe this is a bigger ask (Non-Tensor Basis Support #350).
  • Note: Left for a later PR
• Diagonal assembly: Does anything special need to be done here? CeedOperatorGetBasisPointer? What about for p-multigrid?
  • Fixes to multigrid and operator assembly completed for CPU, left for later for backends
• Update Julia/Python/Rust bindings. Also Fortran?
  • I updated the Julia bindings which seems to have added a few methods which I didn't touch (CeedOperatorContextGetDoubleRead, CeedOperatorContextRestoreDoubleRead, etc.). I think this might be just because the generator hadn't been run for a past PR, but let me know if I'm doing something incorrectly.

[jedbrown]

Thanks! It looks like maybe your editor formatting configuration is different. Can you have your editor use clangd so the auto-formatter works, or use make format before committing to avoid whitespace noise?

[sebastiangrimberg]

Thanks! It looks like maybe your editor formatting configuration is different. Can you have your editor use clangd so the auto-formatter works, or use make format before committing to avoid whitespace noise?

I used make format but it might be an issue with my version of clang-format (Homebrew uses 15 now vs. Ubuntu 22 default is 14). I'll rerun tomorrow with v14!

[jeremylt]

Our pipeline currently uses 14. I haven't had the chance to update yet

[jedbrown]

Done.
#1157

[jeremylt]

Re some of your points at the top

1. I'm fine adding GPU support in a follow-up PR. That way this PR won't get too big. Adding initial [hip,cuda]/ref support won't be too bad, but [hip,cuda]/[shared,gen] is related to the other issue you pointed out

2. Similarly, I think its good to add Diagonal assembly support in a follow-up PR. The multigrid infrastructure should be largely unaffected, though I'm not sure if our process for automatically building a coarse to fine interpolation works for H(div) and H(curl) elements?

3. Python, Rust, and Julia bindings would be nice. Our Fortran interface is not a high priority but I wouldn't object to adding this there as well. Let me know if you get stuck trying to add any of those and I can help.

[sebastiangrimberg]

Re some of your points at the top

1. I'm fine adding GPU support in a follow-up PR. That way this PR won't get too big. Adding initial [hip,cuda]/ref support won't be too bad, but [hip,cuda]/[shared,gen] is related to the other issue you pointed out

OK, sounds good. I am very unfamiliar with what will be required here, but I don't think it will be overly complex because H(div)/H(curl) interpolation contractions should look the same (or nearly) to H1 gradient.

2. Similarly, I think its good to add Diagonal assembly support in a follow-up PR. The multigrid infrastructure should be largely unaffected, though I'm not sure if our process for automatically building a coarse to fine interpolation works for H(div) and H(curl) elements?

For multigrid, the issue I see is that CeedOperatorMultigridLevelCreate calls CeedBasisCreateProjection when prolongation or restriction operators are requested. This method I think needs fixing for H(div) and H(curl) interpolation.

3. Python, Rust, and Julia bindings would be nice. Our Fortran interface is not a high priority but I wouldn't object to adding this there as well. Let me know if you get stuck trying to add any of those and I can help.

I see for Julia that there is an automated generator method in julia/libCEED.jl/gen which makes this easy. Is there something similar for Python and Rust or do I have to do it by hand?

[jeremylt]

The raw bindings are automatic but the user facing bindings need to be hand-coded in all languages

[sebastiangrimberg]

Hi @jedbrown @jeremylt, I just pushed a few commits with a bit of a reworking of this idea. I realized that a general implementation for non-tensor H^1, H(div) and H(curl) basis operations should be possible with the benefit of simplifying backend code a lot. The general idea is that CEED_EVAL_INTERP in H(div) and H(curl) looks exactly like CEED_EVAL_GRAD for H^1, H(div) CEED_EVAL_DIV looks like H^1 CEED_EVAL_INTERP, and CEED_EVAL_CURL looks again looks like H^1 CEED_EVAL_GRAD(in 3D at least). So, I just tried to use the CeedBasisGetQuadratureNumComponents more extensively to make all of these operations look the same, with the output of CeedBasisApply being a q_comp x num_comp x num_qpts x num_elems tensor (where q_comp = 1 for H^1 interpolation, for example). For now you can check out the CPU implementation in CeedBasisApply_Ref and I'd be curious your thoughts.

In this way, the non-tensor cuda-ref and hip-ref backends should be very easy to port to H(div) and H(curl) with only some minor refactoring, and most of the methods in ceed-preconditioning.c should also require only a bit of modification to work with H(div) and H(curl).

[rezgarshakeri]

Hi @jedbrown @jeremylt, I just pushed a few commits with a bit of a reworking of this idea. I realized that a general implementation for non-tensor H^1, H(div) and H(curl) basis operations should be possible with the benefit of simplifying backend code a lot. The general idea is that CEED_EVAL_INTERP in H(div) and H(curl) looks exactly like CEED_EVAL_GRAD for H^1, H(div) CEED_EVAL_DIV looks like H^1 CEED_EVAL_INTERP, and CEED_EVAL_CURL looks again looks like H^1 CEED_EVAL_GRAD(in 3D at least). So, I just tried to use the CeedBasisGetQuadratureNumComponents more extensively to make all of these operations look the same, with the output of CeedBasisApply being a q_comp x num_comp x num_qpts x num_elems tensor (where q_comp = 1 for H^1 interpolation, for example). For now you can check out the CPU implementation in CeedBasisApply_Ref and I'd be curious your thoughts.

In this way, the non-tensor cuda-ref and hip-ref backends should be very easy to port to H(div) and H(curl) with only some minor refactoring, and most of the methods in ceed-preconditioning.c should also require only a bit of modification to work with H(div) and H(curl).

Not always. When the problem is scalar in H^1, yes CEED_EVAL_GRAD is similar to CEED_EVAL_INTERP in H(div).

[sebastiangrimberg]

Not always. When the problem is scalar in H^1, yes CEED_EVAL_GRAD is similar to CEED_EVAL_INTERP in H(div).

Yep that's correct. My point was simply that with num_comp = 1 we should be able to reuse much of the H^1 gradient kernels (for example) and shouldn't need to write much code nor sacrifice performance.

[rezgarshakeri]

Not always. When the problem is scalar in H^1, yes CEED_EVAL_GRAD is similar to CEED_EVAL_INTERP in H(div).

Yep that's correct. My point was simply that with num_comp = 1 we should be able to reuse much of the H^1 gradient kernels (for example) and shouldn't need to write much code nor sacrifice performance.

Yeah, I got your point, I just wanted to mentioned that is true for scalar case and not using it otherwise.

[jeremylt]

PR to update main to non-RELEASE merged into main

[jeremylt]

I think the only thing that we need before merging is create_hdiv_basis and create_hcurl_basis tests for julia/LibCEED.jl/test/rundevtests.jl. Currently there are no dev tests, but there are some examples for create_h1_basis in ``julia/LibCEED.jl/test/runtests.jl`.

Everything else seems to be all set.

[sebastiangrimberg]

I think the only thing that we need before merging is create_hdiv_basis and create_hcurl_basis tests for julia/LibCEED.jl/test/rundevtests.jl. Currently there are no dev tests, but there are some examples for create_h1_basis in ``julia/LibCEED.jl/test/runtests.jl`.

Everything else seems to be all set.

Done! See b2e3f8e.

[jeremylt]

LGTM