I think I discovered incorrect behavior of the MLEBABecLap operator in case of 2D with cylindrical coordinates

[git-user-1p155]

example_1_with_charge_and_no_eb.txt
example_2_with_eb_and_no_charge.txt

Dear AMReX developers,

I think I discovered incorrect behavior of the MLEBABecLap operator. This incorrect behavior is present in the 2D case with cylindrical coordinates, while it is not present when Cartesian coordinates are used. I am providing you with two small programs which manifest this incorrect behavior. I have tested these programs on the CPU, and I do not know whether they would work on a GPU, as I have not yet learned to use GPU functionality. I have tested these programs on the AMReX version from the 21. of April this year.

In the first program, the geometry is regular (no EB is included), and I compare the results of solving the Poisson equation for a simple charge distribution, by using the three linear operators: MLEBAbecLap, MLABecLaplacian, and MLPoisson. All three operators should produce the same result of the Poisson equation. However, in this program the result of using MLEBAbecLap is different from the results of using MLAabecLaplacian and MLPoisson when 2D cylindrical coordinates are used, while all three results are the same when Cartesian coordinates are used (in both 2D and 3D).

In the second program, there is an embedded boundary representing a rod electrode and no charge is included. In this program, the result of using MLEBAbecLap gives different field enhancement in the case of 2D with cylindrical coordinates (axisymmetric case) and the 3D case. This result is incorrect, as the problem is axially symmetric, which means that the axisymmetric case (2D with cylindrical coordinates) should give the same field enhancement as the 3D case. In addition, the obtained resulting electric field in axisymmetric case has practically the same field enhancement around EB as in the 2D case with Cartesian coordinates. This is another indicator that the solution in case of 2D with cylindrical coordinates is incorrect.

I have checked the results of the second program by using the Afivo-streamer open-source code:

https://github.com/MD-CWI/afivo-streamer

The Cartesian 2D result agrees well with the corresponding result from the Afivo-streamer code (when the same embedded boundary and the same applied electric field and domain size are used):

https://github.com/MD-CWI/afivo-streamer/blob/master/programs/standard_2d/streamer_2d_electrode.cfg

The Cartesian 3D result agrees well with the axisymmetric result from the Afivo-streamer code:

https://github.com/MD-CWI/afivo-streamer/blob/master/programs/standard_2d/streamer_cyl_electrode.cfg

I have not done comparison with the 3D Afivo-streamer code, although it can be done.

However, the result of using MLEBAbecLap operator in the second program in 2D with cylindrical coordinates does not agree with the Afivo-streamer axisymmetric result.

Is this incorrect behavior which I have described really a problem with the MLEBABecLap operator or am I doing something wrong?
I expect that the MLEBABecLap operator should read the correct coordinate type from the Geometry objects which are provided in the constructor of the MLEBABecLap object (const Vector< Geometry > &a_geom). I have checked the Doxygen documentation for the MLEBABecLap object, and I have not found a special option for setting the coordinate type.

[WeiqunZhang]

Unfortunately we only support Cartesian coordinates with EB.

[git-user-1p155]

Thank you for the reply! I will have to use 3D Cartesian coordinates then.

[WeiqunZhang]

What kind of embedded boundaries do you have in 2D cyclindrical coordinates?

[git-user-1p155]

Sorry for the late reply, I was not expecting your second message, and I saw the GitHib notification when I was checking my email a few minutes ago.

In the code, which I am developing for work, I have a union of a cylinder and a semisphere. The same embedded boundary is used in the second example that is attached in my first message (example_2_with_eb_and_no_charge). This embedded boundary represents a rod electrode, which creates an electric field enhancement in an ionized gas. An example is shown below:

////////////////////////////////////////////////////////////////////////////////
// defining the embedded boundary representing the rod electrode
amrex::Real radius = 1.0e-3;
amrex::Real x_location = (which_geom == 1)? 0.0:domain_size/2.0;
amrex::Real length = 2.4e-3;
amrex::Real sphere_center = domain_size-length;
amrex::Real cylinder_center = domain_size-0.5*length;

#if(AMREX_SPACEDIM == 3)
EB2::SphereIF sphere(radius, {AMREX_D_DECL(x_location,x_location,sphere_center)}, false);
EB2::CylinderIF cylinder_z(radius, length, 2, {AMREX_D_DECL(x_location,x_location,cylinder_center)}, false);
auto cylinderphere = EB2::makeUnion(sphere, cylinder_z);
#else
EB2::SphereIF sphere(radius, {AMREX_D_DECL(x_location,sphere_center,0.0)}, false);
EB2::CylinderIF cylinder_y(radius, length, 1, {AMREX_D_DECL(x_location,cylinder_center,0.0)}, false);
auto cylinderphere = EB2::makeUnion(sphere, cylinder_y);
#endif
////////////////////////////////////////////////////////
// define EB
auto gshop = EB2::makeShop(cylinderphere);
EB2::Build(gshop, geom, required_coarsening_level, max_coarsening_level);
/////////////////////////////////////////////////////////////////////////////////////////

My solution for the 3D Cartesian case agrees well with the axisymmetric Afivo-streamer solution, and my solution for the 2D Cartesian case agrees well with the 2D Cartesian Afivo-streamer solution. The only incorrect result is obtained when I use 2D AMReX with cylindrical coordinates. I created the two examples that are attached above to check the behavior of the MLEBABecLap operator in all three cases.

I wanted my code to be applicable to 2D Cartesian, 2D axisymmetric and a 3D Cartesian case, just as the Afivo-streamer code. However, if EB is not supported with cylindrical coordinates, I will use only Cartesian coordinates whenever EB is included in my simulations.

[WeiqunZhang]

You might be able to use 2D Cartesian and then give MLEBABecLap variable coefficients with void setEBDirichlet (int amrlev, const MultiFab& phi, const MultiFab& beta) and void void setEBHomogDirichlet (int amrlev, const MultiFab& beta);. Those variable coefficients should be able to mimic the metric terms. You also need to modify the right-hand size.

This is because

$$ \frac{1}{r} \frac{\partial}{\partial r}(r \frac{\partial \phi}{\partial r}) + \frac{\partial^2}{\partial z^2} \phi= \rho $$

can be transformed into

$$ \frac{\partial}{\partial r}(r \frac{\partial \phi}{\partial r}) + \frac{\partial}{\partial z}(r \frac{\partial \phi}{\partial z}) = r \rho. $$

[git-user-1p155]

Thank you very much for your detailed reply! I will implement the solution which you suggested.