Basic simulation capability [WIP] by SolarDrew · Pull Request #33 · PlasmaPy/PlasmaPy

SolarDrew · 2017-05-16T17:01:08Z

Closes #7 and provides a simple-to-use interface for setting up and running simulations. What's here at the moment is based largely on a similar previous project, which can be found here, but I'm happy to deviate from it based on input from others.

Below are some snapshots of density from a demo 2D simulation with a plasma blob initally in the middle and diffusing outwards. It's technically solving the MHD equations but with B = 0. As you can see, I'm getting some instabilities even without magnetic field, and the less said about the boundary conditions the better.
Also speed is obviously an issue at the moment.

Pynite being my previous attempt at a python plasma simulation code.

pep8speaks · 2017-05-16T17:01:12Z

Hello @SolarDrew! Thanks for updating the PR.

In the file plasmapy/classes/simulation.py, following are the PEP8 issues :

Line 22:80: E501 line too long (102 > 79 characters)

In the file plasmapy/numerical/spatial_solvers.py, following are the PEP8 issues :

Line 132:80: E501 line too long (80 > 79 characters)

Comment last updated on August 13, 2017 at 13:34 Hours UTC

StanczakDominik · 2017-05-17T07:07:29Z

+    def density(self, density: u.kg/u.m**3):
+        """Sets the simulation's density profile to the specified array.
+        Other arrays which depend on the density values, such as the kinetic pressure,
+        are then redefined automatically.


That's clever, I like it!

StanczakDominik · 2017-05-17T07:11:00Z

+
+        assert density.shape == self.grid_size, \
+            'Specified density array shape does not match simulation grid.'
+        self._density = density


Doesn't the previous necessitate self._pressure = self.pressure() after this line here?

I think it doesn't, but because the comment is wrong rather than the implementation.

Change pressure to energy, fix grid inconsistancies and pep8ify

Mostly for the purposes of keeping the pep8 bot happy

StanczakDominik

Had a brief look to stay up-to-date on what you're doing here. It looks pretty awesome!

StanczakDominik · 2017-06-01T08:20:53Z

 import astropy.units as u
-
-mu0 = np.pi * 4.0e-7 * (u.newton / (u.amp**2))
+from astropy.utils.console import ProgressBar


Right, because of course astropy would have a built-in progress bar for all things astronomical.

Yeah, astropy is amazing :D

StanczakDominik · 2017-06-01T08:21:12Z

        self._magnetic_field = np.zeros((3, *self.domain_shape)) * u.T

+        # Connect a simulation object for simulating
+        self.simulation_physics = MHDSimulation()


Love that you're keeping this generic!

StanczakDominik · 2017-06-01T08:31:58Z

+            raise ValueError("Either max_time or max_its must be set.")
+
+        if np.isinf(max_time):
+            pb = ProgressBar(max_its)


Perhaps we could raise a warning here if both max_time and max_its are set?

Good idea. But the question then becomes how to procede in a predictable way if there are two possible end 'times'. Stop at whichever one is sooner I suppose, but I feel like that's going to lead to unexpected behaviour if users aren't careful, so maybe raising an exception would be better. I'm happy to be convinced otherwise though.

Sure, an exception does make more sense than a warning.

StanczakDominik · 2017-06-01T08:32:40Z

+
+        # Collect equations into a nice easy-to-use list
+        # self.equations = [self._ddt_density, self._ddt_momentum,
+        #                  self._ddt_energy, self._ddt_magfield]


Wow, I love this idea.

Thanks :) I found in the earlier version of this code that it's an easy and generalisable way to get at the functions you need.

So that it can get at the plasma variables, grid size, etc.

Also associated things

SolarDrew · 2017-06-01T13:29:38Z

+
+
+class MHDSimulation():
+    """Physics class for magnetohydrodynamics.


Cadair · 2017-06-01T13:37:16Z


-mu0 = np.pi * 4.0e-7 * (u.newton / (u.amp**2))

 class Plasma():


You don't need the () in Python 3.

But it looks so strange without them....

SolarDrew · 2017-07-23T15:56:42Z

Update: so before NAM I managed to get some of the hyperviscosity implemented and running. With this in place, I've managed to run the Riemann shock tube test and the output looks like this:

