Add device to FactorizedModel (#99)

* Add device to FactorizedModel

* Update docs

docs/about/changelog.rst

@@ -14,6 +14,9 @@ View release history on `PyPI <https://pypi.org/project/superscreen/#history>`_,
       The public API SHOULD NOT be considered stable.
     - Version 1.0.0 defines the public API.
 
+.. contents::
+    :depth: 2
+
 ----
 
 Version 0.8.1

docs/api/device.rst

@@ -14,6 +14,9 @@ and ``superscreen.parameter`` are used to set up the inputs to a
 - The geometry and penetration depth of all superconducting films
 - The spatial distribution of the applied magnetic field
 
+.. contents::
+    :depth: 2
+
 Devices
 -------
 

docs/api/solve.rst

@@ -10,11 +10,23 @@ Solver
 The ``superscreen.solve`` module contains the actual implementation
 of Brandt's method, as described :ref:`here <background>`.
 
+.. contents::
+    :depth: 2
+
 Solve
 -----
 
 .. autofunction:: superscreen.solve
 
+.. autoclass:: superscreen.Vortex
+
+Pre-factorizing models
+----------------------
+
+.. autofunction:: superscreen.factorize_model
+
+.. autoclass:: superscreen.FactorizedModel
+    :members:
 
 Solution
 --------
@@ -31,8 +43,6 @@ Fluxoid
 .. autoclass:: superscreen.Fluxoid
     :show-inheritance:
 
-.. autoclass:: superscreen.Vortex
-
 .. autofunction:: superscreen.fluxoid.make_fluxoid_polygons
 
 .. autofunction:: superscreen.fluxoid.find_fluxoid_solution
@@ -45,8 +55,6 @@ Supporting functions
 
 .. autofunction:: superscreen.solver.solve_for_terminal_current_stream
 
-.. autofunction:: superscreen.solver.factorize_model
-
 .. autofunction:: superscreen.solver.factorize_linear_systems
 
 .. autofunction:: superscreen.solver.convert_field
@@ -57,9 +65,6 @@ Supporting functions
 Supporting classes
 ------------------
 
-.. autoclass:: superscreen.solver.FactorizedModel
-    :members:
-
 .. autoclass:: superscreen.solver.FilmInfo
     :members:
 

docs/api/sources.rst

@@ -10,6 +10,9 @@ Field Sources
 The ``superscreen.sources`` module provides factory functions for a few convenient
 :class:`superscreen.parameter.Parameter` classes used to define applied magnetic fields. 
 
+.. contents::
+    :depth: 2
+
 ConstantField
 -------------
 

docs/api/visualization.rst

@@ -10,6 +10,9 @@ Visualization
 The ``superscreen.visualization`` module contains functions used to visualize
 the results contained in a :class:`superscreen.Solution`.
 
+.. contents::
+    :depth: 2
+
 Fields
 ------
 

docs/index.rst

@@ -102,8 +102,8 @@ If you use ``SuperScreen`` in your research, please cite the paper linked above.
    :maxdepth: 2
    :caption: Tutorials
 
-   notebooks/polygons.ipynb
    notebooks/terminal-currents.ipynb
+   notebooks/polygons.ipynb
    notebooks/field-sources.ipynb
    notebooks/logo.ipynb
 

superscreen/__init__.py

@@ -4,7 +4,7 @@
 from .fluxoid import find_fluxoid_solution, make_fluxoid_polygons
 from .parameter import Constant, Parameter
 from .solution import FilmSolution, Fluxoid, Solution, Vortex
-from .solver import convert_field, solve
+from .solver import FactorizedModel, convert_field, factorize_model, solve
 from .units import ureg
 from .version import __version__, __version_info__
 from .visualization import (

superscreen/solver/solve.py

@@ -77,6 +77,7 @@ class FactorizedModel:
     """A pre-factorized SuperScreen model.
 
     Args:
+        device: The :class:`superscreen.Device`
         film_info: A dict of ``{film_name: FilmInfo}``
         film_systems: A dict of ``{film_name: LinearSystem}``
         hole_systems: A dict of ``{film_name: {hole_name: LinearSystem}}``
@@ -86,6 +87,7 @@ class FactorizedModel:
         vortices: A dict of ``{film_name: vortices}``
     """
 
+    device: Device
     film_info: Dict[str, FilmInfo]
     film_systems: Dict[str, LinearSystem]
     hole_systems: Dict[str, Dict[str, LinearSystem]]
@@ -95,11 +97,12 @@ class FactorizedModel:
     vortices: Dict[str, Sequence[Vortex]]
 
     def to_hdf5(self, h5group: h5py.Group) -> None:
-        """Save a :class:`superscreen.solver.FactorizedModel` to an :class:`h5py.Group`.
+        """Save a :class:`superscreen.FactorizedModel` to an :class:`h5py.Group`.
 
         Args:
             h5group: The :class:`h5py.Group` in which to save the model
         """
+        self.device.to_hdf5(h5group.create_group("device"))
         film_info_grp = h5group.create_group("film_info")
         for film, info in self.film_info.items():
             info.to_hdf5(film_info_grp.create_group(film))
@@ -126,14 +129,15 @@ def to_hdf5(self, h5group: h5py.Group) -> None:
 
     @staticmethod
     def from_hdf5(h5group: h5py.Group) -> "FactorizedModel":
-        """Load a :class:`superscreen.solver.FactorizedModel` from an :class:`h5py.Group`.
+        """Load a :class:`superscreen.FactorizedModel` from an :class:`h5py.Group`.
 
         Args:
             h5group: The :class:`h5py.Group` from which to load the model
 
         Returns:
-            The loaded :class:`superscreen.solver.FactorizedModel`
+            The loaded :class:`superscreen.FactorizedModel`
         """
+        device = Device.from_hdf5(h5group["device"])
         film_info = {
             film: FilmInfo.from_hdf5(grp) for film, grp in h5group["film_info"].items()
         }
@@ -159,6 +163,7 @@ def from_hdf5(h5group: h5py.Group) -> "FactorizedModel":
             Vortex.from_hdf5(vortex_grp[i]) for i in sorted(vortex_grp, key=int)
         ]
         return FactorizedModel(
+            device=device,
             film_info=film_info,
             film_systems=film_systems,
             hole_systems=hole_systems,
@@ -196,7 +201,7 @@ def factorize_model(
             located in the :class:`superscreen.Device`.
 
     Returns:
-        A :class:`superscreen.solver.FactorizedModel` instance that can be used to solve the model.
+        A :class:`superscreen.FactorizedModel` instance that can be used to solve the model.
     """
     ureg = device.ureg
     circulating_currents = circulating_currents or {}
@@ -224,6 +229,7 @@ def factorize_model(
         device, film_info
     )
     return FactorizedModel(
+        device,
         film_info,
         film_systems,
         hole_systems,
@@ -235,7 +241,7 @@ def factorize_model(
 
 
 def solve(
-    device: Device,
+    device: Optional[Device] = None,
     *,
     model: Optional[FactorizedModel] = None,
     applied_field: Optional[Callable] = None,
@@ -269,7 +275,8 @@ def solve(
 
 
     Args:
-        device: The Device to simulate.
+        device: The Device to simulate. Required if ``model`` is not provided.
+        model: A :class:`superscreen.FactorizedModel` instance.
         applied_field: A callable that computes the applied magnetic field
             as a function of x, y, z coordinates.
         terminal_currents: A dict like ``{film_name: {source_name: source_current}}`` for
@@ -298,28 +305,9 @@ def solve(
     if log_level is not None:
         logging.basicConfig(level=log_level)
 
-    if not device.meshes:
-        raise ValueError(
-            "The device does not have a mesh. Call device.make_mesh() to generate it."
-        )
-
-    dtype = device.solve_dtype
-    ureg = device.ureg
-    length_units = device.length_units
-    meshes = device.meshes
-    applied_field = applied_field or ConstantField(0)
-    field_conversion = field_conversion_factor(
-        field_units,
-        current_units,
-        length_units=length_units,
-        ureg=ureg,
-    )
-    logger.debug(
-        f"Conversion factor from {ureg(field_units).units:~P} to "
-        f"{ureg(current_units) / ureg(length_units):~P}: {field_conversion:~P}."
-    )
-
     if model is None:
+        if device is None:
+            raise ValueError("Either a model or a device must be provided.")
         logger.info("Factorizing model.")
         model = factorize_model(
             device=device,
@@ -329,18 +317,22 @@ def solve(
             vortices=vortices,
         )
     elif (
-        terminal_currents is not None
+        device is not None
+        or terminal_currents is not None
         or circulating_currents is not None
         or vortices is not None
     ):
         raise ValueError(
-            "If model argument is provided, terminal_currents, circulating_currents,"
-            " and vortices must be None."
+            "If model argument is provided, device, terminal_currents,"
+            " circulating_currents, and vortices must be None."
         )
+
     if not isinstance(model, FactorizedModel):
         raise TypeError(
             f"model must be an instance of FactorizedModel (got {type(model)})."
         )
+
+    device = model.device
     film_info = model.film_info
     film_systems = model.film_systems
     hole_systems = model.hole_systems
@@ -349,6 +341,27 @@ def solve(
     circulating_currents = model.circulating_currents
     vortices = model.vortices
 
+    if not device.meshes:
+        raise ValueError(
+            "The device does not have a mesh. Call device.make_mesh() to generate it."
+        )
+
+    dtype = device.solve_dtype
+    ureg = device.ureg
+    length_units = device.length_units
+    meshes = device.meshes
+    applied_field = applied_field or ConstantField(0)
+    field_conversion = field_conversion_factor(
+        field_units,
+        current_units,
+        length_units=length_units,
+        ureg=ureg,
+    )
+    logger.debug(
+        f"Conversion factor from {ureg(field_units).units:~P} to "
+        f"{ureg(current_units) / ureg(length_units):~P}: {field_conversion:~P}."
+    )
+
     applied_fields = {}
     for film, mesh in meshes.items():
         layer = device.layers[film_info[film].layer]
@@ -391,7 +404,7 @@ def solve(
             film_info=film_info[film_name],
             field_conversion=field_conversion.magnitude,
             vortex_flux=vortex_flux,
-            terminal_systems=terminal_systems[film_name],
+            terminal_systems=terminal_systems.get(film_name, None),
             check_inversion=check_inversion,
         )
 
@@ -453,7 +466,7 @@ def solve(
                 film_info=film_info[film_name],
                 field_conversion=field_conversion.magnitude,
                 vortex_flux=vortex_flux,
-                terminal_systems=terminal_systems[film_name],
+                terminal_systems=terminal_systems.get(film_name, None),
                 check_inversion=check_inversion,
             )
 

superscreen/solver/solve_film.py

@@ -151,7 +151,7 @@ def factorize_linear_systems(
 ) -> Tuple[
     Dict[str, LinearSystem],
     Dict[str, Dict[str, LinearSystem]],
-    Dict[str, Optional[TerminalSystems]],
+    Dict[str, TerminalSystems],
 ]:
     """Build and factorize the linear systems for all films, holes, and terminals.
 
@@ -216,7 +216,6 @@ def make_system_2d(indices, Lambda):
                 grad_Lambda_term=grad_Lambda_term,
             )
 
-        terminal_systems[film_name] = None
         if film_name in device.terminals:
             # Make the boundary linear system
             boundary_system = LinearSystem(

superscreen/solver/utils.py

@@ -81,16 +81,14 @@ def from_hdf5(h5group: h5py.Group) -> "LambdaInfo":
         Returns:
             The loaded ``LambdaInfo`` instance
         """
-        london_lambda = thickness = None
+        london_lambda = None
         if "london_lambda" in h5group:
             london_lambda = np.array(h5group["london_lambda"])
-        if "thickness" in h5group:
-            thickness = float(h5group["thickness"])
         return LambdaInfo(
             film=h5group.attrs["film"],
             Lambda=np.array(h5group["Lambda"]),
             london_lambda=london_lambda,
-            thickness=thickness,
+            thickness=h5group.attrs.get("thickness", None),
         )
 
 

superscreen/test/test_solve.py

@@ -1,3 +1,4 @@
+import h5py
 import matplotlib.pyplot as plt
 import numpy as np
 import pint
@@ -119,13 +120,19 @@ def linear(x, y, offset=0):
                 linear, offset=old_lambda
             )
         if pre_factorize:
-            model = sc.solver.factorize_model(
+            model = sc.factorize_model(
                 device=device,
                 circulating_currents=circulating_currents,
                 current_units="uA",
             )
+            model_save_path = tmp_path / "model.h5"
+            with h5py.File(model_save_path, "x") as h5file:
+                model.to_hdf5(h5file)
+
+            with h5py.File(model_save_path, "r") as h5file:
+                model = sc.FactorizedModel.from_hdf5(h5file)
+
             solutions = sc.solve(
-                device=device,
                 model=model,
                 applied_field=applied_field,
                 field_units="mT",