Possible typo in Projector/ProjectorAM2Order.cpp

[Status-Mirror]

Hi all,

I've been going through the source-code trying to understand the cylindrical-PIC algorithm, and there's a line in Projector/ProjectorAM2Order.cpp which looks wrong to me. On lines 140 and 834, we find:

 double r_bar = ((jpo + j_domain_begin_)*dr + deltaold[1*nparts] + rp) * 0.5; // r at t = t0 - dt/2

From the comment, I would guess this is trying to average the radial particle position before the particle push, and the radial position after. My understanding of these terms are as follows:

• (jpo + j_domain_begin_)*dr: Radial position of closest cell-edge to the particle position before the most recent push.
• deltaold[1*nparts]: Radial distance from this cell-edge to the old particle position, divided by dr.
• rp: Radial position of the particle after the push.

I think this line is attempting to work out $\bar{r}=(r_{old} + r_{new})/2$, but I don't believe this is doing it. Should the line instead read:

 double r_bar = ((jpo + j_domain_begin_ + deltaold[1*nparts])*dr + rp) * 0.5; // r at t = t0 - dt/2

I'm still working through the source-code and my understanding is incomplete. There could be a good reason why the line is as it is, but I think it's likely there's a typo/bug here.

Cheers,
Stuart

[Status-Mirror]

Hi developers,

I've been continuing my study of Projector/ProjectorAM2Order.cpp, and I think I've found a second typo - this time in the m=0 $J_\theta$ calculation.

On line 251, we find:

                Jt [linindex] += crt_p*(Sr1[j]*Sl1[i]*e_delta_inv - Sr0[j]*Sl0[i]*( e_delta-1. ));

In this bug report, we will take the correct form of $J_\theta$ from a single macro-particle to be:

$J_\theta = Q W \bar{S} v_\theta / V$

Where $Q$ is the charge of a real particle within the macro-particle, $W$ is the macro-particle weight, $\bar{S}$ is the average fraction of the macro-particle weight within the cell of interest, $v_\theta$ is the macro-particle speed in the $\theta$ direction, and $V$ is the cell volume. While the m=0 $\hat J_\theta$ should include a factor of $1/2\pi$ (see just after (42)), all currents in this function are calculated without this factor, so we will ignore it here for consistency with $J_x$ and $J_r$. The individual variables in the code for imode=0 refer to:

• crt_p: $Q W v_\theta / 2\pi \Delta x \Delta r$, where $\Delta x$, $\Delta r$ are the $x$ and $r$ cell sizes.
• Sr1[j]*Sl1[i]: Total weight fraction in cell $(i,j)$ after the push, divided by radial position of evaluation point ($S_{new}/r$)
• e_delta_inv: 0.5
• Sr0[j]*Sl0[i]: Total weight fraction in cell $(i,j)$ before push, divided by radial position of evaluation point ($S_{old}/r$)
• e_delta: 1.5

Thus, when we combine all terms, we end up with the expression:

$Q W \frac{1}{2}(S_{new} - S_{old}) v_\theta / V$

where cell volume $V = 2\pi r \Delta x \Delta r$. The SMILEI code would agree with my analytic expectation if we calculated $\bar{S} = \frac{1}{2} (S_{new} + S_{old})$, but we seem to have a minus instead.

If I am correct, this typo can be easily fixed by changing e_delta from 1.5 on line 186 to:

    e_delta = 0.5;

This should not affect the subsequent code, as this e_delta is only used for the $m=0$ mode of $J_\theta$, and is overwritten in line 198 for the higher modes.

Cheers,
Stuart

[Z10Frank]

Hello,
thank you very much!
We'll check those lines as soon as possible.

[beck-llr]

Hi Stuart

You are correct on both points.
I have tried a few tests with the corrections you propose and saw very little difference in the results. The most important difference I see is the noise of Jt mode 0 being significantly higher. It was probably underestimated before.

Nevertheless, these corrections should be introduced. I'm just going to run a couple more tests and you should see them appear soon in the development branch.

