Add shock physics

[lemmatum]

In particular, Rankine-Hugoniot equations and other forms of compression (e.g. adiabatic). These are important for many plasmas, but also has overlaps with the dynamic compression research field (squeeze things to high pressures and find new phases of matter).

[namurphy]

I like this idea! I had to cover shocks in a class a couple of years ago, and at the time I was wondering how possible it would be to create an MHD shock solver. At the time I wasn't sure how to go about doing it because the MHD jump conditions are kind of horrible, but one possibility would be to create a Discontinuity class (or metaclass?) and use a web of setters and getters to calculate variables from the available attributes and conservation relations, or say that not enough information is available and be careful that we don't end up in an infinite recursion. Here's some code to illustrate what I'm thinking, with hypothetical functions that can get n from B and vice versa.

class Sample:

    def __init__(self, n=None, B=None):
        self._n = n
        self._B = B

    @property
    def B(self):
        return self._B if self._B is not None else calculate_B_from_n(self._n)

    @B.setter
    def B(self, B_):
        self._B = B_

    @property
    def n(self):
        return self._n if self._n is not None else calculate_n_from_B(self._B)

    @n.setter
    def n(self, n_):
        self._n = n_

There could perhaps be Shock with subclasses IdealHydroShock and IdealMHDShock, and then Discontinuity could also have subclasses TangentialDiscontinuity and RotationalDiscontinuity as well.

Let's brainstorm more about this later, since I want to play video games from the late 1980s before going to bed! 📺 🎮 💤

[lemmatum]

A discontinuity class could be a nice way of handling shock physics as most hydro simulations use an artificial viscosity or something similar to thicken and smooth out the shock such that it can be resolved on the simulation grid scale. This could then be further refined as an upstream and a downstream boundary which hands off to another, more finely resolved, simulation. I imagine we'll have to dig into the literature on these ideas.

For now, I just want a simple 1D function to give me an upstream or downstream shock parameter when i give it enough of the other conditions. For Rankine-Hugoniot I think it is usually giving it 5 parameters and it returns the 6th?

[namurphy]

I'm going to tentatively assign this to myself for v0.2.0 since I'm very interested in doing this, but I will not be able to get to it until probably later this summer. If anyone else is interested and there isn't a pull request that was already created, please go ahead!

[pheuer]

@PlasmaPy/peer-reviewers-shocks @rocco8773 @KhalilBryant

I wanted to resurrect discussion of this issue because I think the sort of RH-solver that @namurphy is proposing here would be a really useful component of the shocks package. Ideally I think such a function/object would take/be initialized with any combination of upstream and downstream conditions sufficient to constrain a solution, then return the full upstream and downstream states as defined by the RH conditions.

I tried using fsolve to numerically solve the RH equations, but that often failed at large or small values of the plasma parameters: perhaps writing the equations in dimensionless terms would improve that. Symbolic math might be another possible approach. Does anyone else have any ideas on how to solve this problem?

Simply hard-coding the one most common solution (find downstream values in terms of upstream values) would be trivial, but not nearly as broadly useful.

[StanczakDominik]

I'm a complete newbie in this particular subject, but: how would you envision symbolic math working here? I've been tinkering with sympy recently and could definitely help implement something, if we can design something to implement.

[pheuer]

@StanczakDominik I have limited experience with sympy, but the task would be to somehow take the RH equations, solve them symbolically for the desired output variables in terms of the input variables, then "lambdify" and evaluate those expressions to get the output variables. That doesn't sound like a very robust process to me, but maybe it's possible?

[StanczakDominik]

If they have closed form solutions, I honestly don't see why not. Do they?

[pheuer]

I think they should? The set of equations to be solved for the full MHD problem is 933-938 here: http://farside.ph.utexas.edu/teaching/plasma/Plasmahtml/node79.html

although a simpler case to start with for development purposes would be the B=0 hydrodynamic case which is equivalent to the parallel MHD case here (eq 942-944):
http://farside.ph.utexas.edu/teaching/plasma/Plasmahtml/node80.html

For the purposes of the solver, gamma could be considered constant, but we would want to solve for an arbitrary (but sufficiently determined) subset of the densities, velocities, pressures, and magnetic fields.

[StanczakDominik]

That looks workable! I suppose we could try to solve the full MHD problem, then use the B=0 as a test case. Could you expand on this part? I'm not 100% sure I get what you mean here:

an arbitrary (but sufficiently determined) subset of the densities, velocities, pressures, and magnetic fields.

[pheuer]

Sure, I just mean that there are six equations and twelve variables (Vx, Vy, Bx, By, P, rho in both the upstream and downstream), so the user must specify the values of exactly six of those variables in order to adequately constrain the problem. It should be able to be any six: users may know the entire upstream and want to know the downstream conditions, or they might know a mix of the two.

[StanczakDominik]

Oh that is entirely workable, then. I should have a draft PR (probably in notebook form for now) soon - the optimistic case is tonight, even.

[pheuer]

That's awesome, looking forward to it!