Fixes and improvements in SOLPSSimulation, SOLPSMesh and related functions

[vsnever]

This PR addresses the issues #8, #20, #31, #33, #34, #36 and #37.

Issue #8 is addressed in the following way:

• Safe to use SOLPSFunction2D and SOLPSVectorFunction2D classes are implemented. They replace SOLPSFunction3D and SOLPSVectorFunction3D. For 3D, use AxisymmetricMapper(SOLPSFunction2D) and VectorAxisymmetricMapper(SOLPSVectorFunction2D).

• SOLPSMesh and SOLPSSimulation can now be initialised without modification of hidden attributes;

• SOLPSMesh has new attributes for basis vectors, cell connection areas, indices of neighbouring cells and new functions to_poloidal() and to_cartesian() for converting vectors defined on a grid from/to (poloidal, radial)/(R, Z);

• species_list in SOLPSSimulation is a list of (name, charge) pairs now instead of a list of strings;

• Setting b_field and velocities vectors in curvilinear poloidal coordinates (e_pol, e_rad, e_tor) automaticaly sets b_field_cylindrical and velocities_cylindrical (e_R, e_φ, e_Z) and vice versa (note: e_tor = -e_φ);

• Incorrect calculation of basis vectors using cell centres is fixed (basis vectors are defined by left and bottom faces of the cell);

• The code for mesh.triangles and mesh.triangle_to_grid_map creation is vectorised (as all other code wherever possible);

• SOLPSSimulation now has setters for its properties;

• load_solps_from_mdsplus(), load_solps_from_raw_output() and load_solps_from_balance() do not brake SOLPSSimulation interface anymore;

• load_solps_from_mdsplus() and load_solps_from_raw_output() are made more robust. Now both support B2 standalone runs. Error handling is added to load_solps_from_mdsplus() for some optional attributes that may not be present on the server;

• Other small improvements.

Issue #20 is fixed. *_f2d() and *_f3d() interpolators are added for plasma parameters. For vectors _carteisan suffix is used instead of _f3d.

Issue #31 is fixed. Now if the neutral densities from Eirene are available, they replace the neutral densities from B2.

Issue #33 is fixed using lookup_isotope() and lookup_element() methods. The plasma composition is defined from nuclear charge, atomic mass number and electrical charge. This is more robust than parsing a string with species names.

Issue #34 is addressed by parsing the following additional quantities in load_solps_from_mdsplus() and load_solps_from_raw_output():

• Ion temperature;
• magnetic field (previously was read only from MDS);
• atom/ion fluxes (B2), neutral atom densities (Eirene, if available),
• neutral atom fluxes (Eirene, if available);
• neutral atom temperatures (Eirene, if available).

The fluxes are not stored in SOLPSSimulation, but used to calculate velocities.

Issue #36 is fixed. The parallel velocity and the B2 fluxes (in s^-1) defined on cell faces are used in b2_flux_to_velocity() to calculate velocities at cell centres. If neutral atom fluxes from Eirene (in m^-2 s^-1) defined at cell centres are available, they are used in eirene_flux_to_velocity() replacing the B2 ones.

Issue #37 is fixed. Now all arrays are row-major. However, the indexing is inverted. The correct index order now is: (radial, poloidal) or (species, radial, poloidal) for scalars and (vector component, radial, poloidal) or (species, vector component, radial, poloidal) for vectors. This should improve the performance in typical workflows.

Remaining issues:

• Electron velocities are still not read from SOLPS data.

• Molecular and total H-alpha emission are still not read from SOLPS data.

[vsnever]

Just updated the code by adding ion and neutral atom temperatures. The latter is not available in stand-alone B2.5 simulations.

[vsnever]

I renamed this PR because now it addresses not just #31 but in some degree #8 too.

Replacing the algorithm for poloidal basis vectors calculation in SOLPSMesh was more a bug fix than an improvement. The SOLPS mesh have break points due to the separatrix. In the old algorithm, which calculated the basis vectors using the cell centres, the basis vectors had wrong values at these break points.

[vsnever]

Just noticed that conversion from poloidal to Cartesian coordinates in mdsplus.py is incorrect. It's for orthogonal basis, but SOLPS uses curvilinear coordinates.

[vsnever]

Just noticed that conversion from poloidal to Cartesian coordinates in mdsplus.py is incorrect. It's for orthogonal basis, but SOLPS uses curvilinear coordinates.

Sorry, I'm wrong here. It converts to Cartesian, so it's OK.

[mattngc]

Hey Vlad, I don't have any strong opinions on these improvements. I don't use this module as much at the moment, so I think opinions from the main users are more important. I'll let @jacklovell and @Mateasek give feedback on these changes.

[vsnever]

Hey Vlad, I don't have any strong opinions on these improvements. I don't use this module as much at the moment, so I think opinions from the main users are more important. I'll let @jacklovell and @Mateasek give feedback on these changes.

Hi Matthew, thanks for reply. I need a couple more days to get this PR ready for review.

[vsnever]

In the first post, I summaraised the changes made in this PR and marked it ready for review.

Hi @jacklovell and @Mateasek, I am awaiting your verdict).

Hello @jrh-ccfe, could you please share any balance.nc file, so I could improve load_solps_from_balance(). I have access to Heimdall cluster, so a path to the file will be enough.