For those who aren't familiar with it, it should look something like this:

(image from here).

So there are still definitely some issues that need fixing, mostly with stability, but it's getting there.

namurphy · 2017-07-24T15:16:48Z

I haven't looked into this much, but sympy.utilities.codegen is probably useful (for some definition of useful). There was a tutorial on automatic code generation with SymPy at SciPy 2017. (These are what I mentioned in our last telecon.)

In the intermediate future we may also wish to have some finite element modeling capabilities in here too, such as with packages like FEniCS. This could be really helpful with non-trivial geometry like tokamaks with non-circular cross sections.

Used for allowing the viscous tensor to interact with things

To more easily see if/when calculations are going wrong

Turns out this was breaking the tests because it should have been using just the dx value for the appropriate axis. Didn't show up in 1D examples.

Was previously default 5/3, now set to match the example I'm trying to replicate

StanczakDominik · 2017-08-13T18:34:54Z

+    dot = np.array([sum([vec[i] * tensor[i, 0, ...] for i in dims]),
+                    sum([vec[i] * tensor[i, 1, ...] for i in dims]),
+                    sum([vec[i] * tensor[i, 2, ...] for i in dims])]) \
+        * vec.unit * tensor.unit


Perhaps, if speed is an issue, this part could be optimized somehow. Two possibilities that come to mind are some crazy numpy broadcasting tricks (I'm not saying I see how to do this just yet, though!) and JITting this part via Numba.

Numba doesn't play nice with Coveragepy, though. 😞

Yeah, I feel like this could be speeded up by doing it the 'proper' numpy way. At the moment most of the viscous tensor stuff is still a bit fluid (boom boom) because I don't really fully understand and I haven't figured out how to do it with broadcasting. Premature optimisation, right...?

Of course! Still, let's keep it in mind for later :)

I think I might have done something like this with numpy.einsum or whatever. I'll try to figure out how to do this via broadcasting a bit later.

Definitely :) Cool, thanks :D

StanczakDominik · 2017-08-13T18:37:12Z

    t0 = datetime.datetime.now()
+    for maxt in np.arange(0, 0.21, 0.02):
+        print('- Running Simulation... ', end='')
+        riemann.simulate(max_time=maxt*u.s)


As I see it, we're running all pictures from 0 to t where t increases in the arange, right? Perhaps it would be easier on the computation to rewrite simulation to avoid that - go from 0 to t_1, take a snapshot.

In all likelihood this is something you've already thought of, but I just thought I'd raise it out in the open if anyone wants to hack away at implementing that.

Huh, this is outdated now. Okay, git trees are magic.

Yep, they are :P You're right, in the long term I'd like to have the Simulation object itself be able to save a snapshot so you don't need this sort of thing, but I haven't got around to it yet so I'm still using this hack.

StanczakDominik · 2017-08-13T18:38:48Z

+        riemann.simulate(max_time=maxt*u.s)
+        print("Simulation complete")
+
+        fig, ax = plt.subplots(2, 2)


I always liked doing it like this:

fig, ((density_ax, velocity_ax), (energy_ax, pressure_ax)) = plt.subplots...

Yeah, this is a bit nicer. Good call, thanks.

StanczakDominik · 2021-03-30T18:55:05Z

Sorry @SolarDrew, I think it's time to mark this one as stale. 😉

SolarDrew · 2021-03-31T08:59:12Z

😭 Fair enough. Some day I'll come back to this and make it work!

Port in a load of stuff from pynite

b63fb6a

Pynite being my previous attempt at a python plasma simulation code.

Small handful of pep8 fixes

51bb91f

StanczakDominik approved these changes May 17, 2017

View reviewed changes

namurphy added the status: in progress label May 18, 2017

SolarDrew added 8 commits May 23, 2017 19:05

Port in more stuff from pynite and pep8 it

dcced0b

Remove some stuff that isn't need yet/is wrong

33c28b6

Update spinxy tag at top of plasma.py

35ddfa4

Reinstate energy as parameter and define dot func

37f69aa

Changes to plasma test functions

d76bf05

Change pressure to energy, fix grid inconsistancies and pep8ify

Add B field strength and alfven speed properties

41ad7bd

Fix tests for < 3D

