Can I specify the potential on embedded boundaries using a MultiFab instead of an analytical function?

[ShuLiu-Clock]

Dear WarpX Community,

I am developing a 2D RZ code with floating potential on embedded boundaries (EB). The potential values differ across the cut cells (cells that intersect the EB). As far as I know, one way to update the EB potential is by using set_potential_on_eb() with an analytical function. Currently, I update the potential by constructing a lengthy function that includes the calculated potential for each cut cell. The function I use to generate this analytical expression is shown below:

   def function_generate(table):

       rows, cols = np.where(table != 0.0)

       if rows.size == 0:
           return "0.0"

       values = table[rows, cols]
       r_lower = rows * dr
       r_upper = (rows + 1) * dr
       z_lower = cols * dz
       z_upper = (cols + 1) * dz

       parts = [
           f"{value:.6g}*("
           f"sqrt(x*x+y*y)>={rl:.6g} and "
           f"sqrt(x*x+y*y)<{ru:.6g} and "
           f"z>={zl:.6g} and z<{zu:.6g}"
           f")"
           for value, rl, ru, zl, zu in zip(values, r_lower, r_upper, z_lower, z_upper)
       ]

       return "+".join(parts)

The generated function consists of parts like c5*(sqrt(x*x+y*y)>=c1 and sqrt(x*x+y*y)<c2 and z>=c3 and z<c4). While this does correctly update the potential, parsing such a long function is very time-consuming. Therefore, I am looking for alternative approaches.
In the Amrex documentation, I found the following information:

This suggests that a MultiFab can be used to specify Dirichlet values on EB faces. My question is: Can I specify the potential on embedded boundaries using a MultiFab in WarpX? Alternatively, is there a Python interface to access the underlying setEBDirichlet routine?
Any guidance would be greatly appreciated.
Thank you.

[WeiqunZhang]

How many parts do you have? Is the potential time dependent?

One thing you could do it to save sqrt(x*x+y*y) in a local variable. Something like

r = sqrt(x*x+y*y); c5*(c1 <= r < c2 and c3 <= z <= c4) + other_parts_using_r + ...

[WeiqunZhang]

Oh, if there are a few thousand terms in the Parser expression, there is not much we can do to make the Parser fast enough. The amrex function you mentioned above is actually from a different class. Nevertheless, it's easy to add a similar function to the class used by WarpX. But WarpX also needs to be modified to use the new amrex function setEBDirichlet(int lev, MultiFab const& eb_potential).

@RemiLehe @ax3l

[ShuLiu-Clock]

OK. Following the same approach, would it also be possible to add a similar function to modify the potential on the domain boundaries ([lo/hi]_[x/y/z])? This would be helpful. (,,・ω・,,)

[WeiqunZhang]

Just wanted to make sure I understand. You really have a few thousand distinct values, not just a few distinct values on a few thousand EB nodes.

[WeiqunZhang]

The Dirichlet domain boundaries are under the application's control. Whatever you set on the Dirichlet boundary nodes will be used the amrex solver. But I am not sure how you fill that in WarpX. That's a question for @RemiLehe.

[ShuLiu-Clock]

Yes.