Skip to content

ag/gradgrad b fix #24

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
andrewgiuliani wants to merge 5 commits into
base: main
Choose a base branch
from
Open

ag/gradgrad b fix #24

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5e298b3
I think that this fixes the grad grad B basis conversion
andrewgiuliani May 29, 2025
e6bbb36
added unit test for B(phi), B'(phi), B''(phi) in the axisymmetric case
andrewgiuliani May 29, 2025
cedac03
fixed Frenet to cylindrical and cylindrical to Cartesian tensor conve…
andrewgiuliani May 29, 2025
6c1a703
added a docstring to the unit test
andrewgiuliani May 29, 2025
b25274f
atol=rtol=1e-8
andrewgiuliani May 29, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
110 changes: 33 additions & 77 deletions qsc/grad_B_tensor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1306,85 +1306,41 @@ def grad_B_tensor_cartesian(self):
return grad_B_vector_cartesian

def grad_grad_B_tensor_cylindrical(self):
'''
Function to calculate the gradient of of the gradient the magnetic field
vector B=(B_R,B_phi,B_Z) at every point along the axis (hence with nphi points)
where R, phi and Z are the standard cylindrical coordinates.
'''
return np.transpose(self.grad_grad_B,(1,2,3,0))
'''
Function to calculate the gradient of of the gradient the magnetic field
vector B=(B_R,B_phi,B_Z) at every point along the axis (hence with nphi points)
where R, phi and Z are the standard cylindrical coordinates.
'''
# return np.transpose(self.grad_grad_B,(1,2,3,0))
t = self.tangent_cylindrical
n = self.normal_cylindrical
b = self.binormal_cylindrical
E = np.stack([n, b, t], axis=0)
grad_grad_B = self.grad_grad_B.transpose(1, 2, 3, 0)
grad_grad_B_tensor_cylindrical = np.einsum('ijkq,iqa,jqb,kqc->abcq', grad_grad_B, E, E, E)
return grad_grad_B_tensor_cylindrical

def grad_grad_B_tensor_cartesian(self):
'''
Function to calculate the gradient of of the gradient the magnetic field
vector B=(B_x,B_y,B_z) at every point along the axis (hence with nphi points)
where x, y and z are the standard cartesian coordinates.
'''
nablanablaB = self.grad_grad_B_tensor_cylindrical()
"""
Transform the second derivative tensor of the magnetic field from cylindrical to Cartesian coordinates.
Shape: (nphi, 3, 3, 3)
"""
grad_grad_B_cyl = self.grad_grad_B_tensor_cylindrical().transpose(3, 0, 1, 2) # shape (nphi, 3, 3, 3)
cosphi = np.cos(self.phi)
sinphi = np.sin(self.phi)

