# Solver¶

The superscreen.solve module contains the actual implementation of Brandt’s method, as described here.

## Solve¶

superscreen.solve(device, *, applied_field=None, terminal_currents=None, circulating_currents=None, vortices=None, field_units='mT', current_units='uA', check_inversion=False, iterations=0, return_solutions=True, directory=None, cache_memory_cutoff=inf, log_level=20, gpu=False, _solver='superscreen.solve')[source]

Computes the stream functions and magnetic fields for all layers in a Device.

The simulation strategy is:

1. Compute the stream functions and fields for each layer given only the applied field.

2. If iterations > 1 and there are multiple layers, then for each layer, calculate the screening field from all other layers and recompute the stream function and fields based on the sum of the applied field and the screening fields from all other layers.

1. Repeat step 2 (iterations - 1) times.

Parameters:
Return type:
Returns:

If return_solutions is True, returns a list of Solutions of length iterations + 1.

## Solve Many¶

superscreen.solve_many(device, *, parallel_method=None, applied_fields=None, circulating_currents=None, vortices=None, layer_updater=None, layer_update_kwargs=None, field_units='mT', current_units='uA', check_inversion=False, iterations=0, product=False, directory=None, return_solutions=False, keep_only_final_solution=False, cache_memory_cutoff=inf, log_level=20, use_shared_memory=True, num_cpus=None, gpu=False)[source]

Solves many models involving the same device, optionally in parallel using multiple processes.

Parameters:
Return type:
Returns:

solutions, paths. If return_solutions is True, solutions is either a list of lists of Solutions (if keep_only_final_solution is False), or a list of Solutions (the final iteration for each setup). If directory is True, paths is a list of paths to the saved solutions, otherwise paths is None.

## Solution¶

class superscreen.Solution(*, device, streams, current_densities, fields, screening_fields, applied_field, field_units, current_units, circulating_currents=None, terminal_currents=None, vortices=None, solver='superscreen.solve')[source]

A container for the calculated stream functions and fields, with some convenient data processing methods.

Parameters:
property field_units: str

The units in which magnetic fields are specified.

Return type:

str

property current_units: str

The units in which currents are specified.

Return type:

str

property solver: str

The solver method that generated the solution.

Return type:

str

property time_created: datetime

The time at which the solution was originally created.

Return type:

datetime

property version_info: Dict[str, str]

A dictionary of dependency versions.

Return type:
grid_data(dataset, *, layers=None, grid_shape=(200, 200), method='linear', with_units=False, **kwargs)[source]

Interpolates results from the triangular mesh to a rectangular grid.

Keyword arguments are passed to scipy.interpolate.griddata().

Parameters:
Return type:
Returns:

x grid, y grid, dict of interpolated data for each layer

grid_current_density(*, layers=None, grid_shape=(200, 200), method='linear', units=None, with_units=False, **kwargs)[source]

Computes the current density J = [dg/dy, -dg/dx] on a rectangular grid.

Keyword arguments are passed to scipy.interpolate.griddata().

Parameters:
Return type:
Returns:

x grid, y grid, dict of interpolated current density for each layer

interp_current_density(positions, *, layers=None, grid_shape=(200, 200), method='linear', units=None, with_units=False, **kwargs)[source]

Computes the current density J = [dg/dy, -dg/dx] at unstructured coordinates via interpolation.

Keyword arguments are passed to scipy.interpolate.griddata().

Parameters:
Returns:

A dict of interpolated current density for each layer.

interp_fields(positions, *, layers=None, method='linear', units=None, with_units=False, **kwargs)[source]

Interpolates the fields in one or more layers.

Additional keyword arguments are passed to the relevant interpolator: :classscipy.interpolate.NearestNDInterpolator, scipy.interpolate.LinearNDInterpolator, or scipy.interpolate.CloughTocher2DInterpolator.

Parameters:
Returns:

A dict of interpolated fields for each layer.

