Add external particle fields ohms law hybrid

[clarkse]

This PR allows for the addition of external fields through the particle fields analytical interface. This is useful for field splitting external vs. self fields in the hybrid ohm's law solver.

[roelof-groenewald]

I've made some progress but have not fully reviewed all the changes yet. Here are some comments I've recorded so far.

[clarkse]

@roelof-groenewald thanks for the partial review. I will start making some of the non-trivial changes and rechecking the tests. Very thorough review as usual. This will make for a better PR.

[roelof-groenewald]

Overall it looks good! I just had a final few small comments and a somewhat bigger question:
In the solve loop $E_{ext} = -\partial A/\partial t$ is subtracted from the E-field (so in total $E = E_{plasma} - E_{ext}$ which is then used to update the B-field, giving $$B_{new} = B_{old} -\nabla \times E_{plasma} dt + \nabla \times E_{ext} dt.$$
But $\nabla \times E_{ext} dt = -\partial/\partial t \nabla \times A dt = -\partial/\partial t B_{ext} dt$, so in effect
$$B_{new} = B_{old} -\nabla \times E_{plasma} dt - B_{ext}^{t+1}.$$
Then in the next half of the field push, $B_{ext}^{t+1}$ is added to $B_{new}$, cancelling the effect of subtracting $E_{ext}$ from the E-field in the first place.
And then at the very end of the field-solve $B_{ext}$ is added to the B-field again. It seems to me it should be sufficient to have $E_{ext}$ added to the E-field at every evaluation of E, and that it would then carry the change in $B_{ext}$ through to Faraday's law (the B-field update). Maybe we can chat about this quickly tomorrow.

Here's what I mean in slightly easier to read font:

[roelof-groenewald]

Why is the external field subtracted rather than added?

[clarkse]

So it is set up this way since the external fields are what they are, they operate in the soak through limit. So because of this since it is not solving A via a Poisson solve with totally self-consistent boundaries. I initially tried to just use the numerical integration to integrate the total B field, but because if the user is not extra careful to make sure the external field is self-consistent with the boundaries this approach was taken. I am doing a subtraction at the beginning of the Hypric PIC loop to remove the external component from the total field solution, then the E and B fields are the plasma response fields only moving through the integration loop. This allows for the Hybrid loop to proceed with integrating the plasma response into B
$\vec{E_{plasma}} + \vec{E_{ext}} = \frac{(\vec{J_{plasma}}-\vec{J_i}) \times (\vec{B_{plasma} + \vec{B_{ext}}})}{ne} - \frac{\nabla p_e}{ne} + \eta \vec{J_{plasma}} - \eta_h \nabla^2 \vec{J_{plasma}} $

Since the E and B fields at this point are the plasma fields in the Ohm's Law solver, however the solution is for a total field in order to get the plasma response only the E_ext needs to be subtracted from both sides of the equation.

In the vacuum region it gets a bit hairier. If the resistivity in that region is high enough or if the holmstrom flag is used then the Generalized Ohms Law asymptotes to:
$\vec{E_{resistive}} \sim \eta \vec{J} $
This can be linearly added to the vacuum field to get the total field since the assumption is that it is a vacuum and not dense enough to shield an external field.
This resistive electric field gets integrated into the B field, then the resistivity term is dropped because of the cancellation of the drag term in the ion momentum update for the particle push.
Since the plasma response in the vacuum region is not well defined the algorithm switches to using this form of Ohm's law to have a smoothly varying E field around the density cutoff to limit vacuum fluctuations integrating into the plasma field. Since in the vacuum field the plasma E field should be zero except for some damping electric fields that show up near the edges of the plasma, then this gets integrated into B with the boundary conditions. E_ext is not subtracted off in this case as the plasma shielding response is just not accurate in the vacuum region. This really cuts down on the vacuum fluctuations for me and makes the simulations more stable.

[roelof-groenewald]

Sorry to get stuck on this more, but looking at it again a few hours after we talked I'm again unsure about the following:
You say the E-field from Ohm's law is the total field (i.e. $E_{plasma} + E_{ext}$), is there a clear argument for why that is the case? I agree Ohm's law of course includes induction (from $v_e \times B$, where $B$ here includes the external field), but it seems it would still only be the E-field due to the plasma response. Shouldn't there be an additional inductive E-field from $-\partial A/\partial t$?
In such a case we'd have:

$$E_{total} = E_{Ohm} + E_{ext} = E_{Ohm} - \partial A/\partial t$$

such that the full Faraday equation would give:

$$\frac{\partial}{\partial t} B_{total} = -\nabla\times E_{total}$$

$$\frac{\partial}{\partial t} (B_{plasma} + B_{ext}) = -\nabla\times E_{Ohm} + \frac{\partial}{\partial t} \nabla \times A$$

giving

$\frac{\partial}{\partial t} B_{plasma} = -\nabla\times E_{Ohm}$ and $\frac{\partial}{\partial t} B_{ext} = \frac{\partial}{\partial t} B_{ext}$.

If the Ohm's law calculated E field includes $-\partial A/\partial t$, there would not be a way to have Faraday's law operate on the full E-field that would result in a B-field with both the plasma and external B-field included. I must still be missing something...

[roelof-groenewald]

In terms of implementation, it seems to me that if the Ohm's law E-field does not already include $E_{ext}$, the formalism here would be that a time varying external field can be mathematically included by replacing the E-field calculation with $E = E_{Ohm} + E_{ext}$ where $E_{ext} = -\partial A / \partial t$. In that case the B-field integration from Faraday would itself add the external B-field. (As you mentioned earlier, the numerics could be complicated by boundary handling, in which case treating $B_{plasma}$ and $B_{ext}$ separately is understandable).

But the point of difference with the current implementation is, for clarity:
Is $E_{total} = E_{Ohm} = E_{plasma} + E_{ext}$ or $E_{total} = E_{Ohm} + E_{ext} = E_{plasma} + E_{ext} $?

[clarkse]

$E_{ohm} = E_{total} = E_{plasma} + E_{external}$. I agree that it would be nicer if the integration worked out, but the boundary handling of the external field is finicky and not consistent with having a PEC at the EB or the domain edge. If we did a vector potential formulation of the integration I could do this all fully self-consistent. I would also like to be able to get rid of one of the sets of MFs. I cannot fix the boundary properly with having a soak through external field that ignores the boundaries and the high-frequency PEC response in this formulation.

I am exploring options to make this better, but I think it would require either having a displacement solver that has slowed down EM waves, or the vector potential formulation with infinite speed propagation. This is more of a solution to quickly add fields to try out. I think this is better suited to deal with in a follow up PR. I think I need to do a bit of work on the Magnetostatic solver, then I can address doing a self-consistent magneto-inductive version of the hybrid algorithm.

[roelof-groenewald]

Why do you say $E_{ohm} = E_{plasma} + E_{external}$? That is just the part I'm stuck on. It seems to me it should be the case that $E_{ohm} = E_{plasma}$.

I think the implementation of the model as you propose is good, so I'll approve this PR and merge now. But I think we should continue a discussion of the accuracy of the underlying physics model (and adjust the implementation in future PRs if we agree it is needed).

[roelof-groenewald]

It's curious to me what the clang-tidy tests focus on. I would think it would prefer this to be

if (holmstrom_vacuum_region && (rho_val < rho_floor))

due to the possibility to short-circuit and avoid the comparison.

[roelof-groenewald]

The code changes look good! Thanks for the all the work to get this implemented!