loganbvh / py-tdgl

2D time-dependent Ginzburg-Landau in Python
MIT License
37 stars 13 forks source link

Traceback on quickstart.ipynb #82

Closed 5hwetabh closed 5 months ago

5hwetabh commented 5 months ago

Firstly, thanks for creating this incredibly helpful package!

I have tried running the 'quickstart.ipynb' notebook as-is (except setting MAKE_ANIMATIONS = 0), but I keep running into the same traceback across both my Windows laptops and various versions of python: On code cell 21, the line fluxoid = zero_current_solution.polygon_fluxoid(polygon, with_units=False) results in an index error:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[21], line 3
      1 for name, (radius, center) in fluxoid_polygons.items():
      2     polygon = circle(radius, center=center, points=201)
----> 3     fluxoid = zero_current_solution.polygon_fluxoid(polygon, with_units=False)
      4     print(
      5         f"{name}:\n\t{fluxoid} Phi_0\n\tTotal fluxoid: {sum(fluxoid):.2f} Phi_0\n"
      6     )

File \.conda\Lib\site-packages\tdgl\solution\solution.py:526, in Solution.polygon_fluxoid(self, polygon_points, interp_method, units, with_units)
    524 dl = np.diff(points, axis=0, prepend=points[:1]) * ureg(device.length_units)
    525 A_units = f"{self.field_units} * {device.length_units}"
--> 526 A_poly = self.vector_potential_at_position(
    527     points,
    528     zs=zs,
    529     units=A_units,
    530     with_units=True,
    531     return_sum=True,
    532 )[:, :2]
    533 # Compute the flux part of the fluxoid:
    534 # \oint_{\\partial poly} \vec{A}\cdot\mathrm{d}\vec{r}
    535 int_A = np.trapz((A_poly * dl).sum(axis=1))

File \.conda\Lib\site-packages\tdgl\solution\solution.py:837, in Solution.vector_potential_at_position(self, positions, zs, units, with_units, return_sum)
    835 A_kwargs = {}
    836 if self.applied_vector_potential.time_dependent:
--> 837     A_kwargs["t"] = self.times[self.solve_step]
    838 applied = self.applied_vector_potential(
    839     positions[:, 0],
    840     positions[:, 1],
    841     zs.squeeze(),
    842     **A_kwargs,
    843 )
    844 if applied.shape[1] == 2:

IndexError: index 308 is out of bounds for axis 0 with size 308

I am not able to find the source of the issue directly, but I did check some relevant values which hopefully will help:

len(zero_current_solution.dynamics.time) returns 30701 len(zero_current_solution.times) returns 308 len(h5py.File(zero_current_solution.path, 'r')["data"]) returns 309

loganbvh commented 5 months ago

Hi, thanks for reporting this. Can you share the output of tdgl.version_dict()?

A fresh install on Google Colab executes the entire quickstart.ipynb notebook successfully:

image
5hwetabh commented 5 months ago

This is the output on my .conda env:

{'tdgl': '0.8.2',
 'Numpy': '1.26.4',
 'SciPy': '1.10.1',
 'matplotlib': '3.9.0',
 'cupy': 'None',
 'numba': '0.60.0',
 'IPython': '8.25.0',
 'Python': '3.11.9 | packaged by Anaconda, Inc. | (main, Apr 19 2024, 16:40:41) [MSC v.1916 64 bit (AMD64)]',
 'OS': 'nt [win32]',
 'Number of CPUs': 'Physical: 8, Logical: 16',
 'BLAS Info': 'Generic'}

And the output on my .venv:

{'tdgl': '0.8.2',
 'Numpy': '1.25.2',
 'SciPy': '1.10.1',
 'matplotlib': '3.7.1',
 'cupy': 'None',
 'numba': '0.58.1',
 'IPython': '8.25.0',
 'Python': '3.11.9 (tags/v3.11.9:de54cf5, Apr  2 2024, 10:12:12) [MSC v.1938 64 bit (AMD64)]',
 'OS': 'nt [win32]',
 'Number of CPUs': 'Physical: 8, Logical: 16',
 'BLAS Info': 'OPENBLAS'}

The same traceback occurs in both.

I think I found the source of the traceback in solution/solution.py: image if len(Solution.dynamics.times)%step == 1 (as on my laptop) the resulting Solution.times array is one shorter than it would otherwise be, which might be what causes this traceback.

loganbvh commented 5 months ago

Does't line 148 in that screenshot handle the case len(solution.dynamics.times)%step == 1?

Can you print out solution.dynamics.time[-3:] and solution.times[-3:]?

loganbvh commented 5 months ago

Actually can you just upload the HDF5 file here so I can take a closer look? Thanks

5hwetabh commented 5 months ago

It seems GitHub doesn't allow the large file, so I shared it via DropBox instead. I have never used h5py before so I hope it is the correct file.

https://www.dropbox.com/scl/fi/wr6ey3jx12nwlw73gjtg5/weak-link-zero-current.h5?rlkey=p5u3qw6oobz62f5k71n8xjcdg&st=asgisgxv&dl=0

I am under the impression lines 145-146 handle cases where len(Solution.dynamics.times)%step == 1 and line 148 handles it when false (i.e. 99% of the time); I apologize if I am incorrect.

Thanks for helping me with this issue!

loganbvh commented 5 months ago

Thanks. I think the issue is that if the total number of solve steps ends up being a multiple of options.save_every, the data for the final time step is saved twice. The data for solve step 30700 is saved in both data[307] and data[308] in the HDF5 file (see below). It is very rare that a simulation with an adaptive time step will end on a step that is a multiple of options.save_every, which I think is why I have never seen this problem before.

I opened a Pull Request that I think will fix the issue: https://github.com/loganbvh/py-tdgl/pull/83. I will merge the Pull Request and release a new version on PyPI either later today or tomorrow. It is a one line fix: https://github.com/loganbvh/py-tdgl/pull/83/commits/1cbeefb2ea200a9294095ea3a91c1e0523a560c4

HDF5 details image image image
5hwetabh commented 5 months ago

Thanks a ton for the prompt fix! I am happy to say that I have learned a lot of Python from this issue.

Edit: I can confirm I no longer get a traceback after manually applying the fix.