Skip to content
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

Jump to bottom

Remove some repetition from Greens #3272

Merged
merged 6 commits into from
Aug 15, 2024
Merged

Remove some repetition from Greens #3272

Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
72238f9
remove some repitition
Apr 18, 2024
efbc9ae
ruff reformat
May 13, 2024
3a95919
remove accidental bracket addtion
May 13, 2024
56dff60
Update bluemira/magnetostatics/greens.py
geograham Aug 15, 2024
41fb36b
Update bluemira/magnetostatics/greens.py
geograham Aug 15, 2024
bd8b76f
Rearrange the 4s
Aug 15, 2024
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
76 changes: 48 additions & 28 deletions bluemira/magnetostatics/greens.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -196,15 +196,9 @@ def greens_psi(
\t:math:`\\mathbf{K} \\equiv` complete elliptic integral of the first kind\n
\t:math:`\\mathbf{E} \\equiv` complete elliptic integral of the second kind
"""
k2 = 4 * x * xc / ((x + xc) ** 2 + (z - zc) ** 2)
# Avoid NaN when coil on grid point
k2 = clip_nb(k2, GREENS_ZERO, 1.0 - GREENS_ZERO)
return (
MU_0_2PI
* np.sqrt(x * xc)
* ((2 - k2) * ellipk_nb(k2) - 2 * ellipe_nb(k2))
/ np.sqrt(k2)
)
_, k2 = calc_a_k2(xc, zc, x, z)
e, k = calc_e_k(k2)
return MU_0_2PI * np.sqrt(x * xc) * ((2 - k2) * k - 2 * e) / np.sqrt(k2)


@nb.jit(nopython=True)
Expand Down Expand Up @@ -256,12 +250,8 @@ def greens_dpsi_dx(

The implementation used here refactors the above to avoid some zero divisions.
"""
a = ((x + xc) ** 2 + (z - zc) ** 2) ** 0.5
k2 = 4 * x * xc / a**2
# Avoid NaN when coil on grid point
k2 = clip_nb(k2, GREENS_ZERO, 1.0 - GREENS_ZERO)
i1 = ellipk_nb(k2) / a
i2 = ellipe_nb(k2) / (a**3 * (1 - k2))
a, k2 = calc_a_k2(xc, zc, x, z)
i1, i2 = calc_i1_i2(a, k2)
return MU_0_2PI * x * ((xc**2 - (z - zc) ** 2 - x**2) * i2 + i1)


Expand Down Expand Up @@ -314,12 +304,44 @@ def greens_dpsi_dz(

The implementation used here refactors the above to avoid some zero divisions.
"""
a = ((x + xc) ** 2 + (z - zc) ** 2) ** 0.5
a, k2 = calc_a_k2(xc, zc, x, z)
i1, i2 = calc_i1_i2(a, k2)
return MU_0_2PI * ((z - zc) * (i1 - i2 * ((z - zc) ** 2 + x**2 + xc**2)))


@nb.jit(nopython=True)
def calc_a_k2(
xc: float | np.ndarray,
zc: float | np.ndarray,
x: float | np.ndarray,
z: float | np.ndarray,
):
a = np.hypot((x + xc), (z - zc))
k2 = 4 * x * xc / a**2
# Avoid NaN when coil on grid point
k2 = clip_nb(k2, GREENS_ZERO, 1.0 - GREENS_ZERO)
i1 = ellipk_nb(k2) / a
i2 = ellipe_nb(k2) / (a**3 * (1 - k2))
return MU_0_2PI * ((z - zc) * (i1 - i2 * ((z - zc) ** 2 + x**2 + xc**2)))
return a, k2


@nb.jit(nopython=True)
def calc_e_k(
k2: float | np.ndarray,
):
return ellipe_nb(k2), ellipk_nb(k2)
geograham marked this conversation as resolved.
Show resolved Hide resolved


@nb.jit(nopython=True)
def calc_i1_i2(
a: float | np.ndarray,
k2: float | np.ndarray,
e: float | np.ndarray | None = None,
k: float | np.ndarray | None = None,
):
if (e is None) or (k is None):
e, k = calc_e_k(k2)
i1 = k / a
i2 = e / (a**3 * (1 - k2))
geograham marked this conversation as resolved.
Show resolved Hide resolved
return i1, i2


@nb.jit(nopython=True)
Expand Down Expand Up @@ -450,16 +472,14 @@ def greens_all(
if xc <= 0
if x <= 0
"""
a = np.hypot((x + xc), (z - zc))
k2 = 4 * x * xc / a**2
# Avoid NaN when coil on grid point
k2 = clip_nb(k2, GREENS_ZERO, 1.0 - GREENS_ZERO)
e, k = ellipe_nb(k2), ellipk_nb(k2)
i_1 = 4 * k / a
i_2 = 4 * e / (a**3 * (1 - k2))
a, k2 = calc_a_k2(xc, zc, x, z)
e, k = calc_e_k(k2)
i1, i2 = calc_i1_i2(a, k2, e, k)
geograham marked this conversation as resolved.
Show resolved Hide resolved
i1 *= 4
i2 *= 4
a_part = (z - zc) ** 2 + x**2 + xc**2
b_part = -2 * x * xc
g_bx = MU_0_4PI * xc * (z - zc) * (i_1 - i_2 * a_part) / b_part
g_bz = MU_0_4PI * xc * ((xc + x * a_part / b_part) * i_2 - i_1 * x / b_part)
g_bx = MU_0_4PI * xc * (z - zc) * (i1 - i2 * a_part) / b_part
g_bz = MU_0_4PI * xc * ((xc + x * a_part / b_part) * i2 - i1 * x / b_part)
g_psi = MU_0_4PI * a * ((2 - k2) * k - 2 * e)
return g_psi, g_bx, g_bz
Loading