Biplanar coil design

Example showing a basic biplanar coil producing homogeneous field in a target region between the two coil planes.

import numpy as np
import matplotlib.pyplot as plt
from mayavi import mlab
import trimesh


from bfieldtools.mesh_conductor import MeshConductor
from bfieldtools.coil_optimize import optimize_streamfunctions
from bfieldtools.viz import plot_cross_section
from bfieldtools.utils import combine_meshes, load_example_mesh


# Load simple plane mesh that is centered on the origin
planemesh = load_example_mesh("10x10_plane_hires")

# Specify coil plane geometry
center_offset = np.array([0, 0, 0])
standoff = np.array([0, 3, 0])

# Create coil plane pairs
coil_plus = trimesh.Trimesh(
    planemesh.vertices + center_offset + standoff, planemesh.faces, process=False
)

coil_minus = trimesh.Trimesh(
    planemesh.vertices + center_offset - standoff, planemesh.faces, process=False
)

joined_planes = combine_meshes((coil_plus, coil_minus))

# Create mesh class object
coil = MeshConductor(
    mesh_obj=joined_planes, fix_normals=True, basis_name="suh", N_suh=100
)

Out:

Calculating surface harmonics expansion...
Computing the laplacian matrix...
Computing the mass matrix...

Set up target and stray field points

# Here, the target points are on a volumetric grid within a sphere

center = np.array([0, 0, 0])

sidelength = 1.5
n = 8
xx = np.linspace(-sidelength / 2, sidelength / 2, n)
yy = np.linspace(-sidelength / 2, sidelength / 2, n)
zz = np.linspace(-sidelength / 2, sidelength / 2, n)
X, Y, Z = np.meshgrid(xx, yy, zz, indexing="ij")

x = X.ravel()
y = Y.ravel()
z = Z.ravel()

target_points = np.array([x, y, z]).T

# Turn cube into sphere by rejecting points "in the corners"
target_points = (
    target_points[np.linalg.norm(target_points, axis=1) < sidelength / 2] + center
)


# Here, the stray field points are on a spherical surface
stray_radius = 20
stray_points_mesh = trimesh.creation.icosphere(subdivisions=3, radius=stray_radius)
stray_points = stray_points_mesh.vertices + center

n_stray_points = len(stray_points)

Create bfield specifications used when optimizing the coil geometry

# The absolute target field amplitude is not of importance,
# and it is scaled to match the C matrix in the optimization function

target_field = np.zeros(target_points.shape)
target_field[:, 0] += 1

target_spec = {
    "coupling": coil.B_coupling(target_points),
    "abs_error": 0.01,
    "target": target_field,
}
stray_spec = {
    "coupling": coil.B_coupling(stray_points),
    "abs_error": 0.01,
    "target": np.zeros((n_stray_points, 3)),
}

bfield_specification = [target_spec, stray_spec]

Out:

Computing magnetic field coupling matrix, 3184 vertices by 160 target points... took 0.29 seconds.
Computing magnetic field coupling matrix, 3184 vertices by 642 target points... took 0.59 seconds.

# Compute the optimal stream function, either using a numerical solver or regularized least squares

import mosek

coil.s, prob = optimize_streamfunctions(
    coil,
    [target_spec, stray_spec],
    objective="minimum_ohmic_power",
    solver="MOSEK",
    solver_opts={"mosek_params": {mosek.iparam.num_threads: 8}},
)

Out:

Computing the resistance matrix...
Pre-existing problem not passed, creating...
Passing parameters to problem...
Passing problem to solver...


Problem
  Name                   :
  Objective sense        : min
  Type                   : CONIC (conic optimization problem)
  Constraints            : 4914
  Cones                  : 1
  Scalar variables       : 203
  Matrix variables       : 0
  Integer variables      : 0

Optimizer started.
Problem
  Name                   :
  Objective sense        : min
  Type                   : CONIC (conic optimization problem)
  Constraints            : 4914
  Cones                  : 1
  Scalar variables       : 203
  Matrix variables       : 0
  Integer variables      : 0

