shredEngineer / MagnetiCalc

MagnetiCalc calculates the magnetic field of arbitrary coils.
https://paulwilhelm.de/magneticalc/
Other
45 stars 4 forks source link

Feature request: loading and saving of data #2

Closed CSChisholm closed 2 years ago

CSChisholm commented 2 years ago

I found the user interface very intuitive and the calculation seems efficient on my CPU (I don't have a GPU in my laptop). To enhance usability, I would suggest two changes:

  1. Allow loading a wire by points defined in a text file.
  2. Allow saving the calculated field to a file (I recommend hdf5 for this use case, see below).

Both of these features would allow computation of complex current geometries. The first allowing one to use numpy (or whatever else) to define the coils numerically and then load in the points without having to fill the table by hand. The second would allow one to take the calculated fields of different sets of coils and combine them with the superposition principle and also calculate quantities such as curvature.

Why I think hdf5 is the most suitable file format: First, the h5py library is excellent for reading and writing these file types. The file is structured like a dictionary in Python. Imagine that the evaluation volume is defined by X, Y, Z = numpy.meshgrid(x,y,z) and we have the calculated fields, A and B, which each have three components: A = (A_x, A_y, A_z), B = (B_x, B_y, B_z). We construct a dictionary fields = {'x': x, 'y': y, 'z': z, 'A_x': A_x, 'A_y': A_y, 'A_z': A_z, 'B_x': B_x, 'B_y': B_y, 'B_z': B_z} We can also take the geometry of the wire from Points (here I write xp, yp, zp, I haven't checked the source code for the actual variable names). Let's assume these are the points after stretch and rotate points = {'x': xp, 'y': yp, 'z': zp} data = {'fields': fields, 'points': points, 'current': wirecurrent} To write it in a .h5 file we do

import h5py

def Dict2Group(h5group,data):
    keys = data.keys()
    for key in keys:
        if (isinstance(data[key],dict)):
            group = h5group.create_group(key)
            Dict2Group(group,data[key])
        else:
            h5group[key] = data[key]
    return

f = h5py.File(filename,'w')
Dict2Group(f,data)
f.close()

Then, all of the information which the user could possibly want is saved in one single file for later processing.

shredEngineer commented 2 years ago

@CSChisholm: Thanks for filing this issue, I really appreciate your suggestions! :) Also, I'm glad that you like the UI and performance.

Indeed, the HDF5 format seems to be very well suited – thanks for pointing me to this! (I've worked with netCDF once, which seems similar, but I did not yet know about HDF.)

Exporting the field, wire and current data is straight-forward and I have already implemented the field export today. The only problem is that in the current release either the A-field or the B-field is calculated, but not both. So, for now only the currently selected field type is exported. A possible solution would be to provide a modal dialog which allows the user to select which points and fields to include in the exported file. If both fields are selected, the "missing" field is calculated before the export, while blocking the UI and having saved the other field beforehand. What do you think?

I will implement the wire import from HDF5 tomorrow. (EDIT: Importing from TXT file, okay.) (In the long run, I think I will switch to the HDF5 format completely for storing project files!)

I also plan on releasing version 1.10 including the changes mentioned above this week, and I will notify you about it in this thread.

As an aside, because you mentioned superpositioning multiple geometries / fields: Indeed, that would be very useful. I already plan on making this possible in the UI by using a list of objetcs (wire / current / transformation pipeline) instead of the good old Wire_Widget. :)

Also, thank you for pointing me to "magnetic curvature"; this seems like a metric I could also implement.

Have a good one!

CSChisholm commented 2 years ago

Hi, I looked at the source code thinking I would try to implement the things I requested but I don't have time to get into it for the next few weeks (I also don't have a GPU to test CUDA) so I appreciate that you are working on it so enthusiastically.

I will implement the wire import from HDF5 tomorrow. (EDIT: Importing from TXT file, okay.)

I wrote text file because it's easy but if the output is already going to be hdf5 the input might as well be too.

(In the long run, I think I will switch to the HDF5 format completely for storing project files!)

What I like about the current format (.ini which is actually plain text) is that I could get around using the table to enter points by using a python script to insert points into a template file. This would be irrelevant when loading wires from a file is implemented anyway. I see an advantage in using hdf5 for everything in that it unifies the code and allows more precision for points than the 6 decimal places currently in the .ini files, on the other hand I see that a text file is easier for a user to look at/manipulate as reference for saved results. I suggest that a python module similar to the one I'm attaching at the bottom be included with the code so that when a user installs from pip they would also get a module installed which they could invoke like from magneticalc import hdf5tools for use in python scripts for post-processing or creating wires to load. Once the configuration files are hdf5 we can add functions to process it back to the current format for readability and to process older .ini files into the updated format so no one has to redo all of their work (in fact that one could run in the background if someone tries to load a .ini file after hdf5 is implemented).