polygon_flux(*, polygons=None, units=None, with_units=True)[source]

Computes the flux through all polygons (films, holes, and flux regions) by integrating the calculated fields.

Parameters:
Return type:

Dict[str, Union[float, Quantity]]

Returns:

A dict of {polygon_name: polygon_flux}

polygon_fluxoid(polygon_points, layers=None, grid_shape=(200, 200), interp_method='linear', units='Phi_0', with_units=True)[source]

Computes the Fluxoid (flux + supercurrent) for a given polygonal region.

The fluxoid for a closed region $$S$$ with boundary $$\partial S$$ is defined as:

$\Phi^f_S = \underbrace{ \int_S \mu_0 H_z(\vec{r})\,\mathrm{d}^2r }_{\text{flux part}} + \underbrace{ \oint_{\partial S} \mu_0\Lambda(\vec{r})\vec{J}(\vec{r})\cdot\mathrm{d}\vec{r} }_{\text{supercurrent part}}$
Parameters:
Return type:
Returns:

A dict of {layer_name: fluxoid} for each specified layer, where fluxoid is an instance of Fluxoid.

hole_fluxoid(hole_name, points=None, grid_shape=(200, 200), interp_method='linear', units='Phi_0', with_units=True)[source]

Calculcates the fluxoid for a polygon enclosing the specified hole.

Parameters:
Return type:

Fluxoid

Returns:

The hole’s Fluxoid.

field_at_position(positions, *, zs=None, vector=False, units=None, with_units=True, return_sum=True)[source]

Calculates the field due to currents in the device at any point(s) in space.

Parameters:
• positions (ndarray) – Shape (m, 2) array of (x, y) coordinates, or (m, 3) array of (x, y, z) coordinates at which to calculate the magnetic field. A single sequence like [x, y] or [x, y, z] is also allowed.

• zs (Union[float, ndarray, None]) – z coordinates at which to calculate the field. If positions has shape (m, 3), then this argument is not allowed. If zs is a scalar, then the fields are calculated in a plane parallel to the x-y plane. If zs is any array, then it must be same length as positions.

• vector (bool) – Whether to return the full vector magnetic field or just the z component.

• units (Optional[str]) – Units to which to convert the fields (can be either magnetic field H or magnetic flux density B = mu0 * H). If not given, then the fields are returned in units of self.field_units.

• with_units (bool) – Whether to return the fields as pint.Quantity with units attached.

• return_sum (bool) – Whether to return the sum of the fields from all layers in the device, or a dict of {layer_name: field_from_layer}.

Return type:
Returns:

An np.ndarray if return_sum is True, otherwise a dict of {layer_name: field_from_layer}. If with_units is True, then the array(s) will contain pint.Quantities. field_from_layer will have shape (m, ) if vector is False, or shape (m, 3) if vector is True.

vector_potential_at_position(positions, *, zs=None, units=None, with_units=True, return_sum=True)[source]

Calculates the vector potential due to currents in the device at any point(s) in space. Note that this only considers the vector potential due to currents in the device, so only represents the total vector potential in cases where the applied field is zero (e.g. models with only vortices and/or circulating currents).

