Out-of-bounds error in 2D EB

[WeiqunZhang]

While working on #2340, I observed an out-of-bounds error in WarpX::ComputeEdgeLengths if WarpX (2D and with EB) is compiled in debug mode. The issue is even in 2D there are 3 components in m_edge_lengths. So the code like below is incorrect in a number of ways.

    auto const &flags = eb_fact.getMultiEBCellFlagFab();
    auto const &edge_centroid = eb_fact.getEdgeCent();   
    for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi){
        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim){
            auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi);                                                                                                                    
            // ...
            } else {
                auto const &edge_cent = edge_centroid[idim]->const_array(mfi);
                // ...
            }
        }
    }

First, only two out of the three components of m_edge_lengths are initialized here. The second issue is that the second dimension in amrex is the third dimension in 2D WarpX. So for idim = 1, edge_cent from amrex and m_edge_lengths[][1] are for different edges.

[WeiqunZhang]

If the m_edge_lengths etc. are only used by the ECT solver, maybe we do not need to allocate and initialize them.

[lgiacome]

Hi Weiqun,

Alright, I see that this won't work in 2D. I also have to say that the ECT, for now, is only implemented in 3D, but it should be easy to adapt it to work also in 2D and it is already on my TODO list. Maybe for the moment, it would be better to do this initialization process only when working in 3D since 2D and electromagnetics are not really supported yet.

If the m_edge_lengths etc. are only used by the ECT solver, maybe we do not need to allocate and initialize them.

These mfabs are also used by the Yee solver with staircased embedded boundaries, which is useful when modeling simple "square" geometries.

[ax3l]

@lge0303 Thank you for working on those! :)

Can you please add runtime asserts for the things that are not supported as long as they are in-development (RZ and 2D EB/RZ algorithms that are not yet implemented, for instance). That way, it is easier for users/other devs to see (self-documenting behavior).

[PhilMiller]

Rather than asserts, could you please add the appropriate no-op returns to the code, so that the built in calls to these routines don't break code that's not relying on them?

[WeiqunZhang]

Since 2d does not currently work anyway, we could return immediately for 2d as a stopgap.

[ax3l]

X-ref: #2401