d48f47f

Variable name shortenings

e474cc1

Mostly for the purposes of keeping the pep8 bot happy

StanczakDominik approved these changes May 29, 2017

View reviewed changes

SolarDrew added 7 commits May 30, 2017 19:16

Merge branch 'master' into simulation

d1054d3

Add a simulate() method to the plasma class

0d327cd

Basic test function for simuation

8e569fb

Minor changes to Plasma to do with introduccing a simulation object

8694b88

Basic shell for a Simulation object

2326931

Tweak imports in plasma.py

4818eca

Move dot function from plasma.py to simulation.py

ecedbd6

StanczakDominik approved these changes Jun 1, 2017

View reviewed changes

SolarDrew added 6 commits June 1, 2017 10:20

Add plasma attribute to simulation object

b0f8d22

So that it can get at the plasma variables, grid size, etc.

A few minor rearrangements

367e7e7

Bring in RK4 solver from pynite as time_stepper

83b2820

Also associated things

Vector operators and minor tweaks

0092012

Cut and paste spatial solver from pynite

cddac44

Tweak simulation test to make it make sense

7c5b1e3

SolarDrew commented Jun 1, 2017

View reviewed changes

Cadair reviewed Jun 1, 2017

View reviewed changes

SolarDrew added 10 commits July 30, 2017 11:57

Minor correction and formatting changes

0065071

Add diffusive component of energy artificial viscosity term

89258de

Function to calculate vector-tensor dot product

d4b8cd5

Used for allowing the viscous tensor to interact with things

Another visc function and tweak c parameters

844c2ea

Add more terms to energy hyperviscosity

57943ff

Re-introduce plotting of snapshots

b795f08

To more easily see if/when calculations are going wrong

Correct how viscosity deals with dx

425d95f

Turns out this was breaking the tests because it should have been using just the dx value for the appropriate axis. Didn't show up in 1D examples.

Correct adiabatic index

2b16500

Was previously default 5/3, now set to match the example I'm trying to replicate

Up grid resolution to match example in paper

c37df82

Merge branch 'master' of github.com:PlasmaPy/PlasmaPy into simulation

791f3c8

StanczakDominik reviewed Aug 13, 2017

View reviewed changes

StanczakDominik mentioned this pull request Oct 25, 2017

Reimplement Plasma numerical solvers via Numba for potential speed increases #137

Closed

StanczakDominik self-assigned this Dec 23, 2017

StanczakDominik force-pushed the master branch 4 times, most recently from 02e758d to 749a6bd Compare April 29, 2018 16:04

namurphy mentioned this pull request Jul 20, 2018

Write PLEP on how to implement plasma simulation capabilities PlasmaPy/PlasmaPy-PLEPs#17

Open

StanczakDominik added this to the Future milestone Oct 16, 2018

StanczakDominik removed their assignment Oct 16, 2018

namurphy mentioned this pull request Dec 18, 2018

PLEP on next generation plasma simulator PlasmaPy/PlasmaPy-PLEPs#24

Open

StanczakDominik removed Enhancement labels Jun 25, 2019

StanczakDominik closed this Mar 30, 2021

namurphy added plasmapy.simulation Related to the plasmapy.simulation subpackage and removed status: in progress labels May 23, 2023



		class MHDSimulation():
		"""Physics class for magnetohydrodynamics.


		mu0 = np.pi * 4.0e-7 * (u.newton / (u.amp**2))

		class Plasma():

Conversation

SolarDrew commented May 16, 2017 • edited by StanczakDominik Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

pep8speaks commented May 16, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Comment last updated on August 13, 2017 at 13:34 Hours UTC

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

StanczakDominik left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

SolarDrew Jun 1, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

SolarDrew commented Jul 23, 2017

Uh oh!

namurphy commented Jul 24, 2017

Uh oh!

StanczakDominik Aug 13, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

StanczakDominik commented Mar 30, 2021

Uh oh!

SolarDrew commented Mar 31, 2021

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

SolarDrew commented May 16, 2017 •

edited by StanczakDominik

Loading

pep8speaks commented May 16, 2017 •

edited

Loading

SolarDrew Jun 1, 2017 •

edited

Loading

StanczakDominik Aug 13, 2017 •

edited

Loading