Add standalone demo.

[jacklovell]

This demo uses the raw files from Andrey Pshenov's PhD thesis
simulations, converted to SOLPS 5 format by Vlad Neverov. The demo
contains similar functionality to the existing AUG/raw, MDSplus and
balance.nc demos. Additionally, it demonstrates how to add emission
models to the plasma and image the plasma with a visible camera,
providing a complete workflow example which should be useful to new
Cherab users and familiar to experienced users new to the SOLPS
package.

N.B. The simulation output files are rather large, at 31 MB, so they're version controlled in a compressed zip file which is much smaller at 5.5 MB. The user will need to extract the files before running the demo: instructions for how to do this are provided in the demo if the files have not already been extracted.

Discussion on this PR can also be used to decide whether we want to modify the other demos. For example, I've replaced the nested for loops for sampling with Cherab's own sampling functions in this PR, which should be considered best practice. Similarly, the other demos don't actually use any emission models: they reproduce the functionality of b2plot but don't exploit Cherab's unique capabilities, so perhaps we should add some imaging to them too.

[jacklovell]

@vsnever Github won't let me request a review from you on this, but I'd be interested in your input too.

[vsnever]

This is a good demonstration. However, I experienced problems with interactive mode of matplotlib. Whether I run this demo with python or ipython, the figures hang at input('Press Return to continue...'). Maybe the problem is specific to my setup, but giving matplotlib the time to render the figures solves the issue:

plt.pause(0.01)
input('Press Return to continue...')

The final image required more time to render for some reason:

plt.pause(0.1)
input('Press Return to finish...')

Also, usually I run scripts from code editors like Sublime Text, and I think other users do too. The request for the user input in the middle of the script requires me to run this demo from console or to install an additional plugin. The rest of the script adds only a single figure, and I think closing other figures before the rendering starts is unnecessary. So, replacing

input('Press Return to continue...')
plt.close('all')

with just a pause

plt.pause(0.01)

would help. Also in this case, the user can explore the figures while the final image is rendering.

Below are some other potential improvements.

[vsnever]

Why ask the user to unzip files manually when the script can do it automatically?

Suggested change

	try:
	print('Loading simulation...')
	sim = load_solps_from_raw_output(simulation_directory)
	except RuntimeError:
	print('ERROR: simulation files not found. Please extract the raw '
	'files from the zip archive before running this demo:')
	print('> unzip {0}/generomak_raw_files.zip -d {0}/generomak_raw_files'.format(demos_directory))
	sys.exit(1)
	try:
	print('Loading simulation...')
	sim = load_solps_from_raw_output(simulation_directory)
	except RuntimeError:
	print('Extracting files...')
	with zipfile.ZipFile('generomak_raw_files.zip', 'r') as zip_ref:
	zip_ref.extractall(simulation_directory)
	sim = load_solps_from_raw_output(simulation_directory)

[Mateasek]

This is a good idea, zipfile is in the python standard library, so why not to use it.

What about doing the extraction in the setup.py script? The similar way as the installation of the OpenADAS repo? We then don't have to care about the files not being extracted from the archive in every demo script.

[jacklovell]

We've had issues in the past with re-installing and re-downloading the OpenADAS files where we get some cryptic error, which put me off adding this step to setup.py. I generally prefer the demos to be "pure": i.e. the demo should not modify the surrounding file system without being asked. But I see that some demos in other packages do write output files, and this isn't too far from that.

As there are two votes for this, I'll happily make the change.

[Mateasek]

I understand that the unpacked files are large and they will be tracked by git if upacked, but isn't it an acceptable cost to avoid all the problems we see with having archives?

[vsnever]

Whether the user unzip files manually, or the script does it, in both cases we have the same problem. Namely, if we update the files in the archive, this demo will use the old versions if they were unzipped previously.

[jacklovell]

Yes, that's true. We have a similar issue with compiling the pyx files: a user must run ./dev/build.sh to ensure the latest versions are used. But I don't imagine we'll be updating these data files very often, if ever, so perhaps a one-time hit of c.30MB added to the repository isn't the end of the world.

@CnlPepper do you have any advice for this?

[CnlPepper]

I'm with Jack on this. I think it would be better to inform the user the files are not available when they run the demo. The error message can then describe the process to unzip the data. I'd only revisit automation if we get a lot of issues raised.

[CnlPepper]

