loganbvh / py-tdgl Goto Github PK

https://docs.cupy.dev/en/stable/reference/generated/cupyx.scipy.sparse.linalg.factorized.html

The LU factorization is done on the CPU, but the linear solve is done on the GPU.

I vaguely remember something about loading the hd5 files into ParaView for visualizing in 3D but can't seem to find any mention of it anywhere. I loaded the file into ParaView and of the various readers that were available only 1 seemed to work however it did not seem to import the mesh appropriately. I used this reader:

For a large geometry (1um square with 3nm coherence length, 6nm hole pattern), the hole fluxoid error is high - suggesting that screening is necessary. I chose parameters such that effective magnetic screening length = 90 microns.

Is this an issue with how I am calculating the fluxoid or do I really need screening for these simulations. Below are images of the device and order parameter at a specific time (with bias).

For long-running simulations, it would be nice to have a separate thread/process that plots the current state of the system in real time. I think this can be implemented safely using https://docs.h5py.org/en/stable/swmr.html

I think this issue may have risen after the snapshot feature was implemented. The command I use to produce animations no longer works as the --output argument is not recognized. Removing the output argument does not work either as it does for snapshot since it fails to produce an animation without a file to dump to.

python -m tdgl.visualize --input "nobiasExp3.h5" animate --output "anim.mp4" --quantities ORDER_PARAMETER

The plotting of data whilst the simulation is running is extremely taxing on my pc. The figures do not even load when simulating geometries 300*xi. This could be due to the density of data and having to output basically 4 different figures, I'm not exactly sure. Perhaps a subsample of the data could be beneficial to improve efficiency? I have not done much testing on this, so it might be working on other machines whilst mine has an error somewhere.

First of all, what a wonderfully simple to use package. Thank you very much for designing this.