grad_grad_B_vector_cartesian = np.array([[
[cosphi**3*nablanablaB[0, 0, 0] - cosphi**2*sinphi*(nablanablaB[0, 0, 1] +
nablanablaB[0, 1, 0] + nablanablaB[1, 0, 0]) +
cosphi*sinphi**2*(nablanablaB[0, 1, 1] + nablanablaB[1, 0, 1] +
nablanablaB[1, 1, 0]) - sinphi**3*nablanablaB[1, 1, 1],
cosphi**3*nablanablaB[0, 0, 1] + cosphi**2*sinphi*(nablanablaB[0, 0, 0] -
nablanablaB[0, 1, 1] - nablanablaB[1, 0, 1]) + sinphi**3*nablanablaB[1, 1, 0] -
cosphi*sinphi**2*(nablanablaB[0, 1, 0] + nablanablaB[1, 0, 0] -
nablanablaB[1, 1, 1]), cosphi**2*nablanablaB[0, 0, 2] -
cosphi*sinphi*(nablanablaB[0, 1, 2] + nablanablaB[1, 0, 2]) +
sinphi**2*nablanablaB[1, 1, 2]], [cosphi**3*nablanablaB[0, 1, 0] +
sinphi**3*nablanablaB[1, 0, 1] + cosphi**2*sinphi*(nablanablaB[0, 0, 0] -
nablanablaB[0, 1, 1] - nablanablaB[1, 1, 0]) -
cosphi*sinphi**2*(nablanablaB[0, 0, 1] + nablanablaB[1, 0, 0] -
nablanablaB[1, 1, 1]), cosphi**3*nablanablaB[0, 1, 1] -
sinphi**3*nablanablaB[1, 0, 0] + cosphi*sinphi**2*(nablanablaB[0, 0, 0] -
nablanablaB[1, 0, 1] - nablanablaB[1, 1, 0]) +
cosphi**2*sinphi*(nablanablaB[0, 0, 1] + nablanablaB[0, 1, 0] -
nablanablaB[1, 1, 1]), cosphi**2*nablanablaB[0, 1, 2] -
sinphi**2*nablanablaB[1, 0, 2] + cosphi*sinphi*(nablanablaB[0, 0, 2] -
nablanablaB[1, 1, 2])], [cosphi**2*nablanablaB[0, 2, 0] -
cosphi*sinphi*(nablanablaB[0, 2, 1] + nablanablaB[1, 2, 0]) +
sinphi**2*nablanablaB[1, 2, 1], cosphi**2*nablanablaB[0, 2, 1] -
sinphi**2*nablanablaB[1, 2, 0] + cosphi*sinphi*(nablanablaB[0, 2, 0] -
nablanablaB[1, 2, 1]), cosphi*nablanablaB[0, 2, 2] -
sinphi*nablanablaB[1, 2, 2]]],
[[sinphi**3*nablanablaB[0, 1, 1] + cosphi**3*nablanablaB[1, 0, 0] +
cosphi**2*sinphi*(nablanablaB[0, 0, 0] - nablanablaB[1, 0, 1] -
nablanablaB[1, 1, 0]) - cosphi*sinphi**2*(nablanablaB[0, 0, 1] +
nablanablaB[0, 1, 0] - nablanablaB[1, 1, 1]), -(sinphi**3*nablanablaB[0, 1, 0]) +
cosphi**3*nablanablaB[1, 0, 1] + cosphi*sinphi**2*(nablanablaB[0, 0, 0] -
nablanablaB[0, 1, 1] - nablanablaB[1, 1, 0]) +
cosphi**2*sinphi*(nablanablaB[0, 0, 1] + nablanablaB[1, 0, 0] -
nablanablaB[1, 1, 1]), -(sinphi**2*nablanablaB[0, 1, 2]) +
cosphi**2*nablanablaB[1, 0, 2] + cosphi*sinphi*(nablanablaB[0, 0, 2] -
nablanablaB[1, 1, 2])], [-(sinphi**3*nablanablaB[0, 0, 1]) +
cosphi*sinphi**2*(nablanablaB[0, 0, 0] - nablanablaB[0, 1, 1] -
nablanablaB[1, 0, 1]) + cosphi**3*nablanablaB[1, 1, 0] +
cosphi**2*sinphi*(nablanablaB[0, 1, 0] + nablanablaB[1, 0, 0] -
nablanablaB[1, 1, 1]), sinphi**3*nablanablaB[0, 0, 0] +
cosphi*sinphi**2*(nablanablaB[0, 0, 1] + nablanablaB[0, 1, 0] +
nablanablaB[1, 0, 0]) + cosphi**2*sinphi*(nablanablaB[0, 1, 1] +
nablanablaB[1, 0, 1] + nablanablaB[1, 1, 0]) + cosphi**3*nablanablaB[1, 1, 1],
sinphi**2*nablanablaB[0, 0, 2] + cosphi*sinphi*(nablanablaB[0, 1, 2] +
nablanablaB[1, 0, 2]) + cosphi**2*nablanablaB[1, 1, 2]],
[-(sinphi**2*nablanablaB[0, 2, 1]) + cosphi**2*nablanablaB[1, 2, 0] +
cosphi*sinphi*(nablanablaB[0, 2, 0] - nablanablaB[1, 2, 1]),
sinphi**2*nablanablaB[0, 2, 0] + cosphi*sinphi*(nablanablaB[0, 2, 1] +
nablanablaB[1, 2, 0]) + cosphi**2*nablanablaB[1, 2, 1],
sinphi*nablanablaB[0, 2, 2] + cosphi*nablanablaB[1, 2, 2]]],
[[cosphi**2*nablanablaB[2, 0, 0] - cosphi*sinphi*(nablanablaB[2, 0, 1] +
nablanablaB[2, 1, 0]) + sinphi**2*nablanablaB[2, 1, 1],
cosphi**2*nablanablaB[2, 0, 1] - sinphi**2*nablanablaB[2, 1, 0] +
cosphi*sinphi*(nablanablaB[2, 0, 0] - nablanablaB[2, 1, 1]),
cosphi*nablanablaB[2, 0, 2] - sinphi*nablanablaB[2, 1, 2]],
[-(sinphi**2*nablanablaB[2, 0, 1]) + cosphi**2*nablanablaB[2, 1, 0] +
cosphi*sinphi*(nablanablaB[2, 0, 0] - nablanablaB[2, 1, 1]),
sinphi**2*nablanablaB[2, 0, 0] + cosphi*sinphi*(nablanablaB[2, 0, 1] +
nablanablaB[2, 1, 0]) + cosphi**2*nablanablaB[2, 1, 1],
sinphi*nablanablaB[2, 0, 2] + cosphi*nablanablaB[2, 1, 2]],
[cosphi*nablanablaB[2, 2, 0] - sinphi*nablanablaB[2, 2, 1],
sinphi*nablanablaB[2, 2, 0] + cosphi*nablanablaB[2, 2, 1], nablanablaB[2, 2, 2]]
]])