That said, 30Mb isn't a killer, you could just include the files directly in the repo. We do that for Raysect. All the demo meshes are in the repo.

[jacklovell]

OK, for simplicity I'll include the uncompressed files and remove all logic related to decompression. That should make the demo easier to follow.

[vsnever]

With the new interpolators, these distributions can be sampled easier.

Suggested change

	xsamp, _, zsamp, ne_samples = sample3d(plasma.electron_distribution.density,
	(xl, xu, nx), (0, 0, 1), (zl, zu, ny))
	ne_samples = ne_samples.squeeze()
	te_samples = sample3d_grid(plasma.electron_distribution.effective_temperature,
	xsamp, [0], zsamp).squeeze()
	h0_samples = sample3d_grid(h0.distribution.density, xsamp, [0], zsamp).squeeze()
	h1_samples = sample3d_grid(h1.distribution.density, xsamp, [0], zsamp).squeeze()
	c0_samples = sample3d_grid(c0.distribution.density, xsamp, [0], zsamp).squeeze()
	c1_samples = sample3d_grid(c1.distribution.density, xsamp, [0], zsamp).squeeze()
	c2_samples = sample3d_grid(c2.distribution.density, xsamp, [0], zsamp).squeeze()
	c3_samples = sample3d_grid(c3.distribution.density, xsamp, [0], zsamp).squeeze()
	c4_samples = sample3d_grid(c4.distribution.density, xsamp, [0], zsamp).squeeze()
	c5_samples = sample3d_grid(c5.distribution.density, xsamp, [0], zsamp).squeeze()
	c6_samples = sample3d_grid(c6.distribution.density, xsamp, [0], zsamp).squeeze()
	# Magnitude of main ion velocity vector.
	h1_velocity = samplevector3d_grid(h1.distribution.bulk_velocity, xsamp, [0], zsamp).squeeze()
	xsamp, zsamp, ne_samples = sample2d(sim.electron_density_f2d, (xl, xu, nx), (zl, zu, ny))
	te_samples = sample2d_grid(sim.electron_temperature_f2d, xsamp, zsamp)
	h0_samples = sample2d_grid(sim.species_density_f2d[('hydrogen', 0)], xsamp, zsamp)
	h1_samples = sample2d_grid(sim.species_density_f2d[('hydrogen', 1)], xsamp, zsamp)
	c0_samples = sample2d_grid(sim.species_density_f2d[('carbon', 0)], xsamp, zsamp)
	c1_samples = sample2d_grid(sim.species_density_f2d[('carbon', 1)], xsamp, zsamp)
	c2_samples = sample2d_grid(sim.species_density_f2d[('carbon', 2)], xsamp, zsamp)
	c3_samples = sample2d_grid(sim.species_density_f2d[('carbon', 3)], xsamp, zsamp)
	c4_samples = sample2d_grid(sim.species_density_f2d[('carbon', 4)], xsamp, zsamp)
	c5_samples = sample2d_grid(sim.species_density_f2d[('carbon', 5)], xsamp, zsamp)
	c6_samples = sample2d_grid(sim.species_density_f2d[('carbon', 6)], xsamp, zsamp)
	# Magnitude of main ion velocity vector.
	h1_velocity = samplevector2d_grid(sim.velocities_cylindrical_f2d[('hydrogen', 1)], xsamp, zsamp)

However, if the goal of this demo is to show how the plasma object is built from the simulation, then it's better to leave it as it is.

[Mateasek]

This really depends on the the aim of the demo. If this is to be aiming at the solps package, then I would use what @vsnever suggests. The option with plasma distribution can be used in a demo more aiming at the plasma class.

[jacklovell]

Yes, this is open to interpretation. Ultimately users are likely to be using Cherab for its plasma object, and not just an alternative to b2plot. But verifying the SOLPS data has been loaded correctly is also important.

Perhaps we should have two similar demos, one which samples from the simulation object and one which samples from the plasma object. What do you think?

[vsnever]

I would show different ways of sampling in the same demo and add comments explaining that there are a few ways to do this.

[jacklovell]

OK, sounds good. I'll put in something like np.all(sim_samples == plasma_samples) to demonstrate, rather than duplicating all the plots. And I won't use both methods for all the carbon ionisation stages, but can show the two methods for ne, Te, nC0 for example.

[vsnever]

This is very optional, but the units can be shown in the colorbar instead of the title.

