Note
Go to the end to download the full example code.
Coil with minimal eddy currents¶
Compact example of design of a cylindrical coil surrounded by a RF shield, i.e. a conductive surface. The effects of eddy currents due to inductive interaction with the shield is minimized
import pyvista as pv
from bfieldtools._example_plots import mesh_actor
from bfieldtools._example_plots import vector_glyphs, point_cloud
pv.set_plot_theme("document")
PLOT = True
SAVE_FIGURES = False
SAVE_PATH = "./Minimal eddy current coil/"
import numpy as np
import trimesh
from bfieldtools.mesh_conductor import MeshConductor
from bfieldtools.coil_optimize import optimize_streamfunctions
from bfieldtools.contour import scalar_contour
from bfieldtools.viz import plot_3d_current_loops, plot_data_on_vertices
from bfieldtools import resources
# Set unit, e.g. meter or millimeter.
# This doesn't matter, the problem is scale-invariant
scaling_factor = 1
# Load example coil mesh that is centered on the origin
coilmesh = trimesh.load(
file_obj=resources.files("bfieldtools", "example_meshes/open_cylinder.stl"),
process=True,
)
angle = np.pi / 2
rotation_matrix = np.array(
[
[np.cos(angle), 0, np.sin(angle), 0],
[0, 1, 0, 0],
[-np.sin(angle), 0, np.cos(angle), 0],
[0, 0, 0, 1],
]
)
coilmesh.apply_transform(rotation_matrix)
coilmesh1 = coilmesh.copy()
# coilmesh1.apply_scale(1.3)
coilmesh2 = coilmesh.copy()
# coilmesh1 = coilmesh.union(coilmesh1)
# coilmesh1 = coilmesh1.subdivide().subdivide()
# coilmesh2 = coilmesh.subdivide()
# Create mesh class object
coil = MeshConductor(
verts=coilmesh1.vertices * 0.75,
tris=coilmesh1.faces,
fix_normals=True,
basis_name="suh",
N_suh=400,
)
def alu_sigma(T):
ref_T = 293 # K
ref_rho = 2.82e-8 # ohm*meter
alpha = 0.0039 # 1/K
rho = alpha * (T - ref_T) * ref_rho + ref_rho
return 1 / rho
resistivity = 1 / alu_sigma(T=293) # room-temp Aluminium
thickness = 0.5e-3 # 0.5 mm thick
# Separate object for shield geometry
shield = MeshConductor(
verts=coilmesh2.vertices.copy() * 1.1,
tris=coilmesh2.faces.copy(),
fix_normals=True,
basis_name="inner",
resistivity=resistivity,
thickness=thickness,
)
# shield.mesh.vertices[:,2] -= 3
# shield.mesh.vertices *= np.array([1.2, 1.2, 1.2])
#
# angle = np.pi/2
# rotation_matrix = np.array([[np.cos(angle), 0, np.sin(angle), 0],
# [0, 1, 0, 0],
# [-np.sin(angle), 0, np.cos(angle), 0],
# [0, 0, 0, 1]
# ])
#
# shield.mesh.apply_transform(rotation_matrix)
#
# shield.mesh = shield.mesh.subdivide()
Calculating surface harmonics expansion...
Computing the laplacian matrix...
Computing the mass matrix...
Set up target points and plot geometry
# Here, the target points are on a volumetric grid within a sphere
center = np.array([0, 0, 0])
sidelength = 0.25 * scaling_factor
n = 12
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
)
# Plot coil, shield and target points
if PLOT:
plotter = f = pv.Plotter(window_size=(800, 800))
plotter = coil.plot_mesh(figure=f, opacity=0.2)
plotter = shield.plot_mesh(figure=f, opacity=0.2)
plotter.add_mesh(
point_cloud(*target_points.T).glyph(
scale=False, geom=pv.Sphere(radius=0.05 / 2)
),
show_scalar_bar=False,
)
/home/runner/work/bfieldtools.github.io/bfieldtools.github.io/examples/publication_software/pub_minimal_eddy_current_cylindrical_coil_design.py:146: UserWarning: No vector-like data to use for orient. orient will be set to False.
point_cloud(*target_points.T).glyph(
Compute C matrices that are used to compute the generated magnetic field
mutual_inductance = coil.mutual_inductance(shield)
# Take into account the field produced by currents induced into the shield
# NB! This expression is for instantaneous step-function switching of coil current, see Eq. 18 in G.N. Peeren, 2003.
shield.M_coupling = np.linalg.solve(-shield.inductance, mutual_inductance.T)
secondary_C = shield.B_coupling(target_points) @ -shield.M_coupling
Estimating 69923 MiB required for 4764 by 4764 vertices...
Computing inductance matrix in 120 chunks (13881 MiB memory free), when approx_far=True using more chunks is faster...
Computing triangle-coupling matrix
Computing the inductance matrix...
Computing self-inductance matrix using rough quadrature (degree=2). For higher accuracy, set quad_degree to 4 or more.
Estimating 69923 MiB required for 4764 by 4764 vertices...
Computing inductance matrix in 120 chunks (13661 MiB memory free), when approx_far=True using more chunks is faster...
Computing triangle-coupling matrix
Inductance matrix computation took 46.31 seconds.
Computing magnetic field coupling matrix, 4764 vertices by 672 target points... took 2.16 seconds.
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[:, 1] = target_field[:, 1] + 1
target_spec = {
"coupling": coil.B_coupling(target_points),
"abs_error": 0.01,
"target": target_field,
}
from scipy.linalg import eigh
l, U = eigh(shield.resistance, shield.inductance, subset_by_index=(0, 500))
#
# U = np.zeros((shield.inductance.shape[0], len(li)))
# U[shield.inner_verts, :] = Ui
#
# plt.figure()
# plt.plot(1/li)
# shield.M_coupling = np.linalg.solve(-shield.inductance, mutual_inductance.T)
# secondary_C = shield.B_coupling(target_points) @ -shield.M_coupling
#
# tmin, tmax = 0.001, 0.001
# Fs=10000
time = [0.001, 0.003, 0.005]
eddy_error = [0.05, 0.01, 0.0025]
# time_decay = U @ np.exp(-l[None, :]*time[:, None]) @ np.pinv(U)
time_decay = np.zeros(
(len(time), shield.inductance.shape[0], shield.inductance.shape[1])
)
induction_spec = []
Uinv = np.linalg.pinv(U)
for idx, t in enumerate(time):
time_decay = U @ np.diag(np.exp(-l * t)) @ Uinv
eddy_coupling = shield.B_coupling(target_points) @ time_decay @ shield.M_coupling
induction_spec.append(
{
"coupling": eddy_coupling,
"abs_error": eddy_error[idx],
"rel_error": 0,
"target": np.zeros_like(target_field),
}
)
Computing magnetic field coupling matrix, 4764 vertices by 672 target points... took 2.16 seconds.
Computing the resistance matrix...
Run QP solver
coil.s, prob = optimize_streamfunctions(
coil,
[target_spec] + induction_spec,
objective="minimum_inductive_energy",
solver="CLARABEL",
)
from bfieldtools.mesh_conductor import StreamFunction
shield.induced_s = StreamFunction(shield.M_coupling @ coil.s, shield)
Computing the inductance matrix...
Computing self-inductance matrix using rough quadrature (degree=2). For higher accuracy, set quad_degree to 4 or more.
Estimating 69923 MiB required for 4764 by 4764 vertices...
Computing inductance matrix in 120 chunks (13300 MiB memory free), when approx_far=True using more chunks is faster...
Computing triangle-coupling matrix
Inductance matrix computation took 46.28 seconds.
Passing problem to solver...
===============================================================================
CVXPY
v1.9.2
===============================================================================
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
-------------------------------------------------------------
Clarabel.rs v0.11.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 400
constraints = 16128
nnz(P) = 80200
nnz(A) = 6451200
cones (total) = 1
: Nonnegative = 1, numel = 16128
settings:
linear algebra: direct / faer, precision: 64 bit (4 threads)
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-4, max_scale = 1.0e4
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 +4.6122e+02 -4.1051e+04 9.00e+01 2.23e-01 2.14e-16 1.00e+00 1.42e+01 ------
1 +6.8935e+02 -1.4477e+04 2.20e+01 8.10e-02 2.59e-16 2.24e+00 8.54e+00 6.30e-01
2 +6.9954e+02 -1.4649e+04 2.19e+01 8.14e-02 2.26e-16 2.34e+00 8.52e+00 1.66e-02
3 +8.7439e+02 -1.4682e+04 1.78e+01 7.56e-02 2.70e-16 2.94e+00 7.94e+00 1.85e-01
4 +1.6926e+03 -8.0929e+03 5.78e+00 3.67e-02 3.58e-16 2.84e+00 4.63e+00 5.37e-01
5 +1.9720e+03 -9.9592e+03 6.05e+00 4.07e-02 4.26e-16 3.86e+00 4.11e+00 3.39e-01
6 +9.7392e+03 -3.5326e+03 3.76e+00 1.05e-02 3.54e-16 4.19e+00 1.26e+00 8.73e-01
7 +2.0783e+04 +1.2539e+04 6.58e-01 1.96e-03 3.14e-16 1.51e+00 3.38e-01 9.90e-01
8 +2.1698e+04 +1.4459e+04 5.01e-01 1.68e-03 6.74e-17 1.35e+00 2.84e-01 3.96e-01
9 +2.5003e+04 +2.1079e+04 1.86e-01 7.96e-04 1.57e-16 7.17e-01 1.43e-01 5.85e-01
10 +2.4866e+04 +2.0586e+04 2.08e-01 8.61e-04 2.31e-16 8.29e-01 1.43e-01 2.11e-02
11 +2.6201e+04 +2.2946e+04 1.42e-01 6.22e-04 1.04e-16 6.49e-01 9.70e-02 4.26e-01
12 +2.6810e+04 +2.4010e+04 1.17e-01 5.12e-04 1.26e-16 5.56e-01 7.29e-02 5.61e-01
13 +2.7538e+04 +2.5158e+04 9.46e-02 4.17e-04 1.97e-16 4.69e-01 5.60e-02 4.29e-01
14 +2.9130e+04 +2.7646e+04 5.37e-02 2.39e-04 1.63e-16 2.85e-01 3.10e-02 5.13e-01
15 +2.9164e+04 +2.7724e+04 5.19e-02 2.31e-04 1.92e-16 2.80e-01 2.92e-02 1.54e-01
16 +2.9483e+04 +2.8209e+04 4.52e-02 1.94e-04 7.84e-17 2.54e-01 2.26e-02 3.84e-01
17 +3.0641e+04 +2.9826e+04 2.73e-02 9.77e-05 1.91e-16 1.50e-01 1.05e-02 6.19e-01
18 +3.0803e+04 +3.0139e+04 2.21e-02 7.77e-05 1.67e-16 1.23e-01 7.83e-03 6.24e-01
19 +3.1566e+04 +3.1262e+04 9.74e-03 2.88e-05 1.29e-16 5.11e-02 2.80e-03 6.64e-01
20 +3.2011e+04 +3.1980e+04 9.70e-04 1.47e-06 5.59e-17 3.54e-03 1.81e-04 9.51e-01
21 +3.2032e+04 +3.2032e+04 2.08e-05 2.96e-08 5.06e-16 7.38e-05 3.78e-06 9.80e-01
22 +3.2033e+04 +3.2033e+04 2.19e-07 3.06e-10 1.88e-15 7.87e-07 4.03e-08 9.90e-01
23 +3.2033e+04 +3.2033e+04 3.36e-09 3.18e-11 6.96e-16 2.17e-08 1.03e-09 9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = Solved
solve time = 69.584892241s
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
Plot coil windings and target points
loops = scalar_contour(coil.mesh, coil.s.vert, N_contours=6)
# loops = [simplify_contour(loop, min_edge=1e-2, angle_threshold=2e-2, smooth=True) for loop in loops]
# loops = [loop for loop in loops if loop is not None]
if PLOT:
plotter = f = pv.Plotter(window_size=(600, 500))
plotter.clear()
plotter = plot_3d_current_loops(loops, colors="auto", figure=f, tube_radius=0.005)
B_target = coil.B_coupling(target_points) @ coil.s
plotter.add_mesh(
vector_glyphs(*target_points.T, *B_target.T), show_scalar_bar=False
)
plotter = shield.plot_mesh(
style="surface",
opacity=0.5,
cull_back=True,
color=(0.8, 0.8, 0.8),
figure=f,
)
plotter = shield.plot_mesh(
style="surface",
opacity=1,
cull_front=True,
color=(0.8, 0.8, 0.8),
figure=f,
)
f.camera.parallel_projection = 1
f.camera.zoom(1.4)
if SAVE_FIGURES:
f.screenshot(SAVE_PATH + "eddy_yes.png", scale=4)
plotter.close()
Processing contour, value: -7548711.296065303
Processing contour, value: -4529258.383986609
Processing contour, value: -1509805.4719079155
Processing contour, value: 1509647.440170778
Processing contour, value: 4529100.3522494715
Processing contour, value: 7548553.264328165
Plot cross-section of magnetic field and magnetic potential of the discretized loops
import matplotlib.pyplot as plt
x = y = np.linspace(-1, 1, 40)
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 = coil.B_coupling(points) @ coil.s
plotter.add_mesh(vector_glyphs(*points.T, *B.T), show_scalar_bar=False)
# U = U.reshape(x.shape[0], y.shape[0])
# B = B.T[:2].reshape(2, x.shape[0], y.shape[0])
# from bfieldtools.viz import plot_cross_section
# 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.tight_layout()
Computing magnetic field coupling matrix, 4764 vertices by 1600 target points... took 5.08 seconds.
Actor (0x7f12701be8c0)
Center: (-1.7881393432617188e-07, 4.231929779052734e-06, -1.1791475117206573e-05)
Pickable: True
Position: (0.0, 0.0, 0.0)
Scale: (1.0, 1.0, 1.0)
Visible: True
X Bounds -1.001E+00, 1.001E+00
Y Bounds -1.000E+00, 1.000E+00
Z Bounds -1.425E-02, 1.423E-02
User matrix: Identity
Has mapper: True
Property (0x7f12997a8a00)
Ambient: 0.0
Ambient color: Color(name='light_blue', hex='#add8e6ff', opacity=255)
Anisotropy: 0.0
Color: Color(name='light_blue', hex='#add8e6ff', opacity=255)
Culling: "none"
Diffuse: 1.0
Diffuse color: Color(name='light_blue', hex='#add8e6ff', opacity=255)
Edge color: Color(name='black', hex='#000000ff', opacity=255)
Edge opacity: 1.0
Interpolation: InterpolationType.FLAT
Lighting: True
Line width: 1.0
Metallic: 0.0
Opacity: 1.0
Point size: 5.0
Render lines as tubes: False
Render points as spheres: False
Roughness: 0.5
Show edges: False
Specular: 0.0
Specular color: Color(name='light_blue', hex='#add8e6ff', opacity=255)
Specular power: 100.0
Style: "Surface"
DataSetMapper (0x7f1299618400)
Scalar visibility: True
Scalar range: (0.028333157300949097, 80.52881622314453)
Interpolate before mapping: True
Scalar map mode: point
Color mode: map
Attached dataset:
PolyData (0x7f1278e810c0)
N Cells: 68800
N Points: 161600
N Strips: 0
X Bounds: -1.001e+00, 1.001e+00
Y Bounds: -1.000e+00, 1.000e+00
Z Bounds: -1.425e-02, 1.423e-02
N Arrays: 2
For comparison, let’s see how the coils look when we ignore the conducting shield
coil.unshielded_s, coil.unshielded_prob = optimize_streamfunctions(
coil,
[target_spec],
objective="minimum_inductive_energy",
solver="CLARABEL",
)
shield.unshielded_induced_s = StreamFunction(
shield.M_coupling @ coil.unshielded_s, shield
)
loops = scalar_contour(coil.mesh, coil.unshielded_s.vert, N_contours=6)
if PLOT:
plotter = f = pv.Plotter(window_size=(600, 500))
plotter.clear()
plotter = plot_3d_current_loops(loops, colors="auto", figure=f, tube_radius=0.005)
B_target_unshielded = coil.B_coupling(target_points) @ coil.unshielded_s
plotter.add_mesh(
vector_glyphs(*target_points.T, *B_target_unshielded.T), show_scalar_bar=False
)
plotter = shield.plot_mesh(
style="surface",
opacity=0.5,
cull_back=True,
color=(0.8, 0.8, 0.8),
figure=f,
)
plotter = shield.plot_mesh(
style="surface",
opacity=1,
cull_front=True,
color=(0.8, 0.8, 0.8),
figure=f,
)
f.camera.parallel_projection = 1
f.camera.zoom(1.4)
if SAVE_FIGURES:
f.screenshot(SAVE_PATH + "eddy_no.png", scale=4)
plotter.close()
Passing problem to solver...
===============================================================================
CVXPY
v1.9.2
===============================================================================
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
-------------------------------------------------------------
Clarabel.rs v0.11.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 400
constraints = 4032
nnz(P) = 80200
nnz(A) = 1612800
cones (total) = 1
: Nonnegative = 1, numel = 4032
settings:
linear algebra: direct / faer, precision: 64 bit (4 threads)
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-4, max_scale = 1.0e4
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 +4.6056e+01 -3.8110e+02 9.27e+00 1.72e-01 1.01e-15 1.00e+00 1.27e+00 ------
1 +4.5453e+01 -4.0183e+01 2.13e+00 4.10e-02 1.59e-15 2.18e+00 9.50e-01 9.90e-01
2 +4.4513e+01 +2.0872e+01 1.13e+00 1.29e-02 1.09e-15 3.03e+00 5.50e-01 5.84e-01
3 +4.4707e+01 +2.1420e+01 1.09e+00 1.27e-02 7.95e-16 2.97e+00 5.43e-01 5.21e-02
4 +4.5740e+01 +2.8456e+01 6.07e-01 9.27e-03 3.86e-16 2.22e+00 4.40e-01 3.36e-01
5 +4.6480e+01 +4.2844e+01 8.49e-02 1.99e-03 1.93e-16 5.91e-01 1.23e-01 7.63e-01
6 +4.6534e+01 +4.3760e+01 6.34e-02 1.43e-03 1.69e-16 2.66e-01 9.01e-02 5.75e-01
7 +4.6687e+01 +4.6432e+01 5.50e-03 1.27e-04 1.34e-16 1.61e-02 8.51e-03 9.42e-01
8 +4.6700e+01 +4.6695e+01 9.65e-05 2.06e-06 1.73e-16 2.77e-04 1.49e-04 9.84e-01
9 +4.6700e+01 +4.6700e+01 9.80e-07 2.06e-08 2.18e-16 2.81e-06 1.52e-06 9.90e-01
10 +4.6700e+01 +4.6700e+01 9.80e-09 2.06e-10 1.45e-16 2.81e-08 1.52e-08 9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = Solved
solve time = 7.20549645s
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
Processing contour, value: -520119.67162245914
Processing contour, value: -312073.776620582
Processing contour, value: -104027.88161870476
Processing contour, value: 104018.01338317245
Processing contour, value: 312063.90838504955
Processing contour, value: 520109.80338692677
import matplotlib.pyplot as plt
x = y = np.linspace(-1, 1, 40)
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 = coil.B_coupling(points) @ coil.unshielded_s
plotter.add_mesh(vector_glyphs(*points.T, *B.T), show_scalar_bar=False)
Actor (0x7f12995ecca0)
Center: (5.960464477539063e-08, 0.004965066909790039, 8.335337042808533e-07)
Pickable: True
Position: (0.0, 0.0, 0.0)
Scale: (1.0, 1.0, 1.0)
Visible: True
X Bounds -1.005E+00, 1.005E+00
Y Bounds -1.000E+00, 1.010E+00
Z Bounds -1.415E-02, 1.415E-02
User matrix: Identity
Has mapper: True
Property (0x7f1299822f20)
Ambient: 0.0
Ambient color: Color(name='light_blue', hex='#add8e6ff', opacity=255)
Anisotropy: 0.0
Color: Color(name='light_blue', hex='#add8e6ff', opacity=255)
Culling: "none"
Diffuse: 1.0
Diffuse color: Color(name='light_blue', hex='#add8e6ff', opacity=255)
Edge color: Color(name='black', hex='#000000ff', opacity=255)
Edge opacity: 1.0
Interpolation: InterpolationType.FLAT
Lighting: True
Line width: 1.0
Metallic: 0.0
Opacity: 1.0
Point size: 5.0
Render lines as tubes: False
Render points as spheres: False
Roughness: 0.5
Show edges: False
Specular: 0.0
Specular color: Color(name='light_blue', hex='#add8e6ff', opacity=255)
Specular power: 100.0
Style: "Surface"
DataSetMapper (0x7f127984ebc0)
Scalar visibility: True
Scalar range: (0.010458321310579777, 1.6912238597869873)
Interpolate before mapping: True
Scalar map mode: point
Color mode: map
Attached dataset:
PolyData (0x7f12998219c0)
N Cells: 68800
N Points: 161600
N Strips: 0
X Bounds: -1.005e+00, 1.005e+00
Y Bounds: -1.000e+00, 1.010e+00
Z Bounds: -1.415e-02, 1.415e-02
N Arrays: 2
Finally, let’s compare the time-courses
tmin, tmax = 0, 0.025
Fs = 2000
time = np.linspace(tmin, tmax, int(Fs * (tmax - tmin) + 1))
# time_decay = U @ np.exp(-l[None, :]*time[:, None]) @ np.pinv(U)
time_decay = np.zeros(
(len(time), shield.inductance.shape[0], shield.inductance.shape[1])
)
Uinv = np.linalg.pinv(U)
for idx, t in enumerate(time):
time_decay[idx] = U @ np.diag(np.exp(-l * t)) @ Uinv
B_t = shield.B_coupling(target_points) @ (time_decay @ shield.induced_s).T
unshieldedB_t = (
shield.B_coupling(target_points) @ (time_decay @ shield.unshielded_induced_s).T
)
import matplotlib.pyplot as plt
if PLOT and SAVE_FIGURES:
fig, ax = plt.subplots(1, 1, sharex=True, figsize=(8, 4))
ax.plot(
time * 1e3,
np.mean(np.linalg.norm(B_t, axis=1), axis=0).T,
"k-",
label="Constrained",
linewidth=1.5,
)
# ax[0].set_title('Eddy currents minimized')
ax.set_ylabel("Transient field amplitude")
ax.semilogy(
time * 1e3,
np.mean(np.linalg.norm(unshieldedB_t, axis=1), axis=0).T,
"k--",
label="Ignored",
linewidth=1.5,
)
# ax[1].set_title('Eddy currents ignored')
ax.set_xlabel("Time (ms)")
# ax[1].set_ylabel('Transient field amplitude')
ax.set_ylim(1e-4, 0.5)
ax.set_xlim(0, 25)
#
# ax.spines['top'].set_visible(False)
# ax.spines['right'].set_visible(False)
plt.grid(which="both", axis="y", alpha=0.1)
plt.legend()
fig.tight_layout()
ax.vlines([1, 5, 10, 20], 1e-4, 0.5, alpha=0.1, linewidth=3, color="r")
plt.savefig(SAVE_PATH + "eddy_transient.pdf")
from bfieldtools.mesh_calculus import gradient
if PLOT and SAVE_FIGURES:
for plot_time_idx in [2, 10, 20, 40]:
plotter = f = pv.Plotter(window_size=(600, 500))
plotter.clear()
plotter = shield.plot_mesh(
style="surface", color=(0.8, 0.8, 0.8), opacity=1, figure=f
)
s = np.zeros((shield.mesh.vertices.shape[0],))
s[shield.inner_vertices] = time_decay[plot_time_idx] @ shield.induced_s
plotter = plot_data_on_vertices(
shield.mesh,
s,
ncolors=256,
figure=f,
opacity=1,
cull_back=False,
colormap="RdBu",
)
surface1 = mesh_actor(plotter)
surface1.mapper.dataset = surface1.mapper.dataset.contour(10)
surface1.prop.line_width = 10.0
f.camera.parallel_projection = 1
f.view_isometric()
f.camera.zoom(1.4)
f.screenshot(
SAVE_PATH + "shield_eddy_yes_time_%.3f.png" % time[plot_time_idx], scale=2
)
plotter.close()
for plot_time_idx in [2, 10, 20, 40]:
plotter = f = pv.Plotter(window_size=(600, 500))
plotter = shield.plot_mesh(
style="surface", color=(0.8, 0.8, 0.8), opacity=1, figure=f
)
s_u = np.zeros((shield.mesh.vertices.shape[0],))
s_u[shield.inner_vertices] = (
time_decay[plot_time_idx] @ shield.unshielded_induced_s
)
plotter = plot_data_on_vertices(
shield.mesh,
s_u,
ncolors=256,
figure=f,
opacity=1,
cull_back=False,
colormap="RdBu",
)
surface1 = mesh_actor(plotter)
surface1.mapper.dataset = surface1.mapper.dataset.contour(10)
surface1.prop.line_width = 10.0
f.camera.parallel_projection = 1
f.view_isometric()
f.camera.zoom(1.4)
f.screenshot(
SAVE_PATH + "shield_eddy_no_time_%.3f.png" % time[plot_time_idx], scale=2
)
plotter.close()
# Display the assembled scene (Sphinx-Gallery captures it during docs builds).
if "plotter" in globals() and not plotter._closed:
plotter.show()
Total running time of the script: (5 minutes 10.292 seconds)