One thing I noticed while playing with writing points into the .ini file is that if the loop is not closed (i.e. if the start point and end point are not the same), the program will crash. I guess this was written in with the intention of enforcing Kirchoff's law for current continuity but it's not required for the Biot-Savart law and there are cases where open coils make sense. For example, if I have a coil composed of several turns of wire, the leads probably go away in a semi-random direction as a twisted pair towards a power supply some distance away and one would very rarely need to include the contribution of the leads in the field calculation.

I see a few other tweaks which could be nice but with the hdf5 implementation and if we could implement open wires I think MagnetiCalc will already be an extremely powerful tool. I hope I can find the time to make more tangible contributions soon.

import h5py

def ReadH5(filename):
    '''Opens hdf5 file and converts every group and subgroup into a dictionary
    where the keys are the group keys and the items are the datasets

    return data'''
    #Open file
    f =  h5py.File(filename,'r')
    data = {}
    Helper_Group2Dict(f,data)
    f.close()
    return data

def WriteH5(filename,data):
    '''Takes dictionary and writes hdf5 file using keys as keys and items as
    datasets if they are not dictionaries or groups if they are dictionaries'''
    if not(isinstance(data,dict)):
        raise TypeError('Data must be a dictionary')
    f = h5py.File(filename,'w')
    Helper_Dict2Group(f,data)
    f.close()
    return

def Helper_Group2Dict(h5Group,data):
    '''Helper function for ReadH5'''
    keys = [key for key in h5Group]
    for key in keys:
        if (isinstance(h5Group[key],h5py.Dataset)):
            data[key] = h5Group[key][()]
        else:
            data[key] = {}
            Helper_Group2Dict(h5Group[key],data[key])
    return

def Helper_Dict2Group(h5group,data):
    '''Helper function for WriteH5'''
    keys = data.keys()
    for key in keys:
        if (isinstance(data[key],dict)):
            group = h5group.create_group(key)
            Helper_Dict2Group(group,data[key])
        else:
            h5group[key] = data[key]
    return
shredEngineer commented 2 years ago

I wrote text file because it's easy but if the output is already going to be hdf5 the input might as well be too.

What I like about the current format (.ini which is actually plain text) is that I could get around using the table to enter points by using a python script to insert points into a template file. This would be irrelevant when loading wires from a file is implemented anyway. I see an advantage in using hdf5 for everything in that it unifies the code and allows more precision for points than the 6 decimal places currently in the .ini files, on the other hand I see that a text file is easier for a user to look at/manipulate as reference for saved results.

Right, both formats have their advantages and disadvantages. I guess for now the TXT import/export for wires is a good step in the right direction, and the HDF5 container export opens many possibilities for post-processing. I think the complete switch-over to HDF5 will happen in a 2.0 release in the future, as it implies some bigger changes to the Config class, while needing to retain the functionality to import old INI files.

I suggest that a python module similar to the one I'm attaching at the bottom be included with the code so that when a user installs from pip they would also get a module installed which they could invoke like from magneticalc import hdf5tools for use in python scripts for post-processing or creating wires to load. Once the configuration files are hdf5 we can add functions to process it back to the current format for readability and to process older .ini files into the updated format so no one has to redo all of their work (in fact that one could run in the background if someone tries to load a .ini file after hdf5 is implemented).

Thank you for the suggestion! An easy-to-use interface like that seems very handy, and I will look into how it would best be structured! :)

One thing I noticed while playing with writing points into the .ini file is that if the loop is not closed (i.e. if the start point and end point are not the same), the program will crash. I guess this was written in with the intention of enforcing Kirchoff's law for current continuity but it's not required for the Biot-Savart law and there are cases where open coils make sense.

Actually, there is no check that the loop is closed; in fact, MagnetiCalc just automatically closes the loop, no matter what. May I guess that you left a stray semicolon at the end of the line in the .ini file? I just checked, and indeed this does crash the program. I have now inserted a check to ignore this kind of error.

For example, if I have a coil composed of several turns of wire, the leads probably go away in a semi-random direction as a twisted pair towards a power supply some distance away and one would very rarely need to include the contribution of the leads in the field calculation.

Thank you for pointing this out! I have now added a checkbox that allows the user to disable "loop closure".

