magnetized shock

[LudwigBoess]

This PR adds the option to run a magnetized shock with a given angle between magnetic field and shock propagation.
The magnetic and electric field is set up in such a way that it exerts no additional Lorentz force on the particles.
I implemented this into the normal shock pgen to avoid duplicate code.
The default behavior can be recovered by setting Bmag = 0.0 in the parameter file. For a magnetized shock Bmag = 1.0 should be set to use the internal magnetic field units.
Please let me know if this option is fine for you, or if I should set up an independent pgen after all.

@jmahlmann and I found that there's an underlying issue with the shock setup where the field boundaries are not treated correctly. This causes an additional wave to propagate from the left and right side of the domain, which distorts the particles.

Please let me know if this problem should also be addressed in this PR, or in a subsequent PR, as it also affects the default shock setup.

[LudwigBoess]

To clarify the issue with the boundary conditions:

In the default setup I can see two waves propagating through the phase diagram from both ends of the shock setup and I don't think this should happen:

In the last commit I added the option to set up an atmosphere boundary condition.

If I set up this boundary condition to reset the magnetic field to the initial configuration, I can get rid of the rightward propagating wave:

However, I'm not sure if that's really the correct solution, but I'm going to hide behind my naivety towards PIC simulations for now...

[haykh]

@LudwigBoess i rebased this to the latest version + changed the Field boundary condition name to FIXED which should be a bit more generic (basically, it drives to a fixed value). one question remains, though. in FIXED (formerly ATMOSPHERE) field BCs we have a buffer region where the fields are set. this is important, when FIXED field BCs is used together with ATMOSPHERE particle BCs (which is why i had them like that in the first place). but now, FIXED fields can be used with, say, REFLECT. the question is -- should the buffer still be there?

[LudwigBoess]

Thanks a lot, @haykh!
I did a quick check with the phase space plot I made last time and the weird waves are gone now!

So at least phenomenologically I'd say the buffer should be fine.
(edit to clarify: this is a 1:100 mass ratio shock with drift_ux = 0.2 and temperature = 1e-5, so slightly different setup from last time.)

[haykh]

after giving it a thought, i actually came up with a slightly better convention. will push shortly, but basically i suggest we keep the atmosphere as is (with buffers), and simply add "FIXED" field boundaries. the difference will be that FIXED boundaries only modify values on the surface (no buffer). also, since in many cases FIXED boundaries may be just constant, I added an option to simply specify a function in the pgen that returns a single value per each component, which is then simply copied (without invoking a redundant kernel).

i'll add all this to the wiki and will push it forward for you to check.

[LudwigBoess]

That sounds like a good option! I'll wait for your push and then re-run my tests to see if the magnetic field setup works the way it should.

[LudwigBoess]

On second thought: would it be possible to add a boundary condition INIT that uses InitFields in the same way that ATMOSPHERE used DriveFields?
That would avoid duplicate code and would allow me to play around with DriveFields independently.
I saw in Tristan they use the drift-field I set up as init field on the right of the shock and then set the ey and ez component of the reflecting wall to 0.
With those independent boundary conditions, I could implement this behavior as well.

[haykh]

i'm a bit confused. can you explain again the exact boundary conditions imposed?

currently, in X-direction, it's using ["FIXED", "ABSORB"] for fields (and ["REFLECT", "ABSORB"] for particles).

• "FIXED" -- set constant fields (initial values) exactly at the boundary
• "ABSORB" -- damp fields to zero over a certain distance (we can change this to damp to any predefined value)

[haykh]

and this is what the function which sets the fields for "FIXED" looks like:

    auto FixField(const em& comp) const -> real_t {
      if (comp == em::ex2) {
        return init_flds.ex2({ ZERO });
      } else if (comp == em::ex3) {
        return init_flds.ex3({ ZERO });
      } else if (comp == em::bx1) {
        return init_flds.bx1({ ZERO });
      } else {
        raise::Error("Other components should not be requested when BC is in X",
                     HERE);
        return ZERO;
      }
    }

it simply returns a number, which it takes from init_flds.

[LudwigBoess]

The fixed BCs look good thanks a lot for implementing this!
There was a bug in the Bfield setup, so please disregard the last plot for my previous test.

I fixed that now and switched the BCs to be consistent with Tristan. This fixes the issue with the wave propagating from the right to the left.
However, there is still the wave propagating from the left to the right, causing the particles to be stopped.
You can see this in the phase diagram here:

whose middle panel corresponds to the second-to-last panel here:

Tristan seems to set the field components to zero instantly on the left side of the tube, instead of damping it to zero as you mentioned it's done in entity. Is there a parameter I can tweak to set the fields to zero instantly or is that hard-coded?

[haykh]

@LudwigBoess can you outline exactly what boundaries are being imposed in tristan (both on E and B field, on what scale, which components etc)? for the scale across which you set things to zero: there is an absorbing width parameter in the input file (in physical units, not cells). alternatively, i can also add a possibility to set two different "FIXED" boundaries, which will only enforce a single cell.

i think ultimately i will add a more general "MATCHING" boundaries (instead of "ABSORB"), which damp fields to specified values across distance, and then the "FIXED" ones will enforce only in a single cell.

[LudwigBoess]

@vanthieg and I just sat down and tested the new BC.

There is still the propagating wave and especially the dipole forming for the theta=90° setup is weird (output is at t=2.0):

To further test where this is coming from we set the charges of the particles to zero to see if it's an issue with the field BCs and we found this strange behavior for the Ez and By components (plotted are t=0-10 with dt=1):

Do you have an idea what's going on here?

You can find the parameter file for the 90° test here: https://1drv.ms/u/s!AnV4mghSO5Lzgc01Fh_TXmg0hHD7BQ?e=POboCu
All other tests are identical, with only Btheta set to the required angle.