Initial RHSolver attempt

[StanczakDominik]

Tackles #440. For now I'm developing this in a notebook.

I've got something that should theoretically work properly, but seems to be running into computational issues - the notebook is still running on my machine! I'm putting this up as an example of where I'm at with the initial idea.

• I have added a changelog entry for this pull request.

• If adding new functionality, I have added tests and
  docstrings.

• I have fixed any newly failing tests.

[review-notebook-app]

Check out this pull request on 

Review Jupyter notebook visual diffs & provide feedback on notebooks.

---

Powered by ReviewNB

[codecov]

Codecov Report

Merging #887 into master will not change coverage.
The diff coverage is n/a.

@@           Coverage Diff           @@
##           master     #887   +/-   ##
=======================================
  Coverage   96.18%   96.18%           
=======================================
  Files          59       59           
  Lines        5367     5367           
=======================================
  Hits         5162     5162           
  Misses        205      205           

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a27ef92...36e2257. Read the comment docs.

[StanczakDominik]

Okay, I've got something that takes a few seconds but works decently. What helped was substituting (fake, for now) values for Gamma and mu_0; I'll be adjusting this to work with physical Quantity inputs, which should help things.

@pheuer could you take a look? This is of course at a very brief stage and I'll be redesigning most of it, but is the general idea sensible?

Also, I'm somehow blanking on what Gamma was supposed to mean here 😅 Still, I take it I can just take it as (the only required?) input, substitute its numerical value into the equations, then happily solve them.

[StanczakDominik]

The interface is taking shape, and that's reviewable; I'm still fighting issues with Sympy's solution and execution speed. It doesn't seem to be handling inputs instead of my initial "set everything to 1!" attemps particularly well. Perhaps an issue is that I'm substituting numerical values in before the solution, to limit the number of variables it has to work with from 14 (6 inputs, 6 outputs + gamma + mu0) to 6 + numbers - which it seems to be representing sometimes as rationals and sometimes as its own Float type.

[pheuer]

@StanczakDominik This looks like a good start! Here are some initial comments:

• Gamma is the adiabatic index or heat capacity ratio /gamma: it could have a default keyword value of 5/3 for a monatomic gas in 3D. I don't know why this author uses uppercase rather than lowercase gamma... it's confusing.

• The evaluation speed is certainly not ideal, but hopefully that can be optimized. If execution time does remain > 1 s, maybe it would ultimately be a good idea to encapsulate this into an object so that the equations were only solved once (at initialization) but could then be re-evaluated for different values of the six chosen variables. That would speed up making a plot where you want to solve the same equations for many different input values.

• I can see how reducing the problem to six symbolic variables is probably necessary for computational efficiency, but I anticipate problems here based on the scale of the numbers involved. In space, say, typical densities, velocities, and magnetic fields would be 10^6 m^-3, 5 x 10^5 m/s, and 1 x 10^-9 T respectively. Left unchecked, I think numerical issues are very likely.

• The existence of multiple solutions is interesting: I suspect one must be non-physical but I'm not sure?

• I seem to get errors sometimes when specifying different sets of variables: another good test case would be specifying rho, Bx, and Vx upstream and downstream and solving for the other components. In particular, I got this error a lot:

TypeError: unsupported operand type(s) for *=: 'Symbol' and 'CompositeUnit'

I will keep thinking about this too.

[StanczakDominik]

Well... I think we may have a runtime problem here. In trying to solve the full 14 variable set of equations (including gamma and mu_0) for the 6 output ones, I'm getting runtimes of over 1 hour - and it's not like I even get the results, I just stop the computation after it's no longer reasonable. It's bonkers. I didn't think this was such a hard problem.

I still have a few ideas for how to do this. There's a 2019 report on Stack that moving the solution out of the jupyter notebook may help a little, but I have mixed feelings on whether that's actually applicable to our use case.

[StanczakDominik]

Nope, I think we can scrap that idea. The runtimes are atrocious. I've been tinkering with different flags to solve through the day; none of them completed.

[pheuer]

@StanczakDominik Well I suppose that's not too surprising: 14 variables is a lot!

It seemed like the main problem with your previous solution (inserting the known variables) was getting sympy to represent them all as double floats so that very small or large coefficients would retain their precision: maybe there is a way to force sympy to always use float representations?

Another approach might be to go back to the numerical technique I tried initially (using fsolve), but I think that would only work if we could re-write the equations in terms of some dimensionless quantities that would not be particularly large or small. Then again, if we could do that, maybe the sympy constants wouldn't be a problem either?