return grad_grad_B_vector_cartesian
# Rotation matrix R from cylindrical (R, φ, Z) to Cartesian (x, y, z)
R = np.array([
[cosphi, -sinphi, np.zeros_like(cosphi)],
[sinphi, cosphi, np.zeros_like(cosphi)],
[np.zeros_like(cosphi), np.zeros_like(cosphi), np.ones_like(cosphi)],
]) # shape (3, 3, nphi)

# Transpose to shape (nphi, 3, 3)
R = np.transpose(R, (2, 0, 1))

# Apply tensor transformation: R[i,a] * R[j,b] * R[k,c] * T[a,b,c]
T_cartesian = np.einsum('pia,pjb,pkc,pabc->pijk',
R, R, R, grad_grad_B_cyl)

return T_cartesian.transpose(1,2,3,0) # shape (nphi, 3, 3, 3)
84 changes: 84 additions & 0 deletions qsc/tests/test_grad_B_tensor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,6 +121,90 @@ def test_axisymmetry(self):
np.testing.assert_allclose(s.grad_grad_B_alt[:,0,2,0], np.full(nphi, val), rtol=rtol, atol=atol)
np.testing.assert_allclose(s.grad_grad_B_alt[:,0,0,2], np.full(nphi, val), rtol=rtol, atol=atol)

class DirectionalDerivativeTests(unittest.TestCase):
"""
This test computes the B(\phi), B'(\phi), and B''(\phi) on the magnetic
axis in two ways, where \phi is the standard cylindrical angle.

The first way uses the fact that B(\phi) = sG*B0*T(\phi), where T(\phi) is
the tangent vector of the Frenet frame associated to the magnetic axis. Since
we have analytical expressions for the Frenet frame, we can differentiate them
with respect to \phi.

The second way relies on the gradB, and gradgradB tensors in pyqsc associated
to the NAE. We know that hatB(\phi) = B(\Gamma(\phi)) where Bhat is the magnetic
field B evaluated on the magnetic axis, given by \Gamma(\phi). Differentiating
this expression with respect to \phi and applying the chain rule, we obtain
a formula for hatB(\phi) that depends on derivative tensors of B.

With these two independent formulas for B, B', and B'', we can verify that
the two results are the same.
"""
def test_B_Bprime_Bprimeprime(self):
rtol = 1e-8
atol = 1e-8

