Pixel Current Artefact

[Kyroba]

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)

[loganbvh]

Hi @Kyroba,

I think this type of artifact happens when the time step is slightly too big, making the numerics slightly unstable. I would try re-running the simulation with options.dt_max set to, for example, 80% of the value used for the simulation shown above. I think it is mostly a cosmetic thing - it shouldn't significantly affect simulation results, e.g. voltage.

[Kyroba]

I varied dt_max from default 0.1 to 0.05, 0.01, 0.001 yet the artefact was still there. Didn't seem to change in size or position at all either.

[loganbvh]

Can you try increasing the smooth parameter in device.make_mesh()?

[loganbvh]

Also, can you paste the results of solution.dynamics.plot_dt()? Thanks

[Kyroba]

I changed the smooth parameter from 100 to 1000 but still the same result:

Here is the output from solution.dynamics.plot_dt():

[loganbvh]

Can you try either options.dt_max = 5e-4 or resampling the film and hole boundary to have ~30% more points?

The basic problem is that the current density is defined on the edges of the mesh. For plotting, I need to sort of interpolate the current density to evaluate it at the mesh sites rather than edges. This is done by averaging the current density over each edge that's connected to a given mesh site, see below. In certain cases (I think when one or more of the triangles incident to a site is very long/narrow), this average produces artifacts.

py-tdgl/tdgl/finite_volume/mesh.py

Lines 203 to 243 in b779261

	def get_quantity_on_site(
	self,
	quantity_on_edge: np.ndarray,
	vector: bool = True,
	use_cupy: bool = False,
	) -> np.ndarray:
	"""Compute the quantity on site by averaging over all edges
	connecting to each site.
	Args:
	quantity_on_edge: Observable on the edges.
	vector: Whether ``quantity_on_edge`` is a vector quantity.
	use_cupy. Whether to use ``cupy`` interface.
	Returns:
	The quantity vector or scalar at each site.
	"""
	normalized_directions = self.edge_mesh.normalized_directions
	edges = self.edge_mesh.edges
	if use_cupy:
	xp = cupy
	normalized_directions = xp.asarray(normalized_directions)
	edges = xp.asarray(edges)
	else:
	xp = np
	if vector:
	flux_x = quantity_on_edge * normalized_directions[:, 0]
	flux_y = quantity_on_edge * normalized_directions[:, 1]
	else:
	flux_x = flux_y = quantity_on_edge
	# Sum x and y components for every edge connecting to the vertex
	vertices = xp.concatenate([edges[:, 0], edges[:, 1]])
	x_values = xp.concatenate([flux_x, flux_x])
	y_values = xp.concatenate([flux_y, flux_y])
	counts = xp.bincount(vertices)
	x_group_values = xp.bincount(vertices, weights=x_values) / counts
	y_group_values = xp.bincount(vertices, weights=y_values) / counts
	vector_val = xp.array([x_group_values, y_group_values]).T / 2
	if vector:
	return vector_val
	return vector_val[:, 0]

[Kyroba]

That dt_max value did the trick! Your explanation also makes a lot of sense. Thank you!

[loganbvh]

So for a rule of thumb, I would say that if you see this type of artifact, look at solution.dynamics.plot_dt() and try setting options.dt_max somewhere near the median of the distribution of dt values.