Skip to content

Fix grad grad b tensor calculation #13

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
rogeriojorge wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Fix grad grad b tensor calculation #13

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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 qic/grad_B_tensor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1910,85 +1910,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))
Copy link

Copilot AI Jun 30, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] If the commented-out code is no longer needed, it would be cleaner to remove it to reduce clutter and improve maintainability.

Suggested change
# return np.transpose(self.grad_grad_B,(1,2,3,0))

Copilot uses AI. Check for mistakes.
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)
Copy link

Copilot AI Jun 30, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The final transpose call in grad_grad_B_tensor_cartesian changes the tensor shape from (nphi, 3, 3, 3) to (3, 3, 3, nphi), which does not match the documented output shape. Please adjust the transpose to ensure the returned tensor has shape (nphi, 3, 3, 3) as specified in the docstring.

Suggested change
return T_cartesian.transpose(1,2,3,0) # shape (nphi, 3, 3, 3)
return T_cartesian # shape (nphi, 3, 3, 3)

Copilot uses AI. Check for mistakes.
Loading