I'd like to thank you very much for the very careful reading of the code and for reporting your findings to us in a crystal clear post. Congratulations to you and I really look forward to reading your next suggestions to improve the code quality :-)
Also if you have questions about some part of the code, do not hesitate to contact us on the chat.

Cheers

P.S: I'll close the issue when the corrections are available.

[Status-Mirror]

Hi @beck-llr

Thanks for the reply, it's good to know that my understanding of the code is correct! I have now finished my review through the code, and I think I understand all the equations which are relevant to the core cylindrical-PIC algorithm. During this review I found a few more points which could also be potential bugs and typos, and I'm happy to share my findings with you.

Initially I tried typing this out here, but I hit some kind of limit on the GitHub markdown service. The response became so long that equations stopped formatting correctly! As a workaround, I'll provide a link to the report on Overleaf, as I can't seem to attach the PDF here.

Smilei bug report: https://www.overleaf.com/read/hptynqbzzrpj

Thanks for all the help,
Stuart

[beck-llr]

Hi @Status-Mirror
Thanks for this massive report. I'll try to address the different points progressively as I find time to read through them all.
About point 2, I think you are getting something wrong.

The mode 0 of a given current density $J$ that we note $\tilde{J}^0$ is defined by $$\tilde{J}^0(x,r)=\frac{1}{2\pi}\int_0^{2\pi}J(x,r,\theta)d\theta.$$
I think we agree on that definition.
Now if $J$ is homogeneous (independent of x, r and $\theta$) and its value is $-en_ec$, we expect
$$\tilde{J}^0(x,r)=\frac{1}{2\pi}\int_0^{2\pi}-en_ecd\theta=-en_ec.$$

Note that the $\frac{1}{2\pi}$ factor is actually not expected. Please let me know if that answers your concern. Feel free to update your overleaf document accordingly. This document is very convenient.

Cheers

[beck-llr]

Let's move to point 3, boundary condition of $\tilde{B_r}^m$ on axis for $m> 1$.

First of all, note that the $\tilde{B}^m$ fields are functions of $x$ and $r$ and not of $\theta$ so your equation 1) does not make sense even though I think I get what you mean.

