Skip to content

Commit

Permalink
Merge pull request #200 from hiddenSymmetries/ScalarPotentialRZMagnet…
Browse files Browse the repository at this point in the history 
…icField_fix

Fix to ScalarPotentialRZMagneticField
  • Loading branch information
@ejpaul
ejpaul authored Feb 22, 2022
2 parents cb377be + 07c24f0 commit c5702a4
Show file tree
Hide file tree
Showing 2 changed files with 96 additions and 17 deletions.
46 changes: 39 additions & 7 deletions src/simsopt/field/magneticfieldclasses.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,15 +245,47 @@ def _B_impl(self, B):
r = np.sqrt(np.square(points[:, 0]) + np.square(points[:, 1]))
z = points[:, 2]
phi = np.arctan2(points[:, 1], points[:, 0])
B[:] = np.array(self.Blambdify(r, z, phi)).T
B_cyl = np.array(self.Blambdify(r, z, phi)).T
# Bx = Br cos(phi) - Bphi sin(phi)
B[:, 0] = B_cyl[:, 0] * np.cos(phi) - B_cyl[:, 1] * np.sin(phi)
# By = Br sin(phi) + Bphi cos(phi)
B[:, 1] = B_cyl[:, 0] * np.sin(phi) + B_cyl[:, 1] * np.cos(phi)
B[:, 2] = B_cyl[:, 2]

def _dB_by_dX_impl(self, dB):
points = self.get_points_cart_ref()
r = np.sqrt(np.square(points[:, 0]) + np.square(points[:, 1]))
r = np.sqrt(np.square(points[:, 0]) + np.square(points[:, 1]))
z = points[:, 2]
phi = np.arctan2(points[:, 1], points[:, 0])
dB[:] = np.array(self.dBlambdify_by_dX(r, z, phi)).transpose((2, 0, 1))
dB_cyl = np.array(self.dBlambdify_by_dX(r, z, phi)).transpose((2, 0, 1))
dBrdx = dB_cyl[:, 0, 0]
dBrdy = dB_cyl[:, 1, 0]
dBrdz = dB_cyl[:, 2, 0]
dBphidx = dB_cyl[:, 0, 1]
dBphidy = dB_cyl[:, 1, 1]
dBphidz = dB_cyl[:, 2, 1]
dB[:, 0, 2] = dB_cyl[:, 0, 2]
dB[:, 1, 2] = dB_cyl[:, 1, 2]
dB[:, 2, 2] = dB_cyl[:, 2, 2]
dcosphidx = -points[:, 0]**2/r**3 + 1/r
dsinphidx = -points[:, 0]*points[:, 1]/r**3
dcosphidy = -points[:, 0]*points[:, 1]/r**3
dsinphidy = -points[:, 1]**2/r**3 + 1/r
B_cyl = np.array(self.Blambdify(r, z, phi)).T
Br = B_cyl[:, 0]
Bphi = B_cyl[:, 1]
# Bx = Br cos(phi) - Bphi sin(phi)
dB[:, 0, 0] = dBrdx * np.cos(phi) + Br * dcosphidx - dBphidx * np.sin(phi) \
- Bphi * dsinphidx
dB[:, 1, 0] = dBrdy * np.cos(phi) + Br * dcosphidy - dBphidy * np.sin(phi) \
- Bphi * dsinphidy
dB[:, 2, 0] = dBrdz * np.cos(phi) - dBphidz * np.sin(phi)
# By = Br sin(phi) + Bphi cos(phi)
dB[:, 0, 1] = dBrdx * np.sin(phi) + Br * dsinphidx + dBphidx * np.cos(phi) \
+ Bphi * dcosphidx
dB[:, 1, 1] = dBrdy * np.sin(phi) + Br * dsinphidy + dBphidy * np.cos(phi) \
+ Bphi * dcosphidy
dB[:, 2, 1] = dBrdz * np.sin(phi) + dBphidz * np.cos(phi)


class CircularCoil(MagneticField):
Expand Down Expand Up @@ -292,13 +324,13 @@ def __init__(self, r0=0.1, center=[0, 0, 0], I=5e5/np.pi, normal=[0, 0]):