Optimizer  - threads                : 8
Optimizer  - solved problem         : the dual
Optimizer  - Constraints            : 101
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 2728              conic                  : 102
Optimizer  - Semi-definite variables: 0                 scalarized             : 0
Factor     - setup time             : 0.01              dense det. time        : 0.00
Factor     - ML order time          : 0.00              GP order time          : 0.00
Factor     - nonzeros before factor : 5151              after factor           : 5151
Factor     - dense dim.             : 0                 flops                  : 1.80e+07
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME
0   6.6e+02  1.0e+00  2.0e+00  0.00e+00   0.000000000e+00   -1.000000000e+00  1.0e+00  0.09
1   4.4e+02  6.6e-01  1.6e+00  -9.87e-01  1.242404686e-01   -3.791408582e-01  6.6e-01  0.10
2   1.2e+02  1.8e-01  7.8e-01  -9.60e-01  1.109693224e+00   3.956499435e+00   1.8e-01  0.11
3   1.6e+01  2.5e-02  1.8e-01  -7.07e-01  2.865215552e+01   4.020725058e+01   2.5e-02  0.12
4   7.0e-01  1.1e-03  3.6e-03  2.13e-01   3.463100175e+01   3.738850118e+01   1.1e-03  0.12
5   5.4e-01  8.2e-04  2.3e-03  1.20e+00   2.186623490e+01   2.373566156e+01   8.2e-04  0.13
6   4.7e-01  7.1e-04  1.8e-03  1.16e+00   2.672175433e+01   2.827057607e+01   7.1e-04  0.13
7   3.5e-01  5.3e-04  1.1e-03  1.13e+00   1.305976521e+01   1.412810407e+01   5.3e-04  0.14
8   6.6e-02  1.0e-04  7.6e-05  1.09e+00   1.135207342e+01   1.148805971e+01   1.0e-04  0.14
9   4.3e-03  6.5e-06  7.8e-07  1.02e+00   8.889593586e+00   8.892801886e+00   6.5e-06  0.15
10  1.4e-03  2.2e-06  1.5e-07  1.00e+00   8.380254067e+00   8.381363950e+00   2.2e-06  0.16
11  1.1e-04  1.6e-07  3.2e-09  1.00e+00   8.297805381e+00   8.297892684e+00   1.6e-07  0.16
12  4.1e-06  6.2e-09  2.4e-11  1.00e+00   8.294700570e+00   8.294703920e+00   6.2e-09  0.17
13  2.0e-07  3.0e-10  2.5e-13  1.00e+00   8.294657265e+00   8.294657426e+00   3.0e-10  0.18
14  3.0e-09  4.5e-12  1.2e-15  1.00e+00   8.294654536e+00   8.294654538e+00   4.5e-12  0.19
Optimizer terminated. Time: 0.19


Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: 8.2946545357e+00    nrm: 2e+01    Viol.  con: 1e-10    var: 0e+00    cones: 0e+00
  Dual.    obj: 8.2946545378e+00    nrm: 4e+01    Viol.  con: 7e-06    var: 2e-11    cones: 0e+00

Plot the optimized stream function, then discretize it and plot coil windings and the resultant magnetic field

coil.s.plot()

loops = coil.s.discretize(N_contours=10)

loops.plot_loops()

B_target = loops.magnetic_field(target_points)
mlab.quiver3d(*target_points.T, *B_target.T)
  • ../../_images/sphx_glr_biplanar_coil_design_001.png
  • ../../_images/sphx_glr_biplanar_coil_design_002.png

Out:

<mayavi.modules.vectors.Vectors object at 0x7f0c5979a470>

Lets also do the same coil optimization using regularized least-squares. Now we can’t specify inequality constraints (e.g. use error margins in the specification).

from bfieldtools.coil_optimize import optimize_lsq

coil.s2 = optimize_lsq(
    coil, [target_spec, stray_spec], objective="minimum_ohmic_power", reg=1e6
)

Out:

Error tolerances in specification will be ignored when using lsq

Plot the optimized stream function, then discretize it and plot coil windings and the resultant magnetic field

coil.s2.plot()

loops2 = coil.s2.discretize(N_contours=10)

loops2.plot_loops()

B_target = loops2.magnetic_field(target_points)
mlab.quiver3d(*target_points.T, *B_target.T)
  • ../../_images/sphx_glr_biplanar_coil_design_003.png
  • ../../_images/sphx_glr_biplanar_coil_design_004.png

Out:

<mayavi.modules.vectors.Vectors object at 0x7f0c5926ae30>

Plot cross-section of magnetic field and magnetic potential of the discretized loops

x = y = np.linspace(-12, 12, 250)
X, Y = np.meshgrid(x, y, indexing="ij")


points = np.zeros((X.flatten().shape[0], 3))
points[:, 0] = X.flatten()
points[:, 1] = Y.flatten()

B = loops2.magnetic_field(points)
U = loops2.scalar_potential(points)

U = U.reshape(x.shape[0], y.shape[0])
B = B.T[:2].reshape(2, x.shape[0], y.shape[0])

lw = np.sqrt(B[0] ** 2 + B[1] ** 2)

lw = 2 * lw / np.max(lw)

plot_cross_section(X, Y, U, log=False, contours=False)

seed_points = points[:, :2] * 0.3

plt.streamplot(
    x,
    y,
    B[0],
    B[1],
    density=2,
    linewidth=lw,
    color="k",
    integration_direction="both",
    start_points=seed_points,
)
plt.axis("equal")
plt.axis("off")

plt.plot([-5, 5], [-3, -3], "k", linewidth=3, alpha=1)
plt.plot([-5, 5], [3, 3], "k", linewidth=3, alpha=1)

plt.tight_layout()
../../_images/sphx_glr_biplanar_coil_design_005.png

Total running time of the script: ( 1 minutes 20.879 seconds)

Estimated memory usage: 1494 MB

Gallery generated by Sphinx-Gallery