Then, it is important to understand that the so called "below axis" (see https://smileipic.github.io/Smilei/Understand/azimuthal_modes_decomposition.html) boundary conditions are not determined by physics but are simple practical recipe so that conditions on axis actually work.

Let's take your example of $\tilde{B_r}^m$ for $m> 1$. Here I will note $B_r^m$ (without ~) the part of $B_r$ that is reconstructed from $\tilde{B_r}^m$ at a macro-particle position by the interpolation.

Let's assume that this macro particle is sitting on axis (or very close to it) at position $(r,\theta)$ with $r << dr$. Then the interpolator gives

$$ B_r^m = \left(0.75\tilde{B_r}^m[2] + 0.125\tilde{B_r}^m[3] + 0.125\tilde{B_r}^m[1] \right)\exp{im\theta}.$$

You agreed that $\tilde{B_r}^m[2]$ must be zero. If you want $B_r^m$ to be zero then you need to have $\tilde{B_r}^m[1]=-\tilde{B_r}^m[3]$ for all $m>1$.

I hope this is clear enough. Again let me know if you see any objections to that argument.
Cheers

[Status-Mirror]

On point 2, I agree - I've discovered a contradiction in my argument. The macro-particles have no "shape" in the $\theta$ direction, as the interpolator just uses $\exp(im\theta)$ without averaging over nearby $\theta$. Hence, when calculating the current density contribution from a single macro-particle, I took $J(x,r,\theta) = J(x,r)\delta(\theta-\theta_p)$ where $\theta_p$ is the particle $\theta$. Solving this integral did not introduce a cancelling $2\pi$ factor. The error in my approach was that if I take a macro-particle shape with a missing dimension, then the macroparticle number density increases by a factor of $2\pi$, which provides the cancellation factor.

I see your approach to point 3 is provided in your documentation, but I would argue there is an inconsistency between the treatment of the fields and the current densities (see my bug 6 discussion). Your conditions are correct for the goal: "If you want $\hat{B}_r^m$ to be zero", but I don't think you should force that. If the macro-particle shape can deposit current below $r=0$ in a physical way, then why can't the fields below $r=0$ influence acceleration on the macro-particle?

Following current density logic, I think the interpolator should give:

$$B_r^m = \left( 0.75 \hat{B}_r^m[2] + 0.125 \hat{B}_r^m[3] + 0.125 \hat{B}_r^m[1]\right) \exp(im\theta)$$

where, from the geometry of the physical space, $\hat{B}_r^m[1] \exp(im\theta) = -\hat{B}_r^m[3] \exp(im(\theta+\pi))$. This would suggest:

$$B_r^m = \left( 0.75 \hat{B}_r^m[2] + 0.125 \left( \hat{B}_r^m[3] - (-1)^m \hat{B}_r^m[3]\right)\right) \exp(im\theta)$$

Using this interpolation method would make everything consistent, but maybe there's a reason you've chosen not to use it? I see that for a particle with exactly $r=0$, then $\theta$ becomes undefined. Perhaps you've chosen the non-physical method to prevent errors in evaluating $\exp(im\theta)$ when $r=0$?

Ultimately I think point 3 is an implementation choice, and I don't think the differences in our approach would radically change the output of the simulation.

Cheers,
Stuart

[beck-llr]

For the point 4.1 you are correct, the minus sign must be replaced by a plus. Congrats for spotting this one out !

[mccoys]

@beck-llr anything left to do in this issue?

[beck-llr]

Yes I'm not done yet. Unfortunately this takes time and I haven't had time to prepare a proper answer to the interesting points raised here.

[beck-llr]

@Status-Mirror About point 4.2: missing currents in the Buneman boundary conditions.

First please note that the boundary conditions on particles are applied before current deposition. Therefore a particle crossing r=rmax would be removed or reflected before it has a chance to deposit currents as far as you suggest.

Moreover the Buneman boundary conditions are not derived with the objective to solve locally Maxwell's equations but to analytically try to force minimum reflections of outgoing waves back into the domain. In order to do so you have to assume some conditions on the state of the plasma "outside" the simulation domain and that assumption is that there is no current. Of course if there is a strong current at the boundary this assumption becomes false and this is a limitation of this method.

[beck-llr]

@Status-Mirror About point 4.3: you are correct the -1 here is kind of arbitrary and not really justified.
The question of boundary conditions in the corners is a very difficult one since several boundary conditions could possibly apply on top of each other and therefore overwriting each other. One has to choose which boundary is the most relevant. It has been decided that longitudinal boundaries are. And that is why the loop for $B_t$ in the very same function starts at 1 and finished at $n_p$ (which is at 1 cell before the end since $B_t$ is dual along x). So that the extreme points of $B_t$ are defined by the longitudinal boundaries and not the radial ones.
A similar logic does not apply for $B_l$ since it is not affected by longitudinal boundaries hence the loop starts at 0 and should finish at $n_p$ as you noticed.

P.S: In the case of the spectral solver the fields do not have exactly the same size and require that this -1 is used. So instead of having a different implementations of the Buneman conditions for Rmax, we decided to keep the spectral version, including this -1, since for the FDTD solver we favor using the PML boundary conditions instead. Agreed this is a lazy implementation. Maybe we'll try to improve this sometime.

[beck-llr]

@Status-Mirror About 5.1 you are absolutely correct. I've run some basic tests and this correction has some impact on the Bl component of the reflected laser on the xmax boundary. Again, congratulation for spotting it and thanks a lot for reporting it.

5.2 is obviously an ugly typo :-)

[charlesprouveur]

@beck-llr is there anything left to do in this issue or we can close it?

[beck-llr]

Somme additional comments have been addressed but I haven't had time to answer here yet. Let's keep it open.