The vector potential $$\vec{A}$$ at position $$\vec{r}$$ due to sheet current density $$\vec{J}(\vec{r}')$$ flowing in a film with lateral geometry $$S$$ is:

$\vec{A}(\vec{r}) = \frac{\mu_0}{4\pi} \int_S\frac{\vec{J}(\vec{r}')}{|\vec{r}-\vec{r}'|}\mathrm{d}^2r'.$
Parameters:
• positions (ndarray) – Shape (m, 2) array of (x, y) coordinates, or (m, 3) array of (x, y, z) coordinates at which to calculate the vector potential. A single list like [x, y] or [x, y, z] is also allowed.

• zs (Union[float, ndarray, None]) – z coordinates at which to calculate the potential. If positions has shape (m, 3), then this argument is not allowed. If zs is a scalar, then the fields are calculated in a plane parallel to the x-y plane. If zs is any array, then it must be same length as positions.

• units (Optional[str]) – Units to which to convert the vector potential.

• with_units (bool) – Whether to return the vector potential as a pint.Quantity with units attached.

• return_sum (bool) – Whether to return the sum of the potential from all layers in the device, or a dict of {layer_name: potential_from_layer}.

Return type:
Returns:

An np.ndarray if return_sum is True, otherwise a dict of {layer_name: potential_from_layer}. If with_units is True, then the array(s) will contain pint.Quantities. potential_from_layer will have shape (m, 3).

to_file(directory, save_mesh=True, compressed=True, to_zip=False)[source]

Saves a Solution to disk.

Parameters:
• directory (str) – The name of the directory in which to save the solution (must either be empty or not yet exist).

• save_mesh (bool) – Whether to save the device mesh.

• compressed (bool) – Whether to use numpy.savez_compressed rather than numpy.savez.

• to_zip (bool) – Whether to save the Solution to a zip file.

Return type:

None

classmethod from_file(directory, compute_matrices=False)[source]

Parameters:
• directory (str) – The directory from which to load the solution.

• compute_matrices (bool) – Whether to compute the field-independent matrices for the device if the mesh already exists.

Return type:

Solution

Returns:

equals(other, require_same_timestamp=False)[source]

Checks whether two solutions are equal.

Parameters:
• other (Any) – The Solution to compare for equality.

• require_same_timestamp (bool) – If True, two solutions are only considered equal if they have the exact same time_created.

Return type:

bool

Returns:

A boolean indicating whether the two solutions are equal

plot_streams(**kwargs)[source]

Plots the stream function for multiple layers in a Device.

Additional keyword arguments are passed to plt.subplots().

Parameters:
• layers – Name(s) of layer(s) for which to plot the stream function. By default, the stream function is plotted for all layers in the Device.

• units – Units in which to plot the stream function. Defaults to solution.current_units.

• max_cols – Maximum number of columns in the grid of subplots.

• cmap – Name of the matplotlib colormap to use.

• levels – Number of contour levels to used.

• colorbar – Whether to add a colorbar to each subplot.

• filled – If True, plots filled contours.

Return type:
Returns:

matplotlib figure and axes

plot_currents(**kwargs)[source]

Plots the current density (sheet current) for each layer in a Device.

Additional keyword arguments are passed to plt.subplots().

Parameters:
• layers – Name(s) of layer(s) for which to plot the sheet current. By default, the stream function is plotted for all layers in the Device.

• units – Units in which to plot the current density. Defaults to solution.current_units / solution.device.length_units.

• grid_shape – Shape of the desired rectangular grid. If a single integer n is given, then the grid will be square, shape = (n, n).

• grid_method – Interpolation method to use (see scipy.interpolate.griddata).

• max_cols – Maximum number of columns in the grid of subplots.

• cmap – Name of the matplotlib colormap to use.

• colorbar – Whether to add a colorbar to each subplot.

• auto_range_cutoff – Cutoff percentile for auto_range_iqr.

• share_color_scale – Whether to force all layers to use the same color scale.

• symmetric_color_scale – Whether to use a symmetric color scale (vmin = -vmax).

• vmin – Color scale minimum to use for all layers (ignored if share_color_scale is True).

• vmax – Color scale maximum to use for all layers (ignored if share_color_scale is True).

• streamplot – Whether to overlay current streamlines on the plot.

• min_stream_amp – Streamlines will not be drawn anywhere the current density is less than min_stream_amp * max(current_density). This avoids streamlines being drawn where there is no current flowing.

• cross_section_coords – Shape (m, 2) array of (x, y) coordinates for a cross-section (or a list of such arrays).

Return type:
Returns:

matplotlib figure and axes

plot_fields(**kwargs)[source]

Plots either the total field or the screening field for multiple layers in a Device.

Additional keyword arguments are passed to plt.subplots().

Parameters:
• layers – Name(s) of layer(s) for which to plot the fields. By default, the stream function is plotted for all layers in the Device.

• dataset – Which set of fields to plot, either “fields” or “screening_fields”.

• normalize – Whether to normalize the fields by the applied field.

• units – Units in which to plot the fields. Defaults to solution.field_units. This argument is ignored if normalize is True.

• shading – May be "flat" or "gouraud". The latter does some interpolation.

• max_cols – Maximum number of columns in the grid of subplots.

• cmap – Name of the matplotlib colormap to use.

• colorbar – Whether to add a colorbar to each subplot.

• auto_range_cutoff – Cutoff percentile for auto_range_iqr.

• share_color_scale – Whether to force all layers to use the same color scale.

• symmetric_color_scale – Whether to use a symmetric color scale (vmin = -vmax).

• vmin – Color scale minimum to use for all layers

• vmax – Color scale maximum to use for all layers

• cross_section_coords – Shape (m, 2) array of (x, y) coordinates for a cross-section (or a list of such arrays).

Return type:
Returns:

matplotlib figure and axes

plot_field_at_positions(points, **kwargs)[source]

Plots the total field (either all three components or just the z component) at a given set of positions (x, y, z) outside of the device.

Additional keyword arguments are passed to plt.subplots().

Parameters:
• positions – Shape (m, 2) array of (x, y) coordinates, or (m, 3) array of (x, y, z) coordinates at which to calculate the magnetic field. A single list like [x, y] or [x, y, z] is also allowed.

• zs – z coordinates at which to calculate the field. If positions has shape (m, 3), then this argument is not allowed. If zs is a scalar, then the fields are calculated in a plane parallel to the x-y plane. If zs is any array, then it must be same length as positions.

• vector – Whether to return the full vector magnetic field or just the z component.

• units – Units in which to plot the fields. Defaults to solution.field_units. This argument is ignored if normalize is True.

• shading – May be "flat" or "gouraud". The latter does some interpolation.

• max_cols – Maximum number of columns in the grid of subplots.

• cmap – Name of the matplotlib colormap to use.

• colorbar – Whether to add a colorbar to each subplot.

• auto_range_cutoff – Cutoff percentile for auto_range_iqr.

• share_color_scale – Whether to force all layers to use the same color scale.

• symmetric_color_scale – Whether to use a symmetric color scale (vmin = -vmax).

• vmin – Color scale minimum to use for all layers

• vmax – Color scale maximum to use for all layers

• cross_section_coords – Shape (m, 2) array of (x, y) coordinates for a cross-section (or a list of such arrays).

• points (ndarray) –

Return type:
Returns:

matplotlib figure and axes

## Fluxoid¶

class superscreen.Fluxoid(flux_part, supercurrent_part)[source]

Bases: NamedTuple

The fluxoid for a closed region $$S$$ with boundary $$\partial S$$ is defined as:

$\Phi^f_S = \underbrace{ \int_S \mu_0 H_z(\vec{r})\,\mathrm{d}^2r }_{\text{flux part}} + \underbrace{ \oint_{\partial S} \mu_0\Lambda(\vec{r})\vec{J}(\vec{r})\cdot\mathrm{d}\vec{r} }_{\text{supercurrent part}}$
Parameters:
• flux_part (float | Quantity) – $$\int_S \mu_0 H_z(\vec{r})\,\mathrm{d}^2r$$.

• supercurrent_part (float | Quantity) – $$\oint_{\partial S}\mu_0\Lambda(\vec{r})\vec{J}(\vec{r})\cdot\mathrm{d}\vec{r}$$.

class superscreen.Vortex(x, y, layer, nPhi0=1)[source]

Bases: NamedTuple

A vortex located at position (x, y) in layer containing a total flux nPhi0 in units of the flux quantum $$\Phi_0$$.

Parameters:
• x (float) – Vortex x-position.

• y (float) – Vortex y-position.

• layer (str) – The layer in which the vortex is pinned.

• nPhi0 (float) – The number of flux quanta contained in the vortex.

superscreen.fluxoid.make_fluxoid_polygons(device, holes=None, join_style='mitre', interp_points=None)[source]

Generates polygons enclosing the given holes to calculate the fluxoid.

Parameters:
Return type:
Returns:

A dict of {hole_name: fluxoid_polygon}.

superscreen.fluxoid.find_fluxoid_solution(device, *, fluxoids=None, **solve_kwargs)[source]

Calculates the current(s) circulating around hole(s) in the device required to realize the specified fluxoid state.

Parameters:
Return type:

Solution

Returns:

The optimized superscreen.Solution.

## IO¶

superscreen.io.save_solutions(solutions, base_directory, save_mesh=True, return_paths=False, to_zip=False)[source]

Saves a list of Solutions to disk.

Parameters:
Return type:
Returns:

If return_paths is True, returns a list of paths where each solution was saved.

Loads a sequence of Solutions from disk.

Parameters:

base_directory (PathLike) – The name of the directory from which to load the solutions.

Return type:
Returns:

A list of Solutions

An iterator that loads a sequence of Solutions from disk.

Parameters:

base_directory (PathLike) – The name of the directory from which to load the solutions.

Yields:

Solution instances loaded from base_directory

Return type:
superscreen.io.zip_solution(solution, directory)[source]

Save a Solution to a zip archive in the given directory.

Parameters:
Return type:

str

Returns:

The absolute path to the created zip file.

superscreen.io.unzip_solution(path)[source]

Load a solution from a zip file.

Parameters:

path (PathLike) – The path to the zip file.

Return type:

Solution

Returns:

## Supporting Functions¶

### Brandt Core¶

superscreen.solve.solve_layer(*, device, layer, applied_field, kernel, weights, Del2, grad, Lambda_info, terminal_currents=None, circulating_currents=None, vortices=None, current_units='uA', check_inversion=False, gpu=False)[source]

Computes the stream function and magnetic field within a single layer of a Device.

Parameters:
Return type:
Returns:

stream function, current density, total field, film screening field

superscreen.solve.q_matrix(points, dtype=None)[source]

Computes the denominator matrix, q:

$q_{ij} = \frac{1}{4\pi|\vec{r}_i-\vec{r}_j|^3}$

See Eq. 7 in [Brandt-PRB-2005], Eq. 8 in [Kirtley-RSI-2016], and Eq. 8 in [Kirtley-SST-2016].

Parameters:
Return type:

ndarray

Returns:

Shape (n, n) array, qij

superscreen.solve.C_vector(points, dtype=None)[source]

Computes the edge vector, C:

$\begin{split}C_i &= \frac{1}{4\pi}\sum_{p,q=\pm1}\sqrt{(\Delta x - px_i)^{-2} + (\Delta y - qy_i)^{-2}}\\ \Delta x &= \frac{1}{2}(\mathrm{max}(x) - \mathrm{min}(x))\\ \Delta y &= \frac{1}{2}(\mathrm{max}(y) - \mathrm{min}(y))\end{split}$

See Eq. 12 in [Brandt-PRB-2005], Eq. 16 in [Kirtley-RSI-2016], and Eq. 15 in [Kirtley-SST-2016].

Parameters:
Return type:

ndarray

Returns:

Shape (n, ) array, Ci

superscreen.solve.Q_matrix(q, C, weights, dtype=None)[source]

Computes the kernel matrix, Q:

$Q_{ij} = (\delta_{ij}-1)q_{ij} + \delta_{ij}\frac{1}{w_{ij}}\left(C_i + \sum_{l\neq i}q_{il}w_{il}\right)$

See Eq. 10 in [Brandt-PRB-2005], Eq. 11 in [Kirtley-RSI-2016], and Eq. 11 in [Kirtley-SST-2016].

Parameters:
Return type:

ndarray

Returns:

Shape (n, n) array, Qij

superscreen.solve.convert_field(value, new_units, old_units=None, ureg=None, with_units=True)[source]

Converts a value between different field units, either magnetic field H [current] / [length] or flux density B = mu0 * H [mass] / ([curret] [time]^2)).

Parameters:
Return type:

Union[Quantity, ndarray, float]

Returns:

The converted value, either a pint.Quantity (scalar or array with units), or an array or float without units, depending on the with_units argument.

superscreen.solve.field_conversion_factor(field_units, current_units, length_units='m', ureg=None)[source]

Returns a conversion factor from field_units to current_units / length_units.

Parameters:
• field_units (str) – Magnetic field/flux unit to convert, having dimensionality either of magnetic field H (e.g. A / m or Oe) or of magnetic flux density B = mu0 * H (e.g. Tesla or Gauss).

• current_units (str) – Current unit to use for the conversion.

• length_units (str) – Lenght/distance unit to use for the conversion.

• ureg (Optional[UnitRegistry]) – pint UnitRegistry to use for the conversion. If None is provided, a new UnitRegistry is created.

Return type:

Quantity

Returns:

Conversion factor as a pint.Quantity. conversion_factor.magnitude` gives you the numerical value of the conversion factor.

### Parallel Processing¶

superscreen.parallel.create_models(device, applied_fields=None, circulating_currents=None, vortices=None, layer_updater=None, layer_update_kwargs=None, product=False)[source]

Generate a list of (device, applied_field, circulating_currents).

Parameters:
Return type:
Returns:

A list of “models”, (device, applied_field, circulating_currents, vortices).

superscreen.parallel.solve_many_serial(*, device, applied_fields=None, circulating_currents=None, vortices=None, layer_updater=None, layer_update_kwargs=None, field_units='mT', current_units='uA', check_inversion=False, iterations=1, product=False, directory=None, return_solutions=False, keep_only_final_solution=False, cache_memory_cutoff=inf, log_level=None, gpu=False, use_shared_memory=True, num_cpus=None)[source]

Solve many models in a single process.

Parameters:
Return type:
Returns:

solutions, paths. If return_solutions is True, solutions is either a list of lists of Solutions (if keep_only_final_solution is False), or a list of Solutions (the final iteration for each setup). If directory is True, paths is a list of paths to the saved solutions, otherwise paths is None.

superscreen.parallel.solve_many_mp(*, device, applied_fields=None, circulating_currents=None, vortices=None, layer_updater=None, layer_update_kwargs=None, field_units='mT', current_units='uA', check_inversion=False, iterations=0, product=False, directory=None, return_solutions=False, keep_only_final_solution=False, cache_memory_cutoff=inf, log_level=None, use_shared_memory=True, num_cpus=None)[source]

Solve many models in parallel using multiprocessing.

Parameters:
Return type:
Returns:

solutions, paths. If return_solutions is True, solutions is either a list of lists of Solutions (if keep_only_final_solution is False), or a list of Solutions (the final iteration for each setup). If directory is True, paths is a list of paths to the saved solutions, otherwise paths is None.

superscreen.parallel.solve_many_ray(*, device, applied_fields=None, circulating_currents=None, vortices=None, layer_updater=None, layer_update_kwargs=None, field_units='mT', current_units='uA', check_inversion=False, iterations=0, product=False, directory=None, return_solutions=False, keep_only_final_solution=False, cache_memory_cutoff=inf, log_level=None, use_shared_memory=True, num_cpus=None)[source]

Solve many models in parallel using ray.

Parameters:
Return type:
Returns:

solutions, paths. If return_solutions is True, solutions is either a list of lists of Solutions (if keep_only_final_solution is False), or a list of Solutions (the final iteration for each setup). If directory is True, paths is a list of paths to the saved solutions, otherwise paths is None.