Added/modified 1D and 3D Maxwellian speed/velocity distributions #152
Conversation
plasmapy/physics/distribution.py
Outdated
|
|
||
| Returns | ||
| ------- | ||
| f : Quantity | ||
| probability in Velocity^-1, normized so that: | ||
| probability in Velocity^-1, normalized so that: |
plasmapy/physics/distribution.py
Outdated
| @@ -60,6 +65,8 @@ def Maxwellian_1D(v: units.m/units.s, | |||
|
|
|||
| .. math:: | |||
| f = \sqrt{\frac{m}{2 \pi k_B T}} \exp(-\frac{m}{2 k_B T} (v-V)^2) | |||
| f = (\pi * v_Th^2)^{-1/2} \exp(-(v - v_{drift})^2 / v_Th^2) | |||
There was a problem hiding this comment.
Could we add something like "or equivalently," or "two equivalent forms:" here?
That might break with the .. math:: .rst tag. We'll have to keep an eye on that.
There was a problem hiding this comment.
Oh, I was intending to remove the first form listed there. I think it makes more sense to look at the function as a comparison of velocities to obtain a probability with respect to the most likely velocity rather than having temperatures, and masses floating around. We can always add the thermal velocity function under "see also"
There was a problem hiding this comment.
It probably wouldn't hurt too much to have both of those available, but it's not critical to have both here, I suppose!
plasmapy/physics/distribution.py
Outdated
| # accounting for thermal velocity in 3D | ||
| vTh3D = np.sqrt(3) * vTh | ||
| vThSq = vTh3D ** 2 | ||
| # Get particle velocity squared |
There was a problem hiding this comment.
This isn't a big deal, but, "relative" particle velocity, maybe, given it's shifted by the drift one?
plasmapy/physics/distribution.py
Outdated
|
|
||
| Example | ||
| ------- | ||
|
|
There was a problem hiding this comment.
Could definitely use some doctests here :)
There was a problem hiding this comment.
How do you guys run doctests?
There was a problem hiding this comment.
That documentation update really can't come soon enough, can it? :D I'll try to automate and simplify the process a bit while I'm at it.
If you clone the repo to a directory called PlasmaPy (so you have a PlasmaPy/plasmapy directory where the actual repository is located, you should be able to run pytest or more specifically pytest plasmapy/physics/distribution.py from the PlasmaPy directory and that should, due to seeing the setup.cgg file, include doctests.
plasmapy/physics/distribution.py
Outdated
|
|
||
| Example | ||
| ------- | ||
|
|
There was a problem hiding this comment.
As previously, a doctest would be nice! We could go with a single particle for simplicity.
| integ = spint.trapz(y=yData1D, x=xData1D) | ||
| print("1D integral {0}".format(integ)) | ||
| exceptStr = "Integral of distribution function should be 1." | ||
| assert np.isclose(integ.value, 1), exceptStr |
There was a problem hiding this comment.
How about this, using the power of python 3.6 (though a .format would of course work as well, except for being less awesome):
# cut the print
exceptStr = f"Integral of distribution function should be 1 and not {integ}."
assert # as previously except it returns the value now
There was a problem hiding this comment.
oops, I'll delete the print line, that was mainly for my own testing. I like the idea of including the value in the exception string though!
| particle=particle) | ||
| return velDist | ||
| # setting up integration from -10*vTh to 10*vTh, which is close to Inf | ||
| infApprox = 10 * vTh |
There was a problem hiding this comment.
I'm trying to think of a way we could check whether this is a good approximation (aside from calculating this again for 20 vTh) and raise a warning if it turns out not to be one.
There was a problem hiding this comment.
I think this ends up being something like 50 or 100 e-foldings of decay, since we have roughly exp(-v**2 / vTh**2) => exp(-(10*vTh)**2 / vTh**2) => exp(-100), times some coefficients. In any case, this test works for the 1D case and passes the isclose to 1 test with default errors.
There was a problem hiding this comment.
... 50 e-foldings? Okay. Nevermind the warning. 😆
| # if I recall correctly, astropy.units may have trouble with | ||
| # accounting for units in integrals, so this step currently fails when | ||
| # tplquad tries to convert a quantity with units to a scalar. | ||
| # integrating, this should be close to 1 |
There was a problem hiding this comment.
Turn this into a simple float via velDist3D.si.value, then multiply the result by the proper units? infApprox would also have to be changed.
| lambda z, y: -infApprox, | ||
| lambda z, y: infApprox, | ||
| epsabs=1.49e-08*u.m/u.m, | ||
| epsrel=1.49e-08*u.m/u.m) |
There was a problem hiding this comment.
I'm curious, why this number in particular? Also, the units seem off, could this perhaps be why tplquad failed here? If we're integrating velocities, don't we need u.m/u.s?
There was a problem hiding this comment.
Those are just the default values in tplquad 😉 , I put them in explicitly to check if adding dimensionless units would help the test behave correctly. Sometimes this works according to the "known issues" with astropy.units pages I linked below.
The final results from integration should be a probability, which is dimensionless, therefore the errors on the final value are dimensionless. u.m/u.m was just a shortcut for adding a dimensionless unit.
Here are some references for why the integration doesn't behave well with astropy.units. Here and here.
There was a problem hiding this comment.
Oh right, it's an error on the final value! I forgot about that, you're quite right.
Well, we can still go for floats on the internal side!
|
I've fixed the issues with astropy.units and tplquad not playing nicely. The fix was relatively simple, but pytest now runs really slowly or stalls. When it does work, just one of the 3D tests takes about 150 seconds. I've tried profiling with cProfile, but it also stalls. I'm attaching a partial cProfile which may give some insights on how to optimize. |
| epsrel=1.49e-08*u.m/u.m) | ||
| exceptStr = "Integral of distribution function should be 1." | ||
| assert np.isclose(integ.value, 1), exceptStr | ||
| epsabs=0.0, |
There was a problem hiding this comment.
probably should change this from 0.0 to some larger value as I'm not sure how QUADPACK handles 0.0 absolute error (whether infinitely small or just ignored).
|
I just tested out a version of Maxwellian_velocity_3D() with all the unit checks and conversions stripped out and passing additional parameters under args=() in tplquad rather than creating a partially-applied function. It runs in 1.6 seconds.... and gives 1.0000000181399045 as the results of the triple integral, which is pretty damn close. So perhaps the speed of unit conversions, mentioned in #131 , is a valid concern. Need to dig in and figure out what was killing the optimizations in tplquad. |
| assert np.isclose(integ.value, 1), exceptStr | ||
| epsabs=0.0, | ||
| epsrel=5.0e-1) | ||
| # value returned from tplquad is (integral, error), we just need the 1st |
There was a problem hiding this comment.
Or we could use the returned error in line 102-105, as atol? Does that make sense?
There was a problem hiding this comment.
I'm not sure. I don't think the error here is the error on the final value, but rather the error in how the function slices up the bounds into individual units to be calculated, or the error on these local approximate integrals.
There was a problem hiding this comment.
In any case using epabs=5.0e-1 gave me quite a reasonable value on the final integral value when I stripped the Maxwellian function of units. This was the aforementioned 1.0000000181399045.
|
I'd have to bet (70% of my confidence) on unit conversions. @Cadair mentioned them being expensive operations, to be used sparingly. Still... it might be the partial function application. |
namurphy
added
the
plasmapy.formulary
|
Here are some ideas on how to resolve this
Some further optimization could be done by passing vThermal as an argument. For any task, such as integration, where the temperature is the same, vThermal will be the same. The way the function is written right now it needlessly calculates vThermal on each call. Some of this will depend on what part of the code is actually slowing down tplquad. My inclination is towards option 3. |
|
A reasonable inclination, that. Avoiding recalculating the thermal velocity
is also a good idea!
|
…on functions. These function also include optional drift velocities/speeds. I've included tests for these new functions, but the 3D tests fails because scipy's QUADPACK based triple integral doesn't play nice with units.
…dded new functions to __init__.py. Noted in comments that velocities are actually relative particle velocities.
…ng nicely with astropy.units. This was done by applying units to inputs and stripping units of output from partially applied distribution function within the test.
…d. Code needs to be cleaned up.
…ng unit and unitless options, and an option to pass thermal velocity to avoid redundant calculations.
…0*u.m/m.s was passed as an argument. This made the units squared, which then caused failed unit conversions. Also fixed docstring values.
lemmatum
force-pushed
the
distributionFunctions
branch
from
709bcec
to
44e885e
Compare
|
Apologies for all the pushes. I was learning how to rebase since upstream's distribution.py was updated with edits to the docstrings. |
plasmapy/physics/distribution.py
Outdated
| coeff = (vThSq * np.pi) ** (-1 / 2) | ||
| expTerm = np.exp(-vSq / vThSq) | ||
| distFunc = coeff * expTerm | ||
| except Exception: |
There was a problem hiding this comment.
Code coverage whines about this one, so either add some test with pytest.raises(ValueError) or simply plot down a # coveralls: ignore after the Exception: to skip testing this one!
There was a problem hiding this comment.
Or just remove the try-except block? I don't think it's doing much here.
There was a problem hiding this comment.
We could do that as well. I'm not 100% sure why they pop up everywhere.
There was a problem hiding this comment.
I just started by imitating other code I've seen in this package, but it seems like these try-except blocks are mostly useless for such simple functions.
| # need to be assigned units | ||
| if Vx_drift == 0: | ||
| if not isinstance(Vx_drift, astropy.units.quantity.Quantity): | ||
| Vx_drift = Vx_drift * u.m/u.s |
There was a problem hiding this comment.
The three Vdirection_drift = ... lines don't get covered by test. Adding another small one with a floating point value would solve that!
There was a problem hiding this comment.
This refers to tests in test_distribution.py, correct?
plasmapy/physics/distribution.py
Outdated
| V_drift = V_drift.to(u.m/u.s) | ||
| # convert temperature to Kelvins | ||
| T = T.to(u.K, equivalencies=u.temperature_energy()) | ||
| if vTh == -1: |
There was a problem hiding this comment.
This is a neat idea! Would adding np.nan instead of -1 make sense here, to signify that it's a nonsensical value?
There was a problem hiding this comment.
Sure. Just tested it and it does work. Just had to switch from == to np.isnan(vTh) as the Truthiness test.
plasmapy/physics/distribution.py
Outdated
| vTh = vTh.to(u.m/u.s) | ||
| elif vTh == -1 and units == "unitless": | ||
| # assuming unitless temperature is in Kelvins | ||
| vTh = (thermal_speed(T*u.K, particle=particle, method="most_probable")) |
There was a problem hiding this comment.
This one doesn't get tested :(
plasmapy/physics/distribution.py
Outdated
| vTh = (thermal_speed(T, particle=particle, method="most_probable")) | ||
| elif vTh != -1: | ||
| # check units of thermal velocity | ||
| vTh = vTh.to(u.m/u.s) |
There was a problem hiding this comment.
If the velocity variable can be used as a velocity (as the tests are passing later), is this check (elif block from line 338) necessary? Not 100% sure.
There was a problem hiding this comment.
Do you mean that if vTh gets added or multiplied to some other value that has units then astropy will throw an error? I was just trying to be explicit, but that might be a better way to go, less clutter.
There was a problem hiding this comment.
Well, if it's convertible to u.m/u.s, it's not going to throw any errors later even if you don't do the explicit check here. But explicit is neat, I was just asking out of curiosity. Let's keep it this way :)
plasmapy/physics/distribution.py
Outdated
| expTerm = np.exp(-vSq / vThSq) | ||
| distFunc = coeff1 * coeff2 * expTerm | ||
| except Exception: | ||
| raise ValueError("Unable to get distribution function.") |
There was a problem hiding this comment.
This block doesn't get tested either. Could ignore it.
| # need to be assigned units | ||
| if Vx_drift == 0: | ||
| if not isinstance(Vx_drift, astropy.units.quantity.Quantity): | ||
| Vx_drift = Vx_drift * u.m/u.s |
There was a problem hiding this comment.
These 3 lines don't get tested either. Same situation as before, it seems!
|
Oh yeah, I forgot to mention. Mad props on getting that test down to 1s! That's really amazing. |
|
It wasn't too difficult to optimize! Thank those poor FORTRAN programmers who wrote QUADPACK 😜 |
|
Their sacrifice shan't be forgotten. [*] |
…s. This may cause breakages in other parts of plasmapy.
…or each Test class.
…ions of units, unitless, passing vTh, default vTh, V_drifts, default V_drifts.
… correct units and use correct particle defined in setup_method.
|
I added a bunch of tests under |
|
Looks like coveralls are happy with the 3D Maxwellian speed function tests. I'll make similar tests for the other functions. |
…s input permutations (vTh, units, drift) for 1D velocity, 1D speed, and 3D velocity distribution functions. These are similar to tests from previous commit for 3D speed distribution function.
|
Added tests for other functions 🤞 |
| T=self.T_e, | ||
| particle=self.particle, | ||
| units="units") | ||
| errStr = (f"Distribution function should be {self.distFuncTrue} " |
There was a problem hiding this comment.
Bonus points for f-string usage
Added and modified 1D and 3D Maxwellian speed and velocity distribution functions. These functions also include optional drift velocities/speeds. I've included tests for these new functions, but the 3D tests fail because scipy's QUADPACK based triple integral doesn't play nice with units. I'm open to suggestions on how to work around this!
Contributions towards #101 and #141