From 1c79c640b7d8aaf19c916c6b4e99cc10ab3692a7 Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:18:42 -0500 Subject: [PATCH 01/18] Synch improvements for helicity, q(r) discritization, and d_phi advance step --- .../magnetic_geometry/equilibrium_field.py | 54 +++++++++---------- synthwave/magnetic_geometry/filaments.py | 42 +++++++++++---- 2 files changed, 60 insertions(+), 36 deletions(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index e9e5924..5bec7c0 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -42,20 +42,29 @@ def biot_savart_cylindrical( class EquilibriumField: - def __init__(self, eqdsk): + def __init__(self, eqdsk,lam=1e-7): self.eqdsk = eqdsk self.psi = RectBivariateSpline( - eqdsk.r_grid[:, 0], eqdsk.z_grid[0, :], eqdsk.psi, kx=3, ky=3, s=0 + eqdsk.r_grid[:, 0], eqdsk.z_grid[0, :], eqdsk.psi, + kx=3, ky=3, s=0 ) # Linear grid of psi for 1D profiles # https://freeqdsk.readthedocs.io/en/stable/geqdsk.html self.psi_grid = np.linspace(eqdsk.simagx, eqdsk.sibdry, eqdsk.nx) - + # q(psi) - self.qpsi = make_interp_spline(self.psi_grid, eqdsk.qpsi, k=3, axis=0) + # RNC EDIT: Switching to smoothing spline to avoid strange "discretization" jumps + # Note: lam value (lower = less smoothing) should be reasonably consistent across q + # profiles, but this is not certain. Lam=1e-7 works for q(psi) and F(psi) so far. + self.qpsi = make_smoothing_spline( + self.psi_grid, eqdsk.qpsi, lam=lam, axis=0 + ) + # F(psi) - self.F = make_interp_spline(self.psi_grid, eqdsk.fpol, k=3, axis=0) + self.F = make_smoothing_spline( + self.psi_grid, eqdsk.fpol, lam=lam, axis=0 + ) def get_field_at_point(self, R, Z) -> np.ndarray: # Bp = Br + Bz = (d(psi)/dZ - d(psi)/dR) / R @@ -68,7 +77,7 @@ def get_field_at_point(self, R, Z) -> np.ndarray: return np.array([Br, Bt, Bz]) - def get_psi_of_q(self, q, tol=1e-3): + def get_psi_of_q(self, q): """Get psi corresponding to a given q""" qpsi_grid = self.qpsi(self.psi_grid) psi_guess = self.psi_grid[np.argmin(np.abs(qpsi_grid - q))] @@ -76,13 +85,12 @@ def get_psi_of_q(self, q, tol=1e-3): func=lambda psi: self.qpsi(psi) - q, x0=psi_guess, fprime=lambda psi: self.qpsi.derivative(1)(psi), - maxiter=400, - tol=tol, + maxiter = 400,tol=1e-3 ) - # RNC EDIT: - # CHECK: For q~<=1, depening on the resolution of the gEQDSK file, - # the interpolation function can request a psi value at or less + # RNC EDIT: + # CHECK: For q~<=1, depening on the resolution of the gEQDSK file, + # the interpolation function can request a psi value at or less # than the minimum in the psi_grid vector # Check to see if the value we want plausibly exists in the final set # of q values from the eqdsk @@ -90,19 +98,11 @@ def get_psi_of_q(self, q, tol=1e-3): # q values are "grouped" in an odd, stepwise fashion if psi <= self.psi_grid[0]: # check if there's a range of possible q-values (e.g.) if this is plausably a resolution issue - if ( - np.argwhere(qpsi_grid > (qpsi_grid[0] + tol)).squeeze()[0] > 1 - ): # multiple identical q values in a row - lin_interp_q = np.polyfit(self.psi_grid[:30], qpsi_grid[:30], 1) - psi = self.psi_grid[ - np.argmin(np.abs(np.polyval(lin_interp_q, self.psi_grid[:30]) - q)) - ] - - if psi > self.psi_grid[0]: - return psi - else: - raise ValueError( - "Error: requested q=%1.3f is outside the gEQDSK range (q_min = %1.3f)" - % (q, qpsi_grid[0]) - ) - return psi + tol = 1e-3 + if np.argwhere(qpsi_grid>(qpsi_grid[0]+1e-3)).squeeze()[0] > 1: # multiple identical q values in a row + lin_interp_q = np.polyfit(self.psi_grid[:30],qpsi_grid[:30], 1) + psi = self.psi_grid[ np.argmin(np.abs(np.polyval(lin_interp_q,self.psi_grid[:30]) - q)) ] + + if psi > self.psi_grid[0]: return psi + else: raise SyntaxError('Error: requested q=%1.3f is outside the gEQDSK range (q_min = %1.3f)'%(q,qpsi_grid[0])) + return psi \ No newline at end of file diff --git a/synthwave/magnetic_geometry/filaments.py b/synthwave/magnetic_geometry/filaments.py index 5cc78d4..8ff8c23 100644 --- a/synthwave/magnetic_geometry/filaments.py +++ b/synthwave/magnetic_geometry/filaments.py @@ -284,6 +284,7 @@ def _trace( self, num_filament_points: Optional[int] = None, method: TraceType = TraceType.SINGLE, + helicity_sign: int = +1 ) -> tuple[np.ndarray, np.ndarray]: """Trace a filament""" if num_filament_points is None: @@ -322,17 +323,17 @@ def _R_a(eta, a): R = self.eq_field.eqdsk.rmagx + (a * np.cos(eta)) return R + # Currently: +m/+n helicity matches empirical C-Mod pickup def _Z_a(eta, a): - Z = self.eq_field.eqdsk.zmagx - (a * np.sin(eta)) - return Z + Z = self.eq_field.eqdsk.zmagx + (helicity_sign) * (a * np.sin(eta)) + return Z def psi_prime_a(eta, a): # Derivative of psi with respect to a at a given eta R = _R_a(eta, a) Z = _Z_a(eta, a) - return self.eq_field.psi.ev(R, Z, dx=1, dy=0) * np.cos( - eta - ) - self.eq_field.psi.ev(R, Z, dx=0, dy=1) * np.sin(eta) + return self.eq_field.psi.ev(R, Z, dx=1, dy=0) * np.cos(eta) + \ + (helicity_sign) * self.eq_field.psi.ev(R, Z, dx=0, dy=1) * np.sin(eta) for i, eta in enumerate(filament_etas): if i == 0: @@ -349,13 +350,19 @@ def psi_prime_a(eta, a): func=lambda a: self.eq_field.psi.ev(_R_a(eta, a), _Z_a(eta, a)) - psi_q, x0=a_guess, fprime=lambda a: psi_prime_a(eta, a), + maxiter=800, + tol=1e-3, ) poloidal_points[i, :] = [_R_a(eta, a_next), _Z_a(eta, a_next), a_next] def _d_phi(r, R, Bp, Bt, d_eta): # https://wiki.fusion.ciemat.es/wiki/Rotational_transform return (Bt * r * d_eta) / (R * Bp) - + #alternative form removing the assumption that dl = r d_eta (that assumption holds only for circular cross-sections) + def _d_phi_dl(dl, R, Bp, Bt): + # https://youjunhu.github.io/research_notes/tokamak_equilibrium_htlatex/tokamak_equilibrium.html + return (Bt * dl) / (R * Bp) + # Finalize filament trace based on method if method == EquilibriumFilamentTracer.TraceType.CYLINDRICAL: # Circular cross section around the magnetic axis @@ -379,8 +386,19 @@ def _d_phi(r, R, Bp, Bt, d_eta): Z = poloidal_points[:, 1] r = poloidal_points[:, 2] B = self.eq_field.get_field_at_point(R, Z) - d_eta = np.mean(np.diff(filament_etas)) - d_phi = _d_phi(r, R, np.sqrt(B[0] ** 2 + B[2] ** 2), B[1], d_eta) + + + # Switching to improved d_phi formula from the below: + # d_eta = np.mean(np.diff(filament_etas)) + # d_phi = _d_phi(r, R, np.sqrt(B[0] ** 2 + B[2] ** 2), B[1], d_eta) + + # Compute segment lengths with wraparound so the last segment goes from + # the final point back to the first + dR = np.roll(R, -1) - R + dZ = np.roll(Z, -1) - Z + dl = np.sqrt(dR**2 + dZ**2) + d_phi = _d_phi_dl(dl, R, np.sqrt(B[0]**2 + B[2]**2), B[1]) + if method == EquilibriumFilamentTracer.TraceType.SINGLE: phi = np.cumsum(d_phi) - d_phi[0] else: @@ -391,12 +409,17 @@ def _d_phi(r, R, Bp, Bt, d_eta): # Numerical correction to ensure final point is at the proper angle # We can do this multiple times, whenever we know for sure that the phi is a multiple of pi * m / n - # RNC EDIT: Sign flip necessary to account for helicity direction + # Sign flip necessary to account for helicity direction + # Note: The improved d_phi method should ensure this is always correct now [e.g. the correction factor is no longer strictly necessary] + known_phis = np.linspace( 0, 2 * np.pi * self.m_local, (2 * self.n_local) + 1 ) * np.sign(phi[-1]) + for i, known_phi_start in enumerate(known_phis[:-1]): + known_phi_end = known_phis[i + 1] + if np.sign(phi[-1]) == 1: phi_indices = np.squeeze( np.where((phi >= known_phi_start) & (phi <= known_phi_end)) @@ -411,6 +434,7 @@ def _d_phi(r, R, Bp, Bt, d_eta): correction_factor = (known_phi_end - known_phi_start) / ( actual_phi_end - actual_phi_start ) + phi[phi_indices] = ( known_phi_start + (phi[phi_indices] - actual_phi_start) * correction_factor From e6b0e1669b8ca10923ab0df022bd133d6d0ae110 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 26 Jan 2026 21:28:16 +0000 Subject: [PATCH 02/18] Initial plan From 9827a47944959023e161d04e120370bd7fdf2506 Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:28:26 -0500 Subject: [PATCH 03/18] Update synthwave/magnetic_geometry/equilibrium_field.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/equilibrium_field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index 5bec7c0..4655a89 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -42,7 +42,7 @@ def biot_savart_cylindrical( class EquilibriumField: - def __init__(self, eqdsk,lam=1e-7): + def __init__(self, eqdsk, lam=1e-7): self.eqdsk = eqdsk self.psi = RectBivariateSpline( eqdsk.r_grid[:, 0], eqdsk.z_grid[0, :], eqdsk.psi, From cd4e5d827642cbb03efe7f2388ce1a47dd43968b Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:28:44 -0500 Subject: [PATCH 04/18] Update synthwave/magnetic_geometry/equilibrium_field.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/equilibrium_field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index 4655a89..0c67359 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -88,7 +88,7 @@ def get_psi_of_q(self, q): maxiter = 400,tol=1e-3 ) - # RNC EDIT: + # RNC EDIT: # CHECK: For q~<=1, depening on the resolution of the gEQDSK file, # the interpolation function can request a psi value at or less # than the minimum in the psi_grid vector From 1359fd5dfbb43121c7ce42a9045ec71f303223f5 Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:28:51 -0500 Subject: [PATCH 05/18] Update synthwave/magnetic_geometry/equilibrium_field.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/equilibrium_field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index 0c67359..2c0642b 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -85,7 +85,7 @@ def get_psi_of_q(self, q): func=lambda psi: self.qpsi(psi) - q, x0=psi_guess, fprime=lambda psi: self.qpsi.derivative(1)(psi), - maxiter = 400,tol=1e-3 + maxiter=400, tol=1e-3 ) # RNC EDIT: From 749b75cf209ae755906c7d534cc08c1e9f77f0a1 Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:29:01 -0500 Subject: [PATCH 06/18] Update synthwave/magnetic_geometry/filaments.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/filaments.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/synthwave/magnetic_geometry/filaments.py b/synthwave/magnetic_geometry/filaments.py index 8ff8c23..65371ee 100644 --- a/synthwave/magnetic_geometry/filaments.py +++ b/synthwave/magnetic_geometry/filaments.py @@ -326,7 +326,7 @@ def _R_a(eta, a): # Currently: +m/+n helicity matches empirical C-Mod pickup def _Z_a(eta, a): Z = self.eq_field.eqdsk.zmagx + (helicity_sign) * (a * np.sin(eta)) - return Z + return Z def psi_prime_a(eta, a): # Derivative of psi with respect to a at a given eta From 9dbea60922c78889bf8ddbe80a122b96e395c27f Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:29:39 -0500 Subject: [PATCH 07/18] Update synthwave/magnetic_geometry/equilibrium_field.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/equilibrium_field.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index 2c0642b..ab279ff 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -103,6 +103,11 @@ def get_psi_of_q(self, q): lin_interp_q = np.polyfit(self.psi_grid[:30],qpsi_grid[:30], 1) psi = self.psi_grid[ np.argmin(np.abs(np.polyval(lin_interp_q,self.psi_grid[:30]) - q)) ] - if psi > self.psi_grid[0]: return psi - else: raise SyntaxError('Error: requested q=%1.3f is outside the gEQDSK range (q_min = %1.3f)'%(q,qpsi_grid[0])) + if psi > self.psi_grid[0]: + return psi + else: + raise SyntaxError( + 'Error: requested q=%1.3f is outside the gEQDSK range (q_min = %1.3f)' + % (q, qpsi_grid[0]) + ) return psi \ No newline at end of file From 1844e51ed6db5fcfa5672cb2cb2f0fb15b829eae Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:29:49 -0500 Subject: [PATCH 08/18] Update synthwave/magnetic_geometry/equilibrium_field.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/equilibrium_field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index ab279ff..0e89173 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -89,7 +89,7 @@ def get_psi_of_q(self, q): ) # RNC EDIT: - # CHECK: For q~<=1, depening on the resolution of the gEQDSK file, + # CHECK: For q~<=1, depending on the resolution of the gEQDSK file, # the interpolation function can request a psi value at or less # than the minimum in the psi_grid vector # Check to see if the value we want plausibly exists in the final set From 0b594113b1755b6a5b9f7785951b42c71c866359 Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:29:59 -0500 Subject: [PATCH 09/18] Update synthwave/magnetic_geometry/filaments.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/filaments.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/synthwave/magnetic_geometry/filaments.py b/synthwave/magnetic_geometry/filaments.py index 65371ee..3ca42b4 100644 --- a/synthwave/magnetic_geometry/filaments.py +++ b/synthwave/magnetic_geometry/filaments.py @@ -284,7 +284,7 @@ def _trace( self, num_filament_points: Optional[int] = None, method: TraceType = TraceType.SINGLE, - helicity_sign: int = +1 + helicity_sign: int = +1, ) -> tuple[np.ndarray, np.ndarray]: """Trace a filament""" if num_filament_points is None: From cd942aa2ba3210cf014e02b841aa6820351c5cc2 Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:30:11 -0500 Subject: [PATCH 10/18] Update synthwave/magnetic_geometry/equilibrium_field.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/equilibrium_field.py | 1 - 1 file changed, 1 deletion(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index 0e89173..d05b799 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -52,7 +52,6 @@ def __init__(self, eqdsk, lam=1e-7): # Linear grid of psi for 1D profiles # https://freeqdsk.readthedocs.io/en/stable/geqdsk.html self.psi_grid = np.linspace(eqdsk.simagx, eqdsk.sibdry, eqdsk.nx) - # q(psi) # RNC EDIT: Switching to smoothing spline to avoid strange "discretization" jumps # Note: lam value (lower = less smoothing) should be reasonably consistent across q From 088fd1abab81f8ec991c4f1e15c4b94c99b0bb35 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 26 Jan 2026 21:30:47 +0000 Subject: [PATCH 11/18] Add missing import for make_smoothing_spline Co-authored-by: chandrarn <8441592+chandrarn@users.noreply.github.com> --- synthwave/magnetic_geometry/equilibrium_field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index 5bec7c0..3df5f14 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -1,7 +1,7 @@ import numpy as np from scipy.constants import mu_0 -from scipy.interpolate import RectBivariateSpline, make_interp_spline +from scipy.interpolate import RectBivariateSpline, make_interp_spline, make_smoothing_spline from scipy.optimize import newton from synthwave.magnetic_geometry.utils import ( From 8d0f4a85891ee8a3c7f0d1a9b09a30eb13d9500f Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 26 Jan 2026 21:38:48 +0000 Subject: [PATCH 12/18] Initial plan From 37766ecaaebac527802cdc65428f53810b462d21 Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:45:29 -0500 Subject: [PATCH 13/18] Update synthwave/magnetic_geometry/equilibrium_field.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/equilibrium_field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index 9ff2c4f..8483522 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -1,7 +1,7 @@ import numpy as np from scipy.constants import mu_0 -from scipy.interpolate import RectBivariateSpline, make_interp_spline, make_smoothing_spline +from scipy.interpolate import RectBivariateSpline, make_smoothing_spline from scipy.optimize import newton from synthwave.magnetic_geometry.utils import ( From 1b50b4797e4c89669e2b252be6b29b67b156d0bd Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:45:54 -0500 Subject: [PATCH 14/18] Update synthwave/magnetic_geometry/filaments.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/filaments.py | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/synthwave/magnetic_geometry/filaments.py b/synthwave/magnetic_geometry/filaments.py index 3ca42b4..263c9fa 100644 --- a/synthwave/magnetic_geometry/filaments.py +++ b/synthwave/magnetic_geometry/filaments.py @@ -286,7 +286,28 @@ def _trace( method: TraceType = TraceType.SINGLE, helicity_sign: int = +1, ) -> tuple[np.ndarray, np.ndarray]: - """Trace a filament""" + """Trace a filament along the equilibrium magnetic field. + + Parameters + ---------- + num_filament_points : int, optional + Number of points to use for tracing a single poloidal turn of the filament. + If ``None``, the value stored in ``self.num_points`` is used. + method : EquilibriumFilamentTracer.TraceType, optional + Tracing strategy to use. This controls how the field is followed when + computing the filament shape (e.g., cylindrical approximation, naive, + single-point, or averaged tracing). + helicity_sign : int, optional + Sign of the helicity used when following the field lines. A value of + ``+1`` traces in the default direction, while ``-1`` reverses the + direction. + + Returns + ------- + tuple of numpy.ndarray + Tuple containing arrays describing the traced filament coordinates. + + """ if num_filament_points is None: num_filament_points = self.num_points From 2065258617b01d32980eee271fc701a28befcac4 Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:46:05 -0500 Subject: [PATCH 15/18] Update synthwave/magnetic_geometry/equilibrium_field.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/equilibrium_field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index 8483522..2214e79 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -99,7 +99,7 @@ def get_psi_of_q(self, q): # check if there's a range of possible q-values (e.g.) if this is plausably a resolution issue tol = 1e-3 if np.argwhere(qpsi_grid>(qpsi_grid[0]+1e-3)).squeeze()[0] > 1: # multiple identical q values in a row - lin_interp_q = np.polyfit(self.psi_grid[:30],qpsi_grid[:30], 1) + lin_interp_q = np.polyfit(self.psi_grid[:30], qpsi_grid[:30], 1) psi = self.psi_grid[ np.argmin(np.abs(np.polyval(lin_interp_q,self.psi_grid[:30]) - q)) ] if psi > self.psi_grid[0]: From ece2b8ed1bab61844738bd14890d1ec6a2bf1556 Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:46:16 -0500 Subject: [PATCH 16/18] Update synthwave/magnetic_geometry/equilibrium_field.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/equilibrium_field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index 2214e79..01a0546 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -96,7 +96,7 @@ def get_psi_of_q(self, q): # The issue appears to be that although psi is continuous and monotonic, # q values are "grouped" in an odd, stepwise fashion if psi <= self.psi_grid[0]: - # check if there's a range of possible q-values (e.g.) if this is plausably a resolution issue + # check if there's a range of possible q-values (e.g.) if this is plausibly a resolution issue tol = 1e-3 if np.argwhere(qpsi_grid>(qpsi_grid[0]+1e-3)).squeeze()[0] > 1: # multiple identical q values in a row lin_interp_q = np.polyfit(self.psi_grid[:30], qpsi_grid[:30], 1) From d5ff5fdef372acacad97b9465581ae453f12c3be Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 16:46:39 -0500 Subject: [PATCH 17/18] Update synthwave/magnetic_geometry/equilibrium_field.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- synthwave/magnetic_geometry/equilibrium_field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index 01a0546..4f7ac10 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -105,7 +105,7 @@ def get_psi_of_q(self, q): if psi > self.psi_grid[0]: return psi else: - raise SyntaxError( + raise ValueError( 'Error: requested q=%1.3f is outside the gEQDSK range (q_min = %1.3f)' % (q, qpsi_grid[0]) ) From a12345dddbb2437fd7f9b84f3dddc61bc832461a Mon Sep 17 00:00:00 2001 From: Rian Chandra Date: Mon, 26 Jan 2026 19:23:39 -0500 Subject: [PATCH 18/18] Synch post pre-commit check --- .../magnetic_geometry/equilibrium_field.py | 36 ++++----- synthwave/magnetic_geometry/filaments.py | 15 ++-- synthwave/mirnov/run_thincurr_model.py | 77 +++++++++++-------- 3 files changed, 72 insertions(+), 56 deletions(-) diff --git a/synthwave/magnetic_geometry/equilibrium_field.py b/synthwave/magnetic_geometry/equilibrium_field.py index 4f7ac10..82c362f 100644 --- a/synthwave/magnetic_geometry/equilibrium_field.py +++ b/synthwave/magnetic_geometry/equilibrium_field.py @@ -45,8 +45,7 @@ class EquilibriumField: def __init__(self, eqdsk, lam=1e-7): self.eqdsk = eqdsk self.psi = RectBivariateSpline( - eqdsk.r_grid[:, 0], eqdsk.z_grid[0, :], eqdsk.psi, - kx=3, ky=3, s=0 + eqdsk.r_grid[:, 0], eqdsk.z_grid[0, :], eqdsk.psi, kx=3, ky=3, s=0 ) # Linear grid of psi for 1D profiles @@ -54,16 +53,12 @@ def __init__(self, eqdsk, lam=1e-7): self.psi_grid = np.linspace(eqdsk.simagx, eqdsk.sibdry, eqdsk.nx) # q(psi) # RNC EDIT: Switching to smoothing spline to avoid strange "discretization" jumps - # Note: lam value (lower = less smoothing) should be reasonably consistent across q + # Note: lam value (lower = less smoothing) should be reasonably consistent across q # profiles, but this is not certain. Lam=1e-7 works for q(psi) and F(psi) so far. - self.qpsi = make_smoothing_spline( - self.psi_grid, eqdsk.qpsi, lam=lam, axis=0 - ) + self.qpsi = make_smoothing_spline(self.psi_grid, eqdsk.qpsi, lam=lam, axis=0) # F(psi) - self.F = make_smoothing_spline( - self.psi_grid, eqdsk.fpol, lam=lam, axis=0 - ) + self.F = make_smoothing_spline(self.psi_grid, eqdsk.fpol, lam=lam, axis=0) def get_field_at_point(self, R, Z) -> np.ndarray: # Bp = Br + Bz = (d(psi)/dZ - d(psi)/dR) / R @@ -84,12 +79,13 @@ def get_psi_of_q(self, q): func=lambda psi: self.qpsi(psi) - q, x0=psi_guess, fprime=lambda psi: self.qpsi.derivative(1)(psi), - maxiter=400, tol=1e-3 + maxiter=400, + tol=1e-3, ) # RNC EDIT: - # CHECK: For q~<=1, depending on the resolution of the gEQDSK file, - # the interpolation function can request a psi value at or less + # CHECK: For q~<=1, depending on the resolution of the gEQDSK file, + # the interpolation function can request a psi value at or less # than the minimum in the psi_grid vector # Check to see if the value we want plausibly exists in the final set # of q values from the eqdsk @@ -97,16 +93,20 @@ def get_psi_of_q(self, q): # q values are "grouped" in an odd, stepwise fashion if psi <= self.psi_grid[0]: # check if there's a range of possible q-values (e.g.) if this is plausibly a resolution issue - tol = 1e-3 - if np.argwhere(qpsi_grid>(qpsi_grid[0]+1e-3)).squeeze()[0] > 1: # multiple identical q values in a row + + if ( + np.argwhere(qpsi_grid > (qpsi_grid[0] + 1e-3)).squeeze()[0] > 1 + ): # multiple identical q values in a row lin_interp_q = np.polyfit(self.psi_grid[:30], qpsi_grid[:30], 1) - psi = self.psi_grid[ np.argmin(np.abs(np.polyval(lin_interp_q,self.psi_grid[:30]) - q)) ] - + psi = self.psi_grid[ + np.argmin(np.abs(np.polyval(lin_interp_q, self.psi_grid[:30]) - q)) + ] + if psi > self.psi_grid[0]: return psi else: raise ValueError( - 'Error: requested q=%1.3f is outside the gEQDSK range (q_min = %1.3f)' + "Error: requested q=%1.3f is outside the gEQDSK range (q_min = %1.3f)" % (q, qpsi_grid[0]) ) - return psi \ No newline at end of file + return psi diff --git a/synthwave/magnetic_geometry/filaments.py b/synthwave/magnetic_geometry/filaments.py index 263c9fa..6c1f573 100644 --- a/synthwave/magnetic_geometry/filaments.py +++ b/synthwave/magnetic_geometry/filaments.py @@ -353,8 +353,9 @@ def psi_prime_a(eta, a): # Derivative of psi with respect to a at a given eta R = _R_a(eta, a) Z = _Z_a(eta, a) - return self.eq_field.psi.ev(R, Z, dx=1, dy=0) * np.cos(eta) + \ - (helicity_sign) * self.eq_field.psi.ev(R, Z, dx=0, dy=1) * np.sin(eta) + return self.eq_field.psi.ev(R, Z, dx=1, dy=0) * np.cos(eta) + ( + helicity_sign + ) * self.eq_field.psi.ev(R, Z, dx=0, dy=1) * np.sin(eta) for i, eta in enumerate(filament_etas): if i == 0: @@ -379,11 +380,12 @@ def psi_prime_a(eta, a): def _d_phi(r, R, Bp, Bt, d_eta): # https://wiki.fusion.ciemat.es/wiki/Rotational_transform return (Bt * r * d_eta) / (R * Bp) - #alternative form removing the assumption that dl = r d_eta (that assumption holds only for circular cross-sections) + + # alternative form removing the assumption that dl = r d_eta (that assumption holds only for circular cross-sections) def _d_phi_dl(dl, R, Bp, Bt): # https://youjunhu.github.io/research_notes/tokamak_equilibrium_htlatex/tokamak_equilibrium.html return (Bt * dl) / (R * Bp) - + # Finalize filament trace based on method if method == EquilibriumFilamentTracer.TraceType.CYLINDRICAL: # Circular cross section around the magnetic axis @@ -405,9 +407,7 @@ def _d_phi_dl(dl, R, Bp, Bt): # determine d(phi)/d(eta) from magnetic field R = poloidal_points[:, 0] Z = poloidal_points[:, 1] - r = poloidal_points[:, 2] B = self.eq_field.get_field_at_point(R, Z) - # Switching to improved d_phi formula from the below: # d_eta = np.mean(np.diff(filament_etas)) @@ -418,7 +418,7 @@ def _d_phi_dl(dl, R, Bp, Bt): dR = np.roll(R, -1) - R dZ = np.roll(Z, -1) - Z dl = np.sqrt(dR**2 + dZ**2) - d_phi = _d_phi_dl(dl, R, np.sqrt(B[0]**2 + B[2]**2), B[1]) + d_phi = _d_phi_dl(dl, R, np.sqrt(B[0] ** 2 + B[2] ** 2), B[1]) if method == EquilibriumFilamentTracer.TraceType.SINGLE: phi = np.cumsum(d_phi) - d_phi[0] @@ -438,7 +438,6 @@ def _d_phi_dl(dl, R, Bp, Bt): ) * np.sign(phi[-1]) for i, known_phi_start in enumerate(known_phis[:-1]): - known_phi_end = known_phis[i + 1] if np.sign(phi[-1]) == 1: diff --git a/synthwave/mirnov/run_thincurr_model.py b/synthwave/mirnov/run_thincurr_model.py index 3834bef..e406d3e 100644 --- a/synthwave/mirnov/run_thincurr_model.py +++ b/synthwave/mirnov/run_thincurr_model.py @@ -377,57 +377,74 @@ def correct_frequency_response( ) ) + ################################################################################################ ################################################################################################ def plot_sensor_output( - working_files_directory, - probe_details, - mode, - freq, - save_Ext, - doSave, + working_files_directory, + probe_details, + mode, + freq, + save_Ext, + doSave, ): - # Load in saved sensor signals from netCDF format and plot - fName = working_files_directory + \ - "probe_signals_%s_m%02d_n%02d_f%1.1ekHz%s.nc" % \ - ( - probe_details.attrs["probe_set_name"], - mode["m"], - mode["n"], - freq / 1e3, - save_Ext, - ) + fName = working_files_directory + "probe_signals_%s_m%02d_n%02d_f%1.1ekHz%s.nc" % ( + probe_details.attrs["probe_set_name"], + mode["m"], + mode["n"], + freq / 1e3, + save_Ext, + ) # Load data sensors_bode = xr.load_dataset(fName) # Generate signal magnitudes - mags = [ np.linalg.norm(sensors_bode.sel(sensor=sig).signal.values.tolist()) \ - for sig in sensors_bode.sensor.values] + mags = [ + np.linalg.norm(sensors_bode.sel(sensor=sig).signal.values.tolist()) + for sig in sensors_bode.sensor.values + ] sensor_names = sensors_bode.sensor.values.tolist() # Plot - plt.close('Sensor_Signals_m%02d_n%02d_f%1.1ekHz'%( mode["m"],mode["n"],freq/1e3)) - fig,ax = plt.subplots(1,1,figsize=(4,3),tight_layout=True,num='Sensor_Signals_m%02d_n%02d_f%1.1ekHz'%( - mode["m"],mode["n"],freq/1e3)) - sc = ax.plot(sensor_names, mags, '-*',label=r'$m/n=%d/%d,\,f=%1.1f$ kHz'%(mode["m"],mode["n"],freq/1e3)) - + plt.close( + "Sensor_Signals_m%02d_n%02d_f%1.1ekHz" % (mode["m"], mode["n"], freq / 1e3) + ) + + fig, ax = plt.subplots( + 1, + 1, + figsize=(4, 3), + tight_layout=True, + num="Sensor_Signals_m%02d_n%02d_f%1.1ekHz" % (mode["m"], mode["n"], freq / 1e3), + ) + + ax.plot( + sensor_names, + mags, + "-*", + label=r"$m/n=%d/%d,\,f=%1.1f$ kHz" % (mode["m"], mode["n"], freq / 1e3), + ) + ax.set_xticks(range(0, len(sensor_names), 3)) ax.set_xticklabels([sensor_names[i] for i in range(0, len(sensor_names), 3)]) - ax.tick_params(axis='x', rotation=90) - ax.legend(fontsize=8,handlelength=1) - ax.set_ylabel('Signal Magnitude [T/s]') + ax.tick_params(axis="x", rotation=90) + ax.legend(fontsize=8, handlelength=1) + ax.set_ylabel("Signal Magnitude [T/s]") ax.grid() if doSave: - fig.savefig(working_files_directory + - "Sensor_Signals_m%02d_n%02d_f%1.1ekHz%s.svg" % - ( + fig.savefig( + working_files_directory + + "Sensor_Signals_m%02d_n%02d_f%1.1ekHz%s.svg" + % ( mode["m"], mode["n"], freq / 1e3, save_Ext, - ), transparent=True) + ), + transparent=True, + )