Suggested change

	plt.colorbar()
	plt.xlim(xl, xu)
	plt.ylim(zl, zu)
	plt.title('Electron density [m-3]')
	plt.colorbar(label='m$^{-3}$')
	plt.xlim(xl, xu)
	plt.ylim(zl, zu)
	plt.title('Electron density')

[jacklovell]

If this were for a publication I'd agree with you. For a simple code demo I don't think it's as important. I've just tried both ways and I prefer being able to quickly see the units right next to the quantity in the title rather than having to look half way down the colour bar for them.

[vsnever]

When and if raysect/source#343 is merged, I'll make a separate demo to show how the rendering time of this can be reduced by ~50 times.

atomic_data = OpenADAS(permit_extrapolation=True)
emission = np.zeros((nx, 1, ny, camera.spectral_bins))
direction = Vector3D(0, 0, 1)
for species in plasma.composition:
    if species.charge < species.element.atomic_number:
        model = TotalRadiatedPower(species.element, species.charge, plasma=plasma, atomic_data=atomic_data)
        for i, x in enumerate(xsamp):
            for j, z in enumerate(zsamp):
                point = Point3D(x, 0, z)
                spectrum = Spectrum(camera.min_wavelength, camera.max_wavelength, camera.spectral_bins)
                emission[i, 0, j] += model.emission(point, direction, spectrum).samples

rgc = RegularGridCylinder(emission, spectrum.wavelengths, xu, zu - zl, radius_inner=xl, continuous=True,
                          parent=world, transform=translate(0, 0, zl))

In case you want to test this now, do not connect the plasma object to the world and increase the camera resolution to at least 512 x 512 to see a 50x speedup.

[Mateasek]

I think we have to decide where to pile the demo scripts. If we want to use the generomak in Cherab core as example data (which I thought was the purpose) then observation examples should be moved to the core package demos. The camera observation @jacklovell made could be kept in the solps package as an example for the sake of completeness.

[jacklovell]

This would be better as one of a series of "performance optimisation" demos. There's a lot of added complexity in the fast emitter approach, which is a price worth paying if your render takes hours or days but is probably going to be overwhelming to a new Cherab user who doesn't mind waiting 1 minute for a render. I think we should keep the SOLPS demos as simple as possible and illustrate how the package fits into the standard Cherab workflow.

[vsnever]

Sure, it is definitely not for this demo.

[vsnever]