I see a few other tweaks which could be nice but with the hdf5 implementation and if we could implement open wires I think MagnetiCalc will already be an extremely powerful tool. I hope I can find the time to make more tangible contributions soon.

That sounds great. :) Feel free to open issues and/or contribute to the code at any time. Any feedback is warmly welcome, and you seem like someone more knowledgable in this field than myself (not actually being a physicist).

CSChisholm commented 2 years ago

May I guess that you left a stray semicolon at the end of the line in the .ini file? I just checked, and indeed this does crash the program. I have now inserted a check to ignore this kind of error.

Indeed, at the beginning I had this problem but I did manage to get this method to work with a closed loop. Actually, the crash is not on opening but on calculation. I found the minimal example so I'll write it up in a new thread.

shredEngineer commented 2 years ago

Wire import/export from/to TXT file and field/wire/current export to HDF5 file is now implemented in 28d43f9fa365e80a67e78d79947b947d23a443b7.

(Oops, clicked "Close issue" one too many times. :sweat_smile:)

shredEngineer commented 2 years ago

@CSChisholm I just added an API class based on your suggestion above: 6999b95 (Thank you again for the code, which I modified just slightly.)

I also added a corresponding section to the README. There I include the basic code for importing the HDF5 file. I want to extend the code to also plot the field using Matplotlib.

But there I have a problem. I'm too stupid to plot the data right now … I get confused over the array shapes and the ordering (reshape, meshgrid, ??????) after having implemented your pull request #5. :)

Could you help me out by any chance?!

This is the basic construct needed for plotting, where X, Y, Z are 1D lists of coordinates and U, V, W are the vector components at these coordinates:

ax = plt.figure().add_subplot(projection="3d")
ax.quiver(X, Y, Z, U, V, W, length=.1, normalize=True)
plt.show()
CSChisholm commented 2 years ago

My pull request #5 was to put the arrays into a shape that felt more natural as someone who does data analysis with numpy (and Matlab a long time ago). I have to say it's a personal preference and this is your project so there's no reason you can't remove my changes.

However, if you want to keep the changes this code should reshape your arrays the way you want

data = {} #Dictionary constructed recursively from .hdf5 file
x = data['fields']['x']
y = data['fields']['y']
z = data['fields']['z']
U = data['fields']['B_x'].ravel()
V = data['fields']['B_y'].ravel()
W = data['fields']['B_z'].ravel()
X, Y, Z = np.meshgrid(x,y,z,indexing='ij')
X = X.ravel()
Y = Y.ravel()
Z = Z.ravel()

ax = plt.figure().add_subplot(projection="3d")
ax.quiver(X, Y, Z, U, V, W, length=.1, normalize=True)
plt.show()
CSChisholm commented 2 years ago

Actually it was me who got the indexing wrong for #5. I suggest removing #5 and I'll work on a function for the API which will create the arrays #5 was supposed to make if they're wanted.

shredEngineer commented 2 years ago

Actually it was me who got the indexing wrong for #5. I suggest removing #5 and I'll work on a function for the API which will create the arrays #5 was supposed to make if they're wanted.

OK, that sounds good! To have the exported data in a format "native" to the application does seem to be the better choice. Then any conversion can happen later when using the API.

Looking forward! :)

CSChisholm commented 2 years ago

According to the NumPy documentation the whole problem was that I forgot that the default behaviour of np.reshape is 'C-like' put I should have put order='F' for 'Fortran-like' behaviour as indicated by SamplingVolume.Recalculate.i_to_xyz(). So I'll be able to add this function very soon.

shredEngineer commented 2 years ago

Hi Craig,

I changed the API one more time, and I hope you will find it useful! :) (And I hope this doesn't break too much for you…) You could take a look at the corresponding README section which I updated.

The API now returns a MagnetiCalc_Data object, which can be accessed as a dictionary. Moreover it provides some handy member functions for accessing the data and reshaping it.

This means I moved your reshape_data() function from API to MagnetiCalc_Data and split it up into smaller ones. I think this is the best solution as the new convenience functions are more granular and easier to use.

Also, I now include the sampling volume dimension in the exported HDF5 container so the fields can be unraveled without having to reduce the axes first just to find out their dimensions.

CSChisholm commented 2 years ago

Hi Paul,

I took a look at the changes and I like this implementation. I just made a change to MagnetiCalc_Data.get_field() because you forgot to reshape field_y and field_z when as_3d = True #7

shredEngineer commented 2 years ago

Oops, that`s what you get for pulling an all-nighter. Thanks! :)