np.random.seed(0)
for sG in [-1, 1]:
for spsi in [-1, 1]:
for config in [1, 2, 3, 4, 5]:
print(sG, spsi, config)
B0 = np.random.rand() * 0.4 + 0.8
nphi = int(np.random.rand() * 50 + 61)
s = Qsc.from_paper(config, sG=sG, spsi=spsi, B0=B0, nphi=nphi)
s.calculate()

# this test checks that the derivatives B'(\phi) and B''(\phi)
# when calculated by from the gradB and gradgradB tensors are correct
axis_B0 = s.Bfield_cartesian(r=0, theta=0.)
axis_B1 = s.grad_B_tensor_cartesian()
axis_B2 = s.grad_grad_B_tensor_cartesian()

from simsopt.geo import CurveRZFourier
nfp=s.nfp
G0 = CurveRZFourier(np.linspace(0, 1/nfp, s.phi.size, endpoint=False), s.rc.size-1, nfp, True)
G0.rc[:] = s.rc[:]
G0.zs[:] = s.zs[1:]
G0.x = G0.get_dofs()

dG0_dphi = G0.gammadash()
d2G0_dphi2 = G0.gammadashdash()

ds_dphi = np.linalg.norm(dG0_dphi, axis=1)
d2s_dphi2 = np.sum(dG0_dphi*d2G0_dphi2, axis=1)/ds_dphi

T, N, B = G0.frenet_frame()
kappa = G0.kappa()
dkappa_dphi = G0.kappadash()
dkappa_ds = dkappa_dphi * (1/ds_dphi)
torsion = G0.torsion()

dT_ds = kappa[:, None] * N
dN_ds = -kappa[:, None] * T + torsion[:, None]*B
d2T_ds2 = dkappa_ds[:, None] * N + dN_ds*kappa[:, None]

dT_dphi = dT_ds * ds_dphi[:, None]
d2T_dphi2 = d2T_ds2*ds_dphi[:, None]**2 + d2s_dphi2[:, None]*dT_ds

Bexact = np.transpose(axis_B0, axes=(1, 0))
Bqsc = sG*B0*T

#check B(\phi)
np.testing.assert_allclose(Bexact, Bqsc, rtol=rtol, atol=atol)

dG0_dphi = np.transpose(dG0_dphi, axes=(1, 0))
dBqsc_dphi = np.einsum('ijk,ik->jk', axis_B1, dG0_dphi)
dBexact_dphi = sG*B0*dT_dphi
dBexact_dphi = np.transpose(dBexact_dphi, axes=(1, 0))
#check B'(\phi)
np.testing.assert_allclose(dBexact_dphi, dBqsc_dphi, rtol=rtol, atol=atol)

d2G0_dphi2 = np.transpose(d2G0_dphi2, axes=(1, 0))
d2Bqsc_dphi2 = np.einsum('ijkl,il,jl->kl', axis_B2, dG0_dphi, dG0_dphi)+np.einsum('ijl,jl->il', axis_B1, d2G0_dphi2)
d2Bexact_dphi2 = sG*B0*d2T_dphi2
d2Bexact_dphi2 = np.transpose(d2Bexact_dphi2, axes=(1, 0))
#check B''(\phi)
np.testing.assert_allclose(d2Bexact_dphi2, d2Bqsc_dphi2, rtol=rtol, atol=atol)

class CylindricalCartesianTensorsTests(unittest.TestCase):

Expand Down