I almost forgot to add that the actual version of fort.44 file in the generomak demo is 20081111. I manually changed it to 20130210, so the fort44_2013 parser could read it. However, the total radiation of neutral atoms is read incorrectly (it's not used in this demo). I'll try to create a more universal fort44 parser soon, and then we will need to change the version of fort.44 file back to 20081111.

[Mateasek]

Hi @jacklovell, thanks for the demo. I would suggest following changes. I would move the zip archive into a folder named for example "data". In this folder we could add a "utility.py" which could contain methods handling extraction of the zip archive as @vsnever sugeests. This method can be called from the setup.py, if we want it. This would I think make a cleaner structure and open the space for us to have more demo simulations in future. Also it would allow to remove the checks about the files existence from every demo we will use the simulations in.

Also I would separate the demo into two actually. I would make separate demos for the camera and the profile plotting.

The structure could look like this (names of folders and files are optional):

-demos/
    -plot_solps_profiles.py
    -visible_camera.py
    -data/
        -utility.py
        -generomak/
            *solps data files*

[Mateasek]

Is the input and closing really needed? If you remove this the user could be already enjoying the profile plots while the observation is running.

[jacklovell]

I put it in because every 15 seconds when the observation updates the plot it shifts focus to the camera plot. I decided to wait until the user is done admiring the profiles before starting the observation. But the suggestion to separate the demo into two is a good one which I'll implement, so closing the plots will become unnecessary.

[vsnever]

Just one more comment. I think this demo will make more sense from a physical point of view if RadiancePipeline2D is used instead of the default RGBPipeline2D. Also it might be worth setting camera.spectral_bins = 1 for better performance.

[jacklovell]

This is a good demonstration. However, I experienced problems with interactive mode of matplotlib. Whether I run this demo with python or ipython, the figures hang at input('Press Return to continue...'). Maybe the problem is specific to my setup, but giving matplotlib the time to render the figures solves the issue:

plt.pause(0.01)
input('Press Return to continue...')

The final image required more time to render for some reason:

plt.pause(0.1)
input('Press Return to finish...')

Also, usually I run scripts from code editors like Sublime Text, and I think other users do too. The request for the user input in the middle of the script requires me to run this demo from console or to install an additional plugin. The rest of the script adds only a single figure, and I think closing other figures before the rendering starts is unnecessary. So, replacing

input('Press Return to continue...')
plt.close('all')

with just a pause

plt.pause(0.01)

would help. Also in this case, the user can explore the figures while the final image is rendering.

Below are some other potential improvements.

This is strange. I haven't noticed this behaviour myself, and am wary of spraying plt.pause commands everywhere as it then breaks Jupyter notebooks. But I agree with your other suggestion to split this into two demos and remove the "Press return to continue" steps, which should hopefully solve this anyway.

[vsnever]

This is strange. I haven't noticed this behaviour myself, and am wary of spraying plt.pause commands everywhere as it then breaks Jupyter notebooks. But I agree with your other suggestion to split this into two demos and remove the "Press return to continue" steps, which should hopefully solve this anyway.

The behaviour of matplotlib interactive mode is kind of buggy. There are several threads on stackoverflow about this and some of them are very recent:
https://stackoverflow.com/questions/50490426/matplotlib-fails-and-hangs-when-plotting-in-interactive-mode
https://stackoverflow.com/questions/64043973/matplotlib-interactive-mode-wont-work-no-matter-what-i-do

However I do not insist on this change, the user can always turn off the interactive mode if it does not work.

[jacklovell]

I almost forgot to add that the actual version of fort.44 file in the generomak demo is 20081111. I manually changed it to 20130210, so the fort44_2013 parser could read it. However, the total radiation of neutral atoms is read incorrectly (it's not used in this demo). I'll try to create a more universal fort44 parser soon, and then we will need to change the version of fort.44 file back to 20081111.

Ah, yes, I've tried plotting neutrad and got some highly suspicious 1e23 W/m3 emissivity. Thanks for pointing this out. Eneutrad will stay out of the demo until the fort.44 parser can handle it.

[jacklovell]

Just one more comment. I think this demo will make more sense from a physical point of view if RadiancePipeline2D is used instead of the default RGBPipeline2D. Also it might be worth setting camera.spectral_bins = 1 for better performance.

I quite like the colour: it's more engaging than a black and white plot. I've just tried adding a PowerPipeline2D to the camera's pipelines, and the resulting images are surprisingly similar (other than the colour). I presume this is because most of the higher charge states are at very low densities in the edge region, and you're not getting the large volumetric VUV emission from the core as it's not present in the simulation. These edge simulations are where visible light tends to be dominant, so an RGB camera does a pretty good job of capturing most of the plasma emission.

I could leave both pipelines in: you'd see a visible image and a total emission that way. It makes no difference to the render time, which is dominated by evaluating the emission models.

[vsnever]

I quite like the colour: it's more engaging than a black and white plot. I've just tried adding a PowerPipeline2D to the camera's pipelines, and the resulting images are surprisingly similar (other than the colour). I presume this is because most of the higher charge states are at very low densities in the edge region, and you're not getting the large volumetric VUV emission from the core as it's not present in the simulation. These edge simulations are where visible light tends to be dominant, so an RGB camera does a pretty good job of capturing most of the plasma emission.

I could leave both pipelines in: you'd see a visible image and a total emission that way. It makes no difference to the render time, which is dominated by evaluating the emission models.

But the total radiated power is not resolved over wavelength, or am I wrong?

[jacklovell]

Ah yes, you're correct. It's spread uniformly over the spectrum. Fair point, I'll remove RGB.

[jacklovell]

I've made changes based on reviewers' suggestions, and pushed a new version.

I'd like to add fort.44 2008 compatibility to the Eirene data reader and change the version number back in the Generomak data file, so that I can include the total radiated power in the profile plots. Otherwise I hope this is ready to merge.

[mattngc]

Looks like a big improvement to me. Happy for you to merge when ready.

[jacklovell]

I've changed the fort.44 file's version header, so it now records that it's the 20081111 format. Since #40 was merged we can at least read all the relevant data for this demo correctly.

I've also made a small change to the plots: using SymLogNorm rather than LogNorm as the colour normalisation means we can show values of 0 in the plots (previously they were treated as NaN and not shown, even though the data did actually exist in the sampled arrays being plotted). To keep things tidy, this meant refactoring the plotting code out into a separate function which is just called once for each plot.