self.rotMatrix = np.array([
[np.cos(phi) * np.cos(theta)**2 + np.sin(theta)**2,
-np.sin(phi / 2)**2 * np.sin(2 * theta),
-np.sin(phi / 2)**2 * np.sin(2 * theta),
np.cos(theta) * np.sin(phi)],
[-np.sin(phi / 2)**2 * np.sin(2 * theta),
np.cos(theta)**2 + np.cos(phi) * np.sin(theta)**2,
[-np.sin(phi / 2)**2 * np.sin(2 * theta),
np.cos(theta)**2 + np.cos(phi) * np.sin(theta)**2,
np.sin(phi) * np.sin(theta)],
[-np.cos(theta) * np.sin(phi),
-np.sin(phi) * np.sin(theta),
-np.sin(phi) * np.sin(theta),
np.cos(phi)]
])

Expand Down
67 changes: 57 additions & 10 deletions tests/field/test_magneticfields.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,16 +116,63 @@ def test_scalarpotential_Bfield(self):
Bscalar.set_points(points)
B1 = np.array(Bscalar.B())
dB1_by_dX = np.array(Bscalar.dB_by_dX())

# Analytical Formula for B
rphiz = [[np.sqrt(np.power(point[0], 2) + np.power(point[1], 2)), np.arctan2(point[1], point[0]), point[2]] for point in points]
B2 = np.array([[0.2*point[2]+0.8*point[0], (0.1+0.3*point[2])/point[0], 0.2*point[0]+0.3*point[1]+point[2]] for point in rphiz])
# Convert to Cartesian coordinates
r = np.sqrt(np.power(points[:, 0], 2) + np.power(points[:, 1], 2))
phi = np.arctan2(points[:, 1], points[:, 0])
z = points[:, 2]
B2_cart = np.zeros_like(B2)
# Bx = Br cos(phi) - Bphi sin(phi)
B2_cart[:, 0] = B2[:, 0] * np.cos(phi) - B2[:, 1] * np.sin(phi)
# By = Br sin(phi) + Bphi cos(phi)
B2_cart[:, 1] = B2[:, 0] * np.sin(phi) + B2[:, 1] * np.cos(phi)
B2_cart[:, 2] = B2[:, 2]
dB2_by_dX = np.array([
[[0.8*np.cos(point[1]), -(np.cos(point[1])/point[0]**2)*(0.1+0.3*point[2]), 0.2*np.cos(point[1])-0.3*np.sin(point[1])/point[0]],
[0.8*np.sin(point[1]), -(np.sin(point[1])/point[0]**2)*(0.1+0.3*point[2]), 0.2*np.sin(point[1])+0.3*np.cos(point[1])/point[0]],
[0.2, 0.3/point[0], 1]] for point in rphiz])
dBxdx = dB1_by_dX[:, 0, 0]
dBxdy = dB1_by_dX[:, 1, 0]
dBxdz = dB1_by_dX[:, 2, 0]
dBydx = dB1_by_dX[:, 0, 1]
dBydy = dB1_by_dX[:, 1, 1]
dBydz = dB1_by_dX[:, 2, 1]
dB1_by_dX_cyl = np.zeros_like(dB2_by_dX)
dcosphidx = -points[:, 0]**2/r**3 + 1/r
dsinphidx = -points[:, 0]*points[:, 1]/r**3
dcosphidy = -points[:, 0]*points[:, 1]/r**3
dsinphidy = -points[:, 1]**2/r**3 + 1/r
Bx = B1[:, 0]
By = B1[:, 1]
# Br = Bx cos(phi) + By sin(phi)
dB1_by_dX_cyl[:, 0, 0] = dBxdx * np.cos(phi) + Bx * dcosphidx + dBydx * np.sin(phi) \
+ By * dsinphidx
dB1_by_dX_cyl[:, 1, 0] = dBxdy * np.cos(phi) + Bx * dcosphidy + dBydy * np.sin(phi) \
+ By * dsinphidy
dB1_by_dX_cyl[:, 2, 0] = dBxdz * np.cos(phi) + dBydz * np.sin(phi)
# Bphi = - sin(phi) Bx + cos(phi) By
dB1_by_dX_cyl[:, 0, 1] = - dBxdx * np.sin(phi) - Bx * dsinphidx + dBydx * np.cos(phi) \
+ By * dcosphidx
dB1_by_dX_cyl[:, 1, 1] = - dBxdy * np.sin(phi) - Bx * dsinphidy + dBydy * np.cos(phi) \
+ By * dcosphidy
dB1_by_dX_cyl[:, 2, 1] = - dBxdz * np.sin(phi) + dBydz * np.cos(phi)
dB1_by_dX_cyl[:, :, 2] = dB1_by_dX[:, :, 2]
# Verify
assert np.allclose(B1, B2)
assert np.allclose(dB1_by_dX, dB2_by_dX)
assert np.allclose(B1, B2_cart)
assert np.allclose(dB1_by_dX_cyl, dB2_by_dX)

# Check for divergence-free condition for dipole field
# Set up magnetic field scalar potential
PhiStr = "Z/(R*R + Z*Z)**(3/2)"
# Set up scalar potential B
Bscalar = ScalarPotentialRZMagneticField(PhiStr)
Bscalar.set_points(points)
dB1_by_dX = np.array(Bscalar.dB_by_dX())
divB = dB1_by_dX[:, 0, 0] + dB1_by_dX[:, 1, 1] + dB1_by_dX[:, 2, 2]
assert np.allclose(np.abs(divB), 0)

def test_circularcoil_Bfield(self):
current = 1.2e7
Expand Down Expand Up @@ -284,8 +331,8 @@ def test_circularcoil_Bfield_toroidal_arrangement(self):
N_coils = 30

N_turns = 3
a1 = 10 / 2 * 0.0254
a2 = 19.983 / 2 * 0.0254
a1 = 10 / 2 * 0.0254
a2 = 19.983 / 2 * 0.0254
r_array = np.linspace(a1, a2, N_turns)
I_amp = 433 * (33/N_turns)

Expand Down Expand Up @@ -420,16 +467,16 @@ def test_Reiman(self):
x = points[:, 0]
y = points[:, 1]
z = points[:, 2]
Bx = (y*np.sqrt(x**2 + y**2) + x*z*(0.15 + 0.38*((-1 + np.sqrt(x**2 + y**2))**2 + z**2) -
0.06*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.cos(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2))))) +
Bx = (y*np.sqrt(x**2 + y**2) + x*z*(0.15 + 0.38*((-1 + np.sqrt(x**2 + y**2))**2 + z**2) -
0.06*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.cos(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2))))) +
0.06*x*(1 - np.sqrt(x**2 + y**2))*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2 *
np.sin(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2)))))/(x**2 + y**2)
By = (-1.*x*np.sqrt(x**2 + y**2) + y*z*(0.15 + 0.38*((-1 + np.sqrt(x**2 + y**2))**2 + z**2) -
0.06*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.cos(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2))))) +
By = (-1.*x*np.sqrt(x**2 + y**2) + y*z*(0.15 + 0.38*((-1 + np.sqrt(x**2 + y**2))**2 + z**2) -
0.06*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.cos(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2))))) +
0.06*y*(1 - np.sqrt(x**2 + y**2))*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2 *
np.sin(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2)))))/(x**2 + y**2)
Bz = (-((-1 + np.sqrt(x**2 + y**2))*(0.15 + 0.38*((-1 + np.sqrt(x**2 + y**2))**2 + z**2) -
0.06*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.cos(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2)))))) -
Bz = (-((-1 + np.sqrt(x**2 + y**2))*(0.15 + 0.38*((-1 + np.sqrt(x**2 + y**2))**2 + z**2) -
0.06*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.cos(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2)))))) -
0.06*z*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.sin(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2)))))/np.sqrt(x**2 + y**2)
B2 = np.array(np.vstack((Bx, By, Bz)).T)
assert np.allclose(B1, B2)
Expand Down

0 comments on commit c5702a4

Please sign in to comment.