Initial fix for terminal currents

superscreen/device/mesh.py

@@ -9,7 +9,7 @@
 from matplotlib.tri import Triangulation
 
 from ..distance import q_matrix
-from ..fem import gradient_vertices, laplace_operator
+from ..fem import gradient_triangles, gradient_vertices, laplace_operator
 from . import utils
 from .edge_mesh import EdgeMesh
 
@@ -27,6 +27,7 @@ class Mesh:
  form a triangle. [[0, 1, 2], [0, 1, 3]] corresponds to a triangle
  connecting vertices 0, 1, and 2 and another triangle connecting vertices
  0, 1, and 3.
+ triangle_centroids: Triangle centroid coordinates.
  boundary_indices: Indices corresponding to the boundary.
  vertex_areas: The areas corresponding to the sites or vertices.
  triangle_areas: The areas of the triangular mesh elements.
@@ -39,6 +40,7 @@ def __init__(
  self,
  sites: Sequence[Tuple[float, float]],
  elements: Sequence[Tuple[int, int, int]],
+ triangle_centroids: Sequence[Tuple[int, int, int]],
  boundary_indices: Sequence[int],
  vertex_areas: Sequence[float],
  triangle_areas: Sequence[float],
@@ -47,6 +49,7 @@ def __init__(
  ):
  self.sites = np.asarray(sites).squeeze()
  self.elements = np.asarray(elements, dtype=np.int64)
+ self.triangle_centroids = np.asarray(triangle_centroids)
  self.boundary_indices = np.asarray(boundary_indices, dtype=np.int64)
  self.vertex_areas = np.asarray(vertex_areas)
  self.triangle_areas = np.asarray(triangle_areas)
@@ -138,10 +141,12 @@ def from_triangulation(
  boundary_indices = Mesh.find_boundary_indices(elements)
  edge_mesh = EdgeMesh.from_mesh(sites, elements)
  triangle_areas = utils.triangle_areas(sites, elements)
+ centroids = sites[elements].mean(axis=1)
  vertex_areas = utils.vertex_areas(sites, elements, tri_areas=triangle_areas)
  return Mesh(
  sites=sites,
  elements=elements,
+ triangle_centroids=centroids,
  boundary_indices=boundary_indices,
  edge_mesh=edge_mesh,
  vertex_areas=vertex_areas,
@@ -252,6 +257,7 @@ def to_hdf5(self, h5group: h5py.Group, compress: bool = True) -> None:
  h5group["sites"] = self.sites
  h5group["elements"] = self.elements
  if not compress:
+ h5group["triangle_centroids"] = self.triangle_centroids
  h5group["boundary_indices"] = self.boundary_indices
  h5group["vertex_areas"] = self.vertex_areas
  h5group["triangle_areas"] = self.triangle_areas
@@ -274,6 +280,7 @@ def from_hdf5(h5group: h5py.Group) -> "Mesh":
  return Mesh(
  sites=np.array(h5group["sites"]),
  elements=np.array(h5group["elements"], dtype=np.int64),
+ triangle_centroids=np.array(h5group["triangle_centroids"]),
  boundary_indices=np.array(h5group["boundary_indices"], dtype=np.int64),
  vertex_areas=np.array(h5group["vertex_areas"]),
  triangle_areas=np.array(h5group["triangle_areas"]),
@@ -299,6 +306,7 @@ def is_restorable(h5group: h5py.Group) -> bool:
  return (
  "sites" in h5group
  and "elements" in h5group
+ and "triangle_centroids" in h5group
  and "boundary_indices" in h5group
  and "vertex_areas" in h5group
  and "triangle_areas" in h5group
@@ -309,6 +317,7 @@ def copy(self) -> "Mesh":
  mesh = Mesh(
  sites=self.sites.copy(),
  elements=self.elements.copy(),
+ triangle_centroids=self.triangle_centroids.copy(),
  boundary_indices=self.boundary_indices.copy(),
  vertex_areas=self.vertex_areas.copy(),
  triangle_areas=self.triangle_areas.copy(),
@@ -337,12 +346,16 @@ def __init__(
  Q: np.ndarray,
  gradient_x: sp.spmatrix,
  gradient_y: sp.spmatrix,
+ gradient_tri_x: sp.spmatrix,
+ gradient_tri_y: sp.spmatrix,
  laplacian: sp.spmatrix,
  ):
  self.weights = weights
  self.Q = Q
  self.gradient_x = gradient_x
  self.gradient_y = gradient_y
+ self.gradient_tri_x = gradient_tri_x
+ self.gradient_tri_y = gradient_tri_y
  self.laplacian = laplacian
 
  @staticmethod
@@ -360,8 +373,11 @@ def from_mesh(mesh: Mesh) -> "MeshOperators":
  elements = mesh.elements
  weights = mesh.vertex_areas
  Q = MeshOperators.Q_matrix(sites, weights)
+ gradient_tri_x, gradient_tri_y = gradient_triangles(
+ sites, elements, mesh.triangle_areas
+ )
  gradient_x, gradient_y = gradient_vertices(
- sites, elements, areas=mesh.triangle_areas
+ sites, elements, gradient_tri=(gradient_tri_x, gradient_tri_y)
  )
  # gradient_edges = gradient_edges(
  # sites, mesh.edge_mesh.edges, mesh.edge_mesh.edge_lengths
@@ -372,6 +388,8 @@ def from_mesh(mesh: Mesh) -> "MeshOperators":
  Q=Q,
  gradient_x=gradient_x,
  gradient_y=gradient_y,
+ gradient_tri_x=gradient_tri_x,
+ gradient_tri_y=gradient_tri_y,
  laplacian=laplacian,
  )
 

superscreen/fem.py

@@ -350,6 +350,7 @@ def gradient_triangles(
 def gradient_vertices(
  points: np.ndarray,
  triangles: np.ndarray,
+ gradient_tri: Optional[Tuple[sp.csr_array, sp.csr_array]] = None,
  areas: Optional[np.ndarray] = None,
 ) -> Tuple[sp.csr_array, sp.csr_array]:
  """Returns the vertex gradient operators ``gx`` and ``gy``.
@@ -365,15 +366,19 @@ def gradient_vertices(
  Args:
  points: Shape (n, 2) array of x, y coordinates of vertices
  triangles: Shape (m, 3) array of triangle indices
+ gradient_tri: Pre-computed triangle gradient operators for x and y
  areas: Shape (m, ) array of triangle areas
 
  Returns:
  x and y gradient matrices, both of which have shape ``(n, n)``
  """
- if areas is None:
- areas = triangle_areas(points, triangles)
+ if gradient_tri is None:
+ if areas is None:
+ areas = triangle_areas(points, triangles)
+ Gx, Gy = gradient_triangles(points, triangles, areas=areas)
+ else:
+ Gx, Gy = gradient_tri
  n = len(points)
- Gx, Gy = gradient_triangles(points, triangles, areas=areas)
  Gx = Gx.tolil()
  Gy = Gy.tolil()
  gx = sp.lil_array((n, n), dtype=float)

superscreen/fluxoid.py

@@ -5,7 +5,7 @@
 
 from .device import Device
 from .solution import Solution
-from .solver import solve
+from .solver import FactorizedModel, solve
 
 logger = logging.getLogger(__name__)
 
@@ -53,16 +53,15 @@ def make_fluxoid_polygons(
 
 
 def find_fluxoid_solution(
- device: Device,
- *,
+ model: FactorizedModel,
  fluxoids: Optional[Dict[str, float]] = None,
  **solve_kwargs,
 ) -> Solution:
  """Calculates the current(s) circulating around hole(s) in the device required
  to realize the specified fluxoid state.
 
  Args:
- device: The Device for which to find the given fluxoid solution.
+ model: The factorized model for which to find the given fluxoid solution.
  fluxoids: A dict of ``{hole_name: fluxoid_value}``, where ``fluxoid_value`` is
  in units of :math:`\\Phi_0`. The fluxoid for any holes not in this dict
  will default to 0.
@@ -72,38 +71,47 @@ def find_fluxoid_solution(
  Returns:
  The optimized :class:`superscreen.Solution`.
  """
+ device = model.device
  fluxoids = fluxoids or {}
  holes = list(device.holes)
+ current_units = model.current_units
  solve_kwargs = solve_kwargs.copy()
- current_units = solve_kwargs.get("current_units", "uA")
- terminal_currents = solve_kwargs.pop("terminal_currents", None)
  applied_field = solve_kwargs.pop("applied_field", None)
- M = device.mutual_inductance_matrix(
- units=f"Phi_0 / {current_units}", **solve_kwargs
- )
  target_fluxoids = np.array([fluxoids.get(name, 0) for name in holes])
 
+ model = model.copy()
+ model.circulating_currents = {}
+
  # Find the hole fluxoids assuming no circulating currents.
  solution_no_circ = solve(
- device,
+ model=model,
  applied_field=applied_field,
- terminal_currents=terminal_currents,
- circulating_currents=None,
  **solve_kwargs,
  )[-1]
+
+ if not holes:
+ if np.any(target_fluxoids):
+ raise ValueError(
+ "Cannot calculate nonzero fluxoid solution for a device with no holes."
+ )
+ return solution_no_circ
+
  fluxoids = [
  sum(solution_no_circ.hole_fluxoid(name)).to("Phi_0").magnitude for name in holes
  ]
  fluxoids = np.array(fluxoids)
+
+ M = device.mutual_inductance_matrix(
+ units=f"Phi_0 / {current_units}", **solve_kwargs
+ )
  # Solve for the circulating currents needed to realize the target_fluxoids.
  I_circ = np.linalg.solve(M.magnitude, target_fluxoids - fluxoids)
  circulating_currents = dict(zip(holes, I_circ))
  # Solve the model with the optimized circulating currents.
+ model.set_circulating_currents(circulating_currents)
  solution = solve(
- device,
+ model=model,
  applied_field=applied_field,
- terminal_currents=terminal_currents,
- circulating_currents=circulating_currents,
  **solve_kwargs,
  )[-1]
  return solution

superscreen/geometry.py

@@ -25,7 +25,6 @@ def path_vectors(path: np.ndarray) -> Tuple[float, np.ndarray]:
  normals = np.cross(dr, [0, 0, 1])
  unit_normals = unit_vector(normals)
  edge_lengths = la.norm(dr, axis=1)
- unit_normals = np.concatenate([unit_normals, unit_normals[-1:]], axis=0)
  return edge_lengths, unit_normals[:, :2]
 
 

superscreen/solution.py

@@ -355,7 +355,7 @@ def current_through_path(
  )
  edge_lengths, unit_normals = path_vectors(path_coords)
  edge_lengths = edge_lengths * device.ureg(device.length_units)
- J_dot_n = (J_edge * unit_normals[:-1]).sum(axis=1)
+ J_dot_n = np.sum(J_edge * unit_normals, axis=1)
  total_current = np.trapz(J_dot_n * edge_lengths).to(units)
  if not with_units:
  total_current = total_current.magnitude

superscreen/solver/solve.py

@@ -1,3 +1,4 @@
+import copy
 import itertools
 import logging
 import os
@@ -200,6 +201,24 @@ def set_circulating_currents(self, circulating_currents: Dict[str, float]) -> No
  if hole in holes:
  film_info.circulating_currents[hole] = current
 
+ def set_vortices(self, vortices: Sequence[Vortex]) -> None:
+ """Set the vortices for the model.
+
+ Args:
+ vortices: A sequence of vortices.
+ """
+ for film in self.film_info.values():
+ film.vortices = []
+ for vortex in vortices:
+ self.film_info[vortex.film].vortices.append(vortex)
+ self.vortices = {}
+ for name, film in self.film_info.items():
+ film.vortices = tuple(film.vortices)
+ self.vortices[name] = film.vortices
+
+ def copy(self) -> "FactorizedModel":
+ return copy.copy(self)
+
 
 def factorize_model(
  *,
@@ -357,12 +376,12 @@ def solve(
  "If model argument is provided, device, terminal_currents,"
  " circulating_currents, and vortices must be None."
  )
- if current_units is not None and current_units != model.current_units:
- logger.warning(
- "Keyword argument 'current_units' is ignored when "
- "a factorized model is provided. "
- f"Using model.current_units = {model.current_units!r}"
- )
+ # if current_units is not None and current_units != model.current_units:
+ #  logger.warning(
+ #  "Keyword argument 'current_units' is ignored when "
+ #  "a factorized model is provided. "
+ #  f"Using model.current_units = {model.current_units!r}"
+ #  )
 
  if not isinstance(model, FactorizedModel):
  raise TypeError(