The issue I am experiencing is that the results are slightly different when running the exact same Quickstart notebook to what is shown in (https://py-tdgl.readthedocs.io/en/latest/notebooks/quickstart.html). I believe this might be a result of the mesh generation, was there a change in how it is done? I noticed that the tutorial uses tdgl version 0.4.0 whereas I am using the latest release 0.3.1.

What hardware was used to achieve those computation times in the quickstart tutorial?

I have attached my notebook results for more information if needed.
PC specs: i9 9900k processor, 2080 ti gpu, 32 gb ram 3200Mhz, 980 PRO NVME SSD.

JupyterLab_Quickstart.pdf

I was wondering if it is possible to introduce a thickness gradient to the device. For instance, thickness of the thin film is largest at the edges and decreases towards the center, similar to a bowl shape. I understand the thickness value in tdgl.device requests a float. Would this kind of change in thickness even be realistic considering how the TDGL is done?

I was also wondering if it is possible to combine several tdgl.layer together to make a Josephson junction for example (in the same plane or even vertically stacked).

Perhaps I could spatially vary disorder_epsilon in tdgl.solve to mimic thickness changes but I am not sure how I would go about doing this.

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

I would like to study the dynamics of the vortex lattice as the field is increased as a step function. This can easily be done if I define a time dependent field such as

def step_function(x, y, z, *, t, tmax, step_size, final = 1.2):
    if t < tmax:
        return final * np.trunc(t / step_size) * step_size / (tmax - step_size)
    return final

tmax = 3000
final = 1.5
step_size=250
t_dependence = tdgl.Parameter(step_function, tmax=tmax, step_size=step_size, final=final, time_dependent=True)
increasing_vector_potential = tdgl.sources.ConstantField(1., field_units=field_units, length_units=length_units) * t_dependence
options = tdgl.SolverOptions(
    skip_time=0,
    solve_time=3000,
    save_every=200,
    field_units=field_units,
    output_file=os.path.join(tempdir.name, "increasing-field.h5"),
)
solution = tdgl.solve(device, options, decreasing_vector_potential)

The problem is that having a time dependent field greatly increases the runtime. My assumption is that at each step the spatial link variables are recomputed. Since the field remains constant within each step (in terms of the increase of the filed, not simulation time steps), I was able to work my way around this issue by solving the equations one step at a time, using the previous solution as the starting point at each new step:

step_size = 250
tmax = 3000
steps = int(np.ceil(tmax / step_size))
final = 1.5
options = tdgl.SolverOptions(
    skip_time=0,
    solve_time=step_size,
    save_every=200,
    field_units=field_units,
    output_file=os.path.join(tempdir.name, "increasing-field-steps.h5"),
)
for step in range(steps):
    b_z = (final / (steps - 1)) * step
    field = tdgl.sources.ConstantField(b_z, field_units=field_units, length_units=length_units)
    if step == 0:
        solution = tdgl.solve(device, options, applied_vector_potential=b_z)
    else:
        solution = tdgl.solve(device, options, applied_vector_potential=b_z, seed_solution=solution)

By doing this, I decreased the runtime from 63 minutes to only 8! The only problem, though, is that if I want to create an animation, it will only do so for the last step of the solution rather than the solution as a whole. I wanted to know if it would be possible to have an option to append the new solution to the previous one for this matter.

Thank you for your time and attention and thank you for creating this package!

Can the current_through_path be calculated while solving and plotted similarly to DynamicsData.plot()?

I'm new to python and tdgl. After installing the software I ran the tdgl.testing program and it generated the following errors.

FAILED tdgl/test/test_solve.py::test_source_drain_current[1-True-True-True-0-5.0] - AssertionError
FAILED tdgl/test/test_solve.py::test_source_drain_current[1-True-True-True-0-] - AssertionError
FAILED tdgl/test/test_solve.py::test_source_drain_current[1-True-True-True-1-5.0] - AssertionError
FAILED tdgl/test/test_solve.py::test_source_drain_current[1-True-True-True-1-] - AssertionError

In scrolling through the code that is printed with the error I found this repeated every time.

      assert cupy is not None

E AssertionError

However, I did install cupy

python -m pip install -U setuptools pip
pip install cupy-cuda12x
pip install cuTENSOR.

but was unable to install the additional libraries like NCCL, etc
However, I'm not sure if that is really the problem. There are instructions to set tdgl.SolverOptions = True, but I'm not sure how to do that. I would assume it would default to the CPU if I didn't set it to TRUE?

Can anyone confirm that not installing those libraries is the problem? Anybody run into this problem and know how to fix it?

I followed the instructions for installation (installing through PyPI) and ran through the quickstart. I ran the testing suite and all of the tests passed. However, I noticed a TypeError when I tried to change the solver type to use PyPardiso (which I installed using conda install -c conda-forge pypardiso):

Here's the stack trace:

Traceback (most recent call last):
  File "[.../lib/python3.10/site-packages/tdgl/solver/solver.py", line 799](http://localhost:8888/.../lib/python3.10/site-packages/tdgl/solver/solver.py#line=798), in solve
    data_was_generated = runner.run()
  File "[.../lib/python3.10/site-packages/tdgl/solver/runner.py", line 305](http://localhost:8888/.../lib/python3.10/site-packages/tdgl/solver/runner.py#line=304), in run
    success = self._run_stage(
  File "[.../lib/python3.10/site-packages/tdgl/solver/runner.py", line 416](http://localhost:8888/.../lib/python3.10/site-packages/tdgl/solver/runner.py#line=415), in _run_stage
    function_result = self.function(
  File "[.../lib/python3.10/site-packages/tdgl/solver/solver.py", line 675](http://localhost:8888/.../lib/python3.10/site-packages/tdgl/solver/solver.py#line=674), in update
    mu, supercurrent, normal_current = self.solve_for_observables(psi, dA_dt)
  File "[.../lib/python3.10/site-packages/tdgl/solver/solver.py", line 509](http://localhost:8888/.../lib/python3.10/site-packages/tdgl/solver/solver.py#line=508), in solve_for_observables
    mu = pypardiso.spsolve(operators.mu_laplacian, rhs)
  File "[.../lib/python3.10/site-packages/pypardiso/scipy_aliases.py", line 44](http://localhost:8888/.../lib/python3.10/site-packages/pypardiso/scipy_aliases.py#line=43), in spsolve
    solver._check_A(A)
  File "[.../lib/python3.10/site-packages/pypardiso/pardiso_wrapper.py", line 227](http://localhost:8888/.../lib/python3.10/site-packages/pypardiso/pardiso_wrapper.py#line=226), in _check_A
    raise TypeError(msg)
TypeError: PyPardiso requires matrix A to be in CSR or CSC format, but matrix A is: <class 'scipy.sparse._csc.csc_array'>

It looks like when mu_laplacian is generated, it gets generated as a csc_array, but PyPardiso requires a csc_matrix.

If the following section is modified:

if self.sparse_solver is SparseSolver.CUPY:
...
elif self.sparse_solver is SparseSolver.PARDISO:
+ self.mu_laplacian = sp._csc.csc_matrix(self.mu_laplacian)
  self.mu_laplacian_lu = None

Then the simulation runs fine (although it's not any faster than using the default SparseLU solver, but perhaps that's just because of the structure of the example simulation geometry).

My understanding of superconducting vortices is that they nucleate from edges and defects in a sample in the presence of a magnetic field. This translates well in py-tdgl when looking at small geometries. However, I was investigating a large shape with dimensions roughly 40x100 um. In the initialization of vortices using high field, it appears that the vortices are generated all over the sample instead of just from the edges. Is this due to imprecision in numerical handling causing 'defects' to appear or something else? Below is the output animation for zero applied current (non-zero vector potential and default epsilon):

Full Sim

Sim Start Slowed

Could this be the result of a rapid change in applied field?

When attempting a nanoSQUID simulation, I noticed a strange current spike stuck in a specific location of the loop. This current spike does not move or react to applied fields. You can see it just below the left weak-link (25mT applied field no bias):

The first thing I went to check was the mesh to see if there is some kind of dislocation:

But it doesn't seem like anything is wrong. Changing mesh parameters moves the artefact around, but I can never get rid of it. I have also iteratively removed each polygon to see if something was causing the issue, and to my surprise, I saw this artefact in a plain square film with no holes. The only way I was able to get rid of it was to use entirely different xi and lambda values.

length_units = "nm"
# Material parameters
xi = 25
london_lambda = 100
d = 5
layer = tdgl.Layer(coherence_length=xi, london_lambda=london_lambda, thickness=d, gamma=10)

Is it possible to have an option to output the 'total magnetic moment' of the thin film using the number of vortices? The end goal would be to plot moment against applied field for direct comparison with experimental results.

I think having moment vs field, Jc vs field, and IV curves as easily usable functions would improve prototyping.

In "pyTDGL: Time-dependent Ginzburg-Landau in Python" there is a section on the nanoSQUID figure 5. A current-voltage-flux plot is provided but there is no clear way of how this was done in the documents. Would you need to have run hundreds of simulations with varying fields in order to obtain that plot? That would result in an extremely large file, im hoping there is a more elegant solution.

Running quickstart.ipynb on Colab results in the error:

AssertionError                            Traceback (most recent call last)

[<ipython-input-24-3e6395c7c2a1>](https://localhost:8080/#) in <cell line: 11>()
      9 )
     10 # If you do not provide an applied_vector_potential, tdgl defaults to zero applied field.
---> 11 zero_field_solution = tdgl.solve(
     12     device,
     13     options,

1 frames

[/usr/local/lib/python3.10/dist-packages/tdgl/solver/solver.py](https://localhost:8080/#) in __init__(self, device, options, applied_vector_potential, terminal_currents, disorder_epsilon, seed_solution)
    250 
    251         J_scale = 4 * ((ureg(current_units) / length_units) / K0).to_base_units()
--> 252         assert "dimensionless" in str(J_scale.units), str(J_scale.units)
    253         J_scale = J_scale.magnitude
    254         self.current_func = lambda t: {

AssertionError:

The error also occurred on two separate devices with pint freshly installed. Clearly the error has something to do with units, so I tried installing the previous version (0.23), which seemed to fix the problem.

First off, thanks so much for creating this!! It is exactly the kind of simulation tool that the field needs.

Is it possible to have a progress bar added for the mesh creation? I have been trying to set up some simulations of thin film NbN (xi = 0.005, london_lambda = .200, d = .010) with a total_width=1 and total_length=2 [µm]. For those dimensions and a max_edge_length=xi/2, the mesh is quite large and it is unclear if the mesh will be generated or if there is some "freezing" error. Additionally, there is a scenario where the simulation seems to run (thermalization and simulation) when creating the mesh.

I have attached two console logs of the simulation error freezing error.

ipython_mesh_error2.txt
ipython_mesh_error3.txt

In previous issues, it was suggested that the Adapative library (https://github.com/python-adaptive/adaptive) could be used to sample the most important parameters when finding Jc across various fields. I have been thinking of applying oscillating fields (without changing polarity) to a superconducting film and was wondering what considerations to take when using Adapative to scan the through the frequency range. What could be the output that the algorithm learns from? I was thinking average vortex motion, but not sure how to implement this. Maybe the peak to peak voltage between probes is easier?

I was wondering if anyone had any details on the simulation parameters for YBCO. I understand almost all of them must be derived experimentally but I was just looking for a few estimates. Coherence length and penetration depth are easily available but gamma is a lot more difficult to find if at all. I have found a few papers where they provide their value of gamma but it is the ratio of anisotropies and not the same gamma as in PyTDGL.

I found a paper which had a phonon scattering time t_in approximately 20ps and a superconducting gap del_0 of 20meV which provides a gamma of approximately 1210. I found another resource with GL relaxation time t_0 of 0.03 ps which I believe relates to the scattering time. Although, I am not sure how reliable this information is as I cannot access the paper it is using.

This yields a gamma of roughly 0.29

So, I am left qualitatively watching the dynamics of the simulation to see which one closely resembles how vortices behave as I apply a current in an YBCO sample. I don't have the equipment at the moment to experimentally verify it myself.

I was wondering if there could be a method to investigate time dependent external fields?

You could potentially solve the system for 1 field, then use the final state as the seed for the next field. However, this can quickly become tedious. Is there a simpler solution available to keep it all in 1 file?

Thank you!

This was touched on before in several other issues but I thought I would bring it up again in case any development has happened. It is difficult to generate a mesh for extremely large geometries (> 1million mesh sites) since the mesh malforms if the film has not been sampled correctly. This is basically an arbitrary value which sometimes I can get right on the first several tries, and other times I am spending hours trying to obtain a good value. I am currently attempting a method of just repeatedly changing the resample to whatever the first malformed mesh site is at and seeing if this works.

Any ideas? I have tried just increasing resamples by 50% but this does not work for the really large geometries that I am simulating. I also want to not go overboard with the resample number since it can also drastically increase the number of mesh sites unnecessarily.

First of all, thank you very much for developing this package!

I am new to simulating SC devices, and one area of my interest is in SIS / SNS junctions.

For this, is it possible to add functionality to make a device layer with different material instead of a constriction, to serve as an insulator or normal metal interface at the junction?

Thank you!

is it possible to have an option for a time-dependent disorder_epsilon callable function? In the same way current_func(t) is called?

something like disorder_epsilon(r, t)

I was attempting to make an array of holes in the thin film. I finally managed to do it and here is the output:

I essentially made an empty list and appended holes to the list using a loop.

The issue arises right at the end of the simulation. Here is the code:

from tdgl.sources import ConstantField, LinearRamp

options = tdgl.SolverOptions(
    gpu = True,
    solve_time=400,
    sparse_solver="superlu",
    output_file=os.path.join("Z:/Archive/TDGL Sims/Paper 1", "nobiasTriUP.h5"),
    field_units = "mT",
    current_units="uA",
)

uniform = ConstantField(
    0.6,  
    field_units="mT",
    length_units="um",   # Set your length unit, default is "um"
)

A_applied = (LinearRamp(tmin=0, tmax=200) * uniform)

zero_current_solution = tdgl.solve(
    device,
    options,
    # If applied_vector_potential is given as a single number,
    # it is interpreted to mean the vector potential associated with a
    # uniform out-of-plane magnetic field with the specified strength.
    applied_vector_potential = A_applied,
)

and here is the output:

Simulating: 100%|█████████████████████████████████████████████████████████████████▉| 400/400 [39:31<00:00,  5.93s/tau ]
Ignoring the following exception in DataHandler.__exit__():
Traceback (most recent call last):
  File "C:\Users\Ca\AppData\Local\Programs\Python\Python310\lib\site-packages\tdgl\solver\solve.py", line 491, in solve
    solution.to_hdf5()
  File "C:\Users\Ca\AppData\Local\Programs\Python\Python310\lib\site-packages\tdgl\solution\solution.py", line 914, in to_hdf5
    self._save_to_hdf5_file(self.path, save_mesh=save_mesh)
  File "C:\Users\Ca\AppData\Local\Programs\Python\Python310\lib\site-packages\tdgl\solution\solution.py", line 901, in _save_to_hdf5_file
    self.device.to_hdf5(group.create_group("device"), save_mesh=save_mesh)
  File "C:\Users\Ca\AppData\Local\Programs\Python\Python310\lib\site-packages\tdgl\device\device.py", line 805, in to_hdf5
    hole.to_hdf5(group.create_group(hole.name))
  File "C:\Users\Ca\AppData\Local\Programs\Python\Python310\lib\site-packages\h5py\_hl\group.py", line 64, in create_group
    gid = h5g.create(self.id, name, lcpl=lcpl, gcpl=gcpl)
  File "h5py\_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
  File "h5py\_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
  File "h5py\h5g.pyx", line 166, in h5py.h5g.create
ValueError: Unable to create group (name already exists)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[40], line 20
     12 uniform = ConstantField(
     13     0.6,  
     14     field_units="mT",
     15     length_units="um",   # Set your length unit, default is "um"
     16 )
     18 A_applied = (LinearRamp(tmin=0, tmax=200) * uniform)
---> 20 zero_current_solution = tdgl.solve(
     21     device,
     22     options,
     23     # If applied_vector_potential is given as a single number,
     24     # it is interpreted to mean the vector potential associated with a
     25     # uniform out-of-plane magnetic field with the specified strength.
     26     applied_vector_potential = A_applied,
     27 )

File ~\AppData\Local\Programs\Python\Python310\lib\site-packages\tdgl\solver\solve.py:491, in solve(device, options, applied_vector_potential, terminal_currents, disorder_epsilon, seed_solution)
    481 if result:
    482     solution = Solution(
    483         device=device,
    484         path=data_handler.output_path,
   (...)
    489         total_seconds=(end_time - start_time).total_seconds(),
    490     )
--> 491     solution.to_hdf5()
    492     return solution
    493 return None

File ~\AppData\Local\Programs\Python\Python310\lib\site-packages\tdgl\solution\solution.py:914, in Solution.to_hdf5(self, h5path, save_mesh)
    912 if self.saved_on_disk:
    913     if h5path is None:
--> 914         self._save_to_hdf5_file(self.path, save_mesh=save_mesh)
    915     else:
    916         shutil.copy(self.path, h5path)

File ~\AppData\Local\Programs\Python\Python310\lib\site-packages\tdgl\solution\solution.py:901, in Solution._save_to_hdf5_file(self, h5file, save_tdgl_data, save_mesh)
    895 serialize_func(
    896     self.disorder_epsilon,
    897     "disorder_epsilon",
    898     group,
    899 )
    900 group.attrs["total_seconds"] = self.total_seconds
--> 901 self.device.to_hdf5(group.create_group("device"), save_mesh=save_mesh)

File ~\AppData\Local\Programs\Python\Python310\lib\site-packages\tdgl\device\device.py:805, in Device.to_hdf5(self, path_or_group, save_mesh)
    803 for hole in sorted(self.holes, key=attrgetter("name")):
    804     group = f.require_group("holes")
--> 805     hole.to_hdf5(group.create_group(hole.name))
    806 if save_mesh and self.mesh is not None:
    807     self.mesh.to_hdf5(f.create_group("mesh"))

File ~\AppData\Local\Programs\Python\Python310\lib\site-packages\h5py\_hl\group.py:64, in Group.create_group(self, name, track_order)
     62 name, lcpl = self._e(name, lcpl=True)
     63 gcpl = Group._gcpl_crt_order if track_order else None
---> 64 gid = h5g.create(self.id, name, lcpl=lcpl, gcpl=gcpl)
     65 return Group(gid)

File h5py\_objects.pyx:54, in h5py._objects.with_phil.wrapper()

File h5py\_objects.pyx:55, in h5py._objects.with_phil.wrapper()

File h5py\h5g.pyx:166, in h5py.h5g.create()

ValueError: Unable to create group (name already exists)

The XY axis scales are incorrect and do not accurately portray the real size of the geometry. Defining the geometry:

# Device geometry
total_width = 100
total_length = 40
ch_width = 5
ch_length = 20

channel1 = (
    tdgl.Polygon(points=box(ch_width, ch_length))
    .translate(dy=-total_length/2-ch_length/2)
    .translate(dx=-total_width/4)
)

channel2 = channel1.scale(xfact=-1)

channel3 = (
    tdgl.Polygon(points=box(ch_width, ch_length))
    .translate(dy=total_length/2 + ch_length/2)
)


film = (
    tdgl.Polygon("film", points=box(total_width, total_length))
    .union(channel1, channel2, channel3)
    .resample(5001)
    .buffer(0)
)

The following image is displayed accurately:

And this is the output from using CLI (python -m tdgl.visualize --input "weak-link-zero-current.h5" animate --output "anim.mp4" --quantities ORDER_PARAMETER):

Is it possible for probe points to be either moved/added after the simulation is done? I would assume this would only require the recalculation of the electric potential over the points at each solve step?

When dealing with large geometries, small holes can sometimes not be visible in the final animation. Zooming in by changing the axis limits would be a nice work around or even a way to highlight holes in the final image?

Here you can see it initializes with the effect from the tiny holes, but then quickly becomes invisible.

Is it possible in this particular realization of the TDGL to vary the temperature from above Tc to below Tc in the presence of a magnetic field so as to field-cool the superconductor? I would like to trap magnetic flux inside a hole in my superconducting device.