Modify solver

[Aratopn]

Hello,

I am trying to simulate T.Charoy et al benchmark (DOI 10.1088/1361-6595/ab46c5).
Among other things, it involves recalculating the potential after solving the Poisson equation, taking into account the emission plane.

The mentioned benchmark uses the following approach:

1. First, the Poisson equation is solved with the boundary conditions $U(0,y)=200V$ and $U(Lx,y)=0V$.
2. Next, the potential on the emission plane $\overline{U_e}$, averaged in the $y$ direction, is determined.
3. Then, the potential distribution is recalculated according to the expression: $\phi(x,y)'=\phi(x,y)-\frac{x}{x_e}\overline{U_e}$.

Here, $x_e$ is the emission plane position.

I implemented this using a callback function:

def updatePhi():
    phi = sim.fields.get("phi_fp", level=0)
    phi_old = phi[...]
    phi_emi_mean = np.mean(phi[i_emi])
    phi_new = phi_old - (X_grid / x_emi) * phi_emi_mean
    phi[...] = phi_new
    phi_x_hi = np.mean(phi_new[-1])
    sim.extension.warpx.set_potential_on_domain_boundary(potential_hi_x=f"{phi_x_hi:.6f}")
    
"""
some code
"""

callbacks.installbeforeInitEsolve(updatePhi)
callbacks.installbeforeEsolve(updatePhi)

X_grid is precalculated numpy meshgrid with x positions for all nodes.
x_emi is the emission plane position.
i_emi is the index corresponding to x_emi.

I also added a function that should return back the potentials on the boundaries before a new step, so that the Poisson equation in this step is solved with the boundary conditions: $U(0,y)=200V$ and $U(Lx,y)=0V$.

def setZeroPotential():
    sim.extension.warpx.set_potential_on_domain_boundary(potential_hi_x="0.0")
    
"""
some code
"""

callbacks.installbeforestep(setZeroPotential)

However, this approach doesn't seem to work

When I use the function updatePhi in the next step, the calculated potential at the callback location "beforeEsolve" seems to use the updated potentials at the boundaries, not the original ones: $U(0,y)=200V$ and $U(Lx,y)=0V$.

Briefly about what I'm trying to achieve:
==step start==

1. set $U(0,y)=200V$ and $U(Lx,y)=0V$
2. solve Poisson Equation with $U(0,y)=200V$ and $U(Lx,y)=0V$
3. update phi: $\phi(x,y)'=\phi(x,y)-\frac{x}{x_e}\overline{U_e}$
4. calculate electric field and move particles with new phi
5. do the rest

==step end==

If you have a moment, I would be grateful for any insight you could share on how this can be achieved. Any details—such as the general approach, specific functions used, or a pointer to a relevant code section if it is publicly available—would be immensely helpful for my work

Thank you for your time.

[roelof-groenewald]

Hi @Aratopn. Thank you for your question. What you are trying to do can definitely be accomplished in WarpX and you are on the right track.
The main problem with your current approach is that your callbacks are installed before the Poisson solve occurs, and during the solve (after your callbacks) phi is overwritten with the new values. You should use:

callbacks.installbeforeEsolve(setZeroPotential)

and

callbacks.installafterEsolve(updatePhi)

One complication is that the E-field is calculated from phi within the Esolve step (this is because with EBs the E-field accounts for cut-cells), so you will need to update the Efield_fp multifab after you update phi. This just requires a simply finite differencing pass though. Let me know if you need guidance on how to do that.

[Aratopn]

Thank you very much for the detailed and helpful answer.

I thought that the potential is calculated at the callback location "poissonsolver", and at the callback location "beforeEsolve", the calculation of the electric field begins, taking the new potential into account, and ends at the callback location "afterEsolve".

I do need some guidance on the finite-differencing part to update Efield_fp. Could you please provide a small example or hint on how to approach it?
When I use: phi = sim.fields.get("phi_fp", level=0) phi[...] has shape (501, 257). I set number_of_cells=[500, 256], so as far as I understand these are the values ​​on the nodes.
But then I use Ex = sim.fields.get("Efield_fp", dir="x", level=0) for some reason Ex[...] has shape (500,257,1).

I also had the idea that if I use MPI sometimes, then guard cells should appear and I should somehow overwrite them too?

Thank you for your time.

[roelof-groenewald]

I can definitely understand the confusion there. The poissonsolver callback actually replaces the Poisson solver. It is used in cases where an external Poisson solver is used (see for example https://github.com/BLAST-WarpX/warpx/blob/development/Examples/Physics_applications/capacitive_discharge/inputs_test_2d_background_mcc_picmi.py).

The difference in sizes between phi and Ex is due to the staggering of the Yee mesh (when using the energy-conserving PIC scheme) where phi is defined on the grid nodes (as you correctly said) but Ex is on the edges along the x-direction. In order to get E from phi you can simply do:

        phi_full = phi[()]
        Ex[()] = -(phi_full[1:] - phi_full[:-1]) / dx
        Ey[()] = -(phi_full[:, 1:] - phi_full[:, :-1]) / dy
        Ez[()] = -(phi_full[:, :, 1:] - phi_full[:, :, :-1]) / dz

where dx, dy, and dz are the grid spacing values. The line phi_full = phi[()] grabs the values of phi including the guard cells, and similarly writing to Ex[()] populates the regular and guard cells for Ex.

Please give this a try and let us know if you run into any problems.