Fix bug when calling find_fluxoid_solution multiple times with the same model (#117)

Author: loganbvh

superscreen/device/device.py

@@ -574,6 +574,7 @@ def mutual_inductance_matrix(
         from ..solver import factorize_model, solve
 
         holes = self.holes
+        hole_names = list(self.holes)
         if hole_polygon_mapping is None:
             from ..fluxoid import make_fluxoid_polygons
 
@@ -607,11 +608,11 @@ def mutual_inductance_matrix(
                 films_by_hole[hole.name] = film
         model = None
         for j, hole_name in enumerate(
-            tqdm(hole_polygon_mapping, desc="Holes", disable=(not progress_bar))
+            tqdm(hole_names, desc="Holes", disable=(not progress_bar))
         ):
             logger.info(
                 f"Evaluating {self.name!r} mutual inductance matrix "
-                f"column ({j+1}/{len(hole_polygon_mapping)}), "
+                f"column ({j+1}/{len(hole_names)}), "
                 f"source = {hole_name!r}."
             )
             if model is None:
@@ -623,15 +624,15 @@ def mutual_inductance_matrix(
                 I_circ_val = model.circulating_currents[hole_name]
             else:
                 model.set_circulating_currents({hole_name: I_circ_val})
-            solutions = solve(device=None, model=model, **solve_kwargs)[solution_slice]
+            solutions = solve(model=model, **solve_kwargs)[solution_slice]
 
             for n, solution in enumerate(solutions):
                 logger.info(
                     f"Evaluating fluxoids for solution {n + 1}/{len(solutions)}."
                 )
-                for i, (name, polygon) in enumerate(hole_polygon_mapping.items()):
+                for i, name in enumerate(hole_names):
                     fluxoid = solution.polygon_fluxoid(
-                        polygon, film=films_by_hole[name]
+                        hole_polygon_mapping[name], film=films_by_hole[name]
                     )
                     mutual_inductance[n, i, j] = (
                         (sum(fluxoid) / I_circ).to(units).magnitude

superscreen/fluxoid.py

@@ -73,45 +73,47 @@ def find_fluxoid_solution(
     """
     device = model.device
     fluxoids = fluxoids or {}
-    holes = list(device.holes)
+    hole_names = list(device.holes)
     current_units = model.current_units
+    inductance_units = f"Phi_0 / {current_units}"
     solve_kwargs = solve_kwargs.copy()
     applied_field = solve_kwargs.pop("applied_field", None)
-    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(
-        model=model,
-        applied_field=applied_field,
-        **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(
-        model=model,
-        applied_field=applied_field,
-        **solve_kwargs,
-    )[-1]
+    target_fluxoids = np.array([fluxoids.get(name, 0) for name in hole_names])
+
+    orig_circulating_currents = model.circulating_currents
+    try:
+        # Find the hole fluxoids assuming no circulating currents.
+        model.set_circulating_currents({name: 0 for name in hole_names})
+        solution_no_circ = solve(
+            model=model,
+            applied_field=applied_field,
+            **solve_kwargs,
+        )[-1]
+
+        if not hole_names:
+            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 hole_names
+        ]
+        fluxoids = np.array(fluxoids)
+        M = device.mutual_inductance_matrix(units=inductance_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(hole_names, I_circ))
+        # Solve the model with the optimized circulating currents.
+        model.set_circulating_currents(circulating_currents)
+        solution = solve(
+            model=model,
+            applied_field=applied_field,
+            **solve_kwargs,
+        )[-1]
+    finally:
+        model.set_circulating_currents(orig_circulating_currents)
     return solution

superscreen/version.py

@@ -1,2 +1,2 @@
-__version_info__ = (0, 12, 0)
+__version_info__ = (0, 12, 1)
 __version__ = ".".join(map(str, __version_info__))