[vsnever]

I updated balance.py and revert the change in the fort44 parser.

[jacklovell]

Thanks for this Vlad: it's a huge amount of great work and a big improvement to the module overall. I've made some in-line comments on the code, but generally it is of a very good standard.

I've run the 3 demo scripts we currently have on a limited set of runs that I have access to and they seem to produce mostly the same results as before. I expect any small differences are a result of bug fixes. Have you done regression testing on any other SOLPS cases you have?

There's a few places in the code with large commented-out blocks: these can be removed if they won't be added back in later. I'm thinking in particular in solps_plasma.py, where there are several commented-out plotting methods which should live in demo scripts instead.

Also, considering the amount of change here it would be worth adding to the changelog while the changes are still fresh. I'll defer to @CnlPepper on whether this is a big enough change to warrant a 1.2 (or even a 2.0) release, or whether a 1.1.1 would be suitable.

[jacklovell]

else is not necessary after a raise. This applies to all the total_radiation*, b_field* and eirene_simulation properties.

[vsnever]

This code was before me, some properties threw errors when not initialized and some did not. I think, they should not throw errors. If a property is required for a function to work, that function must check if the value is set and throw a RuntimeError if not. I've removed these raise statements from properties and instead added them to functions that don't work without these properties.

[jacklovell]

This should be a ConstantVector3D object, not a lambda function.

[jacklovell]

This should be a ConstantVector3D object, not a lambda function.

[jacklovell]

Trailing whitespace.

[jacklovell]

One other comment: if Axisymmetric wrappers of SOLPSFunction2D and SOLPSVectorFunction2D replace the corresponding 3D functions, perhaps we should remove the latter from the module entirely. The 3D functions are still present in this PR.

[CnlPepper]

One other comment: if Axisymmetric wrappers of SOLPSFunction2D and SOLPSVectorFunction2D replace the corresponding 3D functions, perhaps we should remove the latter from the module entirely. The 3D functions are still present in this PR.

I agree, it would be best to remove redundant code.

[CnlPepper]

Also, considering the amount of change here it would be worth adding to the changelog while the changes are still fresh. I'll defer to @CnlPepper on whether this is a big enough change to warrant a 1.2 (or even a 2.0) release, or whether a 1.1.1 would be suitable.

While there are large API changes to the SOLPSSimulation API, the core functionality for loading and generating Plasma objects
and the object layout are essentially the same. It is definitely a new minor release, possibly a major release.

The arguments for a major/minor version increment would be the degree of impact to the main usage of the package. This one lies on the border I suspect. As the active users of this package, you would be best placed to assess the impact on users and then decide if you want to go for v2.0.0 or v1.2.0.

Please do add a changelog. I would recommend getting into the habit of checking if it needs updating with every pull request.

[Mateasek]

shouldn't it be neighbiy?

[Mateasek]

Wouldn't it be better to make this method a standalone instead of class method, like it is in the case of EFITEquilibrium class? This class is becoming quite big, it would help decreasing it. Also the method itself could be useful/used for various purposes. The same could be done with eirene_flux_to_velocity.

[Mateasek]

That is a large amount of work @vsnever, very nice. Thanks for all the work you put into it, especially when SOLPS is not the easiest code to understand and read! I would like to propose one more change, since this is quite major rewrite. I would suggest also removing the commented plotting methods from the SOLPSimulation class. I think the plots would be better done outside in a separate file, as it is with equilibrium plots in Cherab. I don't mean that you should be also adding it now.

[vsnever]

@jacklovell, @Mateasek, @CnlPepper, thank you very much for the code review. I hope I addressed all the issues you raised. For three of them I left the comments (see above). If I did not leave a comment, it means that the issue is fixed exactly as you suggested.

I've run the 3 demo scripts we currently have on a limited set of runs that I have access to and they seem to produce mostly the same results as before. I expect any small differences are a result of bug fixes. Have you done regression testing on any other SOLPS cases you have?

When the data is read from the MDSPlus, the only difference should be in species velocities, because of the bug fix. When the data is read from the balance.nc file, then in addition to velocities, the density of impurity neutral atoms is different, because now it's obtained from the Eirene data. If the data is read from raw simulation files then the density of main plasma neutral atoms is also different because previously it was obtained from B2.

Also, considering the amount of change here it would be worth adding to the changelog while the changes are still fresh. I'll defer to @CnlPepper on whether this is a big enough change to warrant a 1.2 (or even a 2.0) release, or whether a 1.1.1 would be suitable.

I updated the changelog suggesting a 1.2 version.

In addition to the improvements you suggested, I removed redundant check for inside/outside mesh in SOLPSTotalRadiatedPower emitter, and because now this emitter uses only a total_radiation_f3d from SOLPSSimulation, I think that passing the simulation object to this emitter is no longer needed, a Function3D instance is enough.

[jacklovell]

Looks good to me.

[jacklovell]

@CnlPepper @Mateasek are you happy for these changes to be merged?

[Mateasek]

Yes, sorry I am new to this, forgot to click on approve. Very nice piece of work!

[CnlPepper]

Happy here, nice work!

[vsnever]

@jacklovell, @Mateasek, @CnlPepper, thanks a lot!

@jacklovell, please do not close #34 when this is merged. There is still plenty of work to do.