2D EM solver with EB

[lgiacome]

With this PR we adapt the EB initialization routines to make them work also in 2D (XZ).

Todo

• Modify the initialization routines for edge lengths
• Modify the initialization routines for face areas
• Add a test for the staircased solver
• Add the analysis for the staircased solver
• Add a test for the ECT solver
• Add the analysis for the ECT solver
• Fix particles scraping

This PR is supposed to fix #2390

[lgiacome]

@ax3l it's done :)

[ax3l]

X-ref: particle scraping in 2D, 1D and RZ: #2536
Must have been lost, sorry for the delay.

[ax3l]

Just cc-ing @prkkumar on this.
You are not adding embedded boundary updates yet I think, but I want to make sure we don't collide since this is also improving dimensionality of the code.

[ax3l]

For this large block: is there a way to reuse more code?
It looks to me like we might duplicate a lot of math here.

That duplication makes it hard to maintain in the future, e.g., imagine a fix or extension of any section (which needs to be then remembered to be copied to 4 different geometric implementations that are mostly copied).

[lgiacome]

I'm afraid this won't be possible. The problem is that in 3D the z-direction is stored in the third dimension of the arrays, while in XZ it is stored in the second dimension.
For example, let's say we want to access the top-right neighbor's area. In the 3D code, it is going to be something like Sy(i+1, j, k+1), while in XZ it's going to be Sy(i+1, j+1, k).

I couldn't see an easy way to work around this issue but I am open if you have ideas 😄

[ax3l]

Could you use constants for the index calculation/access and only put those in #ifdefs?

[lgiacome]

You are making a really good point and, in general, I think I should try to reduce the code duplication in this file. We also discussed this with Remi in the past because in 3D we repeat the same operations 3 times and we just change the order of the indices.
The problem is mainly that lines like the following are difficult to generalize

amrex::Real denom = local_avail(0, 1) * Sy(i - 1, j, k) +
                                     local_avail(2, 1) * Sy(i + 1, j, k) +
                                     local_avail(1, 0) * Sy(i, j - 1, k ) +
                                     local_avail(1, 2) * Sy(i, j + 1, k) +
                                     local_avail(0, 0) * Sy(i - 1, j - 1, k) +
                                     local_avail(2, 0) * Sy(i + 1, j - 1, k) +
                                     local_avail(0, 2) * Sy(i - 1, j + 1, k) +
                                     local_avail(2, 2) * Sy(i + 1, j + 1, k);

In principle, I could try to write a sort of decorator for the multifab class that depending on the dimension in consideration accesses the right indices. That might help to not expose this part of code to the choice of indices based on the dimension.

[ax3l]

Yes, we could generally write for instance a constexpr function that does the index permutation at compile time for the chosen dimensionality. That would come with zero runtime overhead.
Additional terms in sums like this can also be #ifdefd out.

[lgiacome]

Sounds like a good strategy. I actually had already started that, but I missed that it should be constexpr
Thanks for the suggestion!

[ax3l]

Thanks a lot!

Formula-wise, I only glimpsed for obvious mistakes and cannot spot any.
Since @RemiLehe is still outstanding, if you like to explore a way to make the 2D/3D code more compact then feel free to push another update. Otherwise ok with me for now :)

[ax3l]

Discussed on Slack: we will already merge this and potential revisions will be addressed in the coming week.