From 5274d22ae0b51d07ce14683a25fbc244f4c2da16 Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Tue, 9 Jul 2024 11:38:15 -0400 Subject: [PATCH 01/19] add from_values to FourierRZCoil --- desc/coils.py | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/desc/coils.py b/desc/coils.py index 1629941789..7431e18058 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -401,6 +401,41 @@ def __init__( ): super().__init__(current, R_n, Z_n, modes_R, modes_Z, NFP, sym, name) + @classmethod + def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", name="", sym=False): + """Fit coordinates to FourierRZCoil representation. + + Parameters + ---------- + current : float + Current through the coil, in Amps. + coords: ndarray, shape (num_coords,3) + coordinates to fit a FourierRZCurve object with each column + corresponding to xyz or rpz depending on the basis argument. + N : int + Fourier resolution of the new R,Z representation. + NFP : int + Number of field periods, the curve will have a discrete toroidal symmetry + according to NFP. + basis : {"rpz", "xyz"} + basis for input coordinates. Defaults to "rpz" + + Returns + ------- + coili : FourierRZCoil + New representation of the coil parameterized by Fourier series for R,Z. + + """ + curve = super().from_values(coords, N=N, NFP=NFP, basis=basis, sym=sym) + return cls( + current, + R_n=curve.R_n, + Z_n=curve.Z_n, + modes_R=curve.R_basis.modes[:, 2], + modes_Z=curve.Z_basis.modes[:, 2], + name=name, + ) + class FourierXYZCoil(_Coil, FourierXYZCurve): """Coil parameterized by fourier series for X,Y,Z in terms of arbitrary angle s. From 4a870a452b44299afad67fc58fa8f1ab8c21c3e7 Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Tue, 9 Jul 2024 11:41:27 -0400 Subject: [PATCH 02/19] add missing attrs --- desc/coils.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/desc/coils.py b/desc/coils.py index 7431e18058..bdfb2363ac 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -433,6 +433,8 @@ def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", name="", sym=Fal Z_n=curve.Z_n, modes_R=curve.R_basis.modes[:, 2], modes_Z=curve.Z_basis.modes[:, 2], + NFP=NFP, + sym=curve.sym, name=name, ) From 01892431f7824171441870324d70fb8f3ddfca0f Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Tue, 9 Jul 2024 18:45:15 -0400 Subject: [PATCH 03/19] fix docstring and reorder args to match init --- desc/coils.py | 9 +++++++-- desc/geometry/curve.py | 6 +++++- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/desc/coils.py b/desc/coils.py index bdfb2363ac..fff929ca23 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -402,7 +402,7 @@ def __init__( super().__init__(current, R_n, Z_n, modes_R, modes_Z, NFP, sym, name) @classmethod - def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", name="", sym=False): + def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", sym=False, name=""): """Fit coordinates to FourierRZCoil representation. Parameters @@ -419,10 +419,15 @@ def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", name="", sym=Fal according to NFP. basis : {"rpz", "xyz"} basis for input coordinates. Defaults to "rpz" + sym : bool + Whether to enforce stellarator symmetry. + name : str + name for this coil + Returns ------- - coili : FourierRZCoil + coil : FourierRZCoil New representation of the coil parameterized by Fourier series for R,Z. """ diff --git a/desc/geometry/curve.py b/desc/geometry/curve.py index b928818cdc..6dc7e55af2 100644 --- a/desc/geometry/curve.py +++ b/desc/geometry/curve.py @@ -230,7 +230,7 @@ def from_input_file(cls, path): return curve @classmethod - def from_values(cls, coords, N=10, NFP=1, basis="rpz", name="", sym=False): + def from_values(cls, coords, N=10, NFP=1, basis="rpz", sym=False, name=""): """Fit coordinates to FourierRZCurve representation. Parameters @@ -245,6 +245,10 @@ def from_values(cls, coords, N=10, NFP=1, basis="rpz", name="", sym=False): according to NFP. basis : {"rpz", "xyz"} basis for input coordinates. Defaults to "rpz" + sym : bool + Whether to enforce stellarator symmetry. + name : str + Name for this curve. Returns ------- From 4ecae78edb3ebc68054434000c35cfb3a7a53c09 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Wed, 17 Jul 2024 16:46:36 -0600 Subject: [PATCH 04/19] add from_values for FourierPlanarCurve --- desc/coils.py | 64 ++++++++++++++---- desc/geometry/curve.py | 144 +++++++++++++++++++++++++++++++++++------ 2 files changed, 176 insertions(+), 32 deletions(-) diff --git a/desc/coils.py b/desc/coils.py index f42f40bb19..1200cc21dc 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -352,7 +352,7 @@ class FourierRZCoil(_Coil, FourierRZCurve): Parameters ---------- current : float - current through coil, in Amperes + Current through the coil, in Amperes. R_n, Z_n: array-like fourier coefficients for R, Z modes_R : array-like @@ -420,7 +420,7 @@ def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", sym=False, name= Parameters ---------- current : float - Current through the coil, in Amps. + Current through the coil, in Amperes. coords: ndarray, shape (num_coords,3) coordinates to fit a FourierRZCurve object with each column corresponding to xyz or rpz depending on the basis argument. @@ -443,9 +443,11 @@ def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", sym=False, name= New representation of the coil parameterized by Fourier series for R,Z. """ - curve = super().from_values(coords, N=N, NFP=NFP, basis=basis, sym=sym) + curve = super().from_values( + coords=coords, N=N, NFP=NFP, basis=basis, sym=sym, name=name + ) return cls( - current, + current=current, R_n=curve.R_n, Z_n=curve.Z_n, modes_R=curve.R_basis.modes[:, 2], @@ -462,7 +464,7 @@ class FourierXYZCoil(_Coil, FourierXYZCurve): Parameters ---------- current : float - current through coil, in Amperes + Current through the coil, in Amperes. X_n, Y_n, Z_n: array-like fourier coefficients for X, Y, Z modes : array-like @@ -525,7 +527,7 @@ def from_values(cls, current, coords, N=10, s=None, basis="xyz", name=""): Parameters ---------- current : float - Current through the coil, in Amps. + Current through the coil, in Amperes. coords: ndarray Coordinates to fit a FourierXYZCoil object with. N : int @@ -544,12 +546,13 @@ def from_values(cls, current, coords, N=10, s=None, basis="xyz", name=""): New representation of the coil parameterized by Fourier series for X,Y,Z. """ - curve = super().from_values(coords, N, s, basis) + curve = super().from_values(coords=coords, N=N, s=s, basis=basis, name=name) return cls( - current, + current=current, X_n=curve.X_n, Y_n=curve.Y_n, Z_n=curve.Z_n, + modes=curve.X_basis.modes[:, 2], name=name, ) @@ -627,6 +630,41 @@ def __init__( ): super().__init__(current, center, normal, r_n, modes, basis, name) + @classmethod + def from_values(cls, current, coords, N=10, basis="xyz", name=""): + """Fit coordinates to FourierPlanarCoil representation. + + Parameters + ---------- + current : float + Current through the coil, in Amperes. + coords: ndarray, shape (num_coords,3) + Coordinates to fit a FourierPlanarCurve object with each column + corresponding to xyz or rpz depending on the basis argument. + N : int + Fourier resolution of the new r representation. + basis : {"rpz", "xyz"} + Basis for input coordinates. Defaults to "xyz". + name : str + Name for this curve. + + Returns + ------- + curve : FourierPlanarCoil + New representation of the coil parameterized by a Fourier series for r. + + """ + curve = super().from_values(coords=coords, N=N, basis=basis, name=name) + return cls( + current=current, + center=curve.center, + normal=curve.normal, + r_n=curve.r_n, + modes=curve.r_basis.modes[:, 2], + basis="xyz", + name=name, + ) + class SplineXYZCoil(_Coil, SplineXYZCurve): """Coil parameterized by spline points in X,Y,Z. @@ -634,7 +672,7 @@ class SplineXYZCoil(_Coil, SplineXYZCurve): Parameters ---------- current : float - current through coil, in Amperes + Current through the coil, in Amperes. X, Y, Z: array-like Points for X, Y, Z describing the curve. If the endpoint is included (ie, X[0] == X[-1]), then the final point will be dropped. @@ -746,7 +784,7 @@ def from_values( Parameters ---------- current : float - Current through the coil, in Amps. + Current through the coil, in Amperes. coords: ndarray Points for X, Y, Z describing the curve. If the endpoint is included (ie, X[0] == X[-1]), then the final point will be dropped. @@ -778,9 +816,11 @@ def from_values( New representation of the coil parameterized by splines in X,Y,Z. """ - curve = super().from_values(coords, knots, method, basis=basis) + curve = super().from_values( + coords=coords, knots=knots, method=method, basis=basis, name=name + ) return cls( - current, + current=current, X=curve.X, Y=curve.Y, Z=curve.Z, diff --git a/desc/geometry/curve.py b/desc/geometry/curve.py index 0db7aa66ec..f347ec59b2 100644 --- a/desc/geometry/curve.py +++ b/desc/geometry/curve.py @@ -6,12 +6,13 @@ from desc.backend import jnp, put from desc.basis import FourierSeries -from desc.compute import rpz2xyz, xyz2rpz +from desc.compute import rpz2xyz, rpz2xyz_vec, xyz2rpz, xyz2rpz_vec +from desc.compute.geom_utils import rotation_matrix from desc.grid import LinearGrid from desc.io import InputReader from desc.optimizable import optimizable_parameter from desc.transform import Transform -from desc.utils import check_nonnegint, check_posint, copy_coeffs, errorif +from desc.utils import check_nonnegint, check_posint, copy_coeffs, errorif, warnif from .core import Curve @@ -230,7 +231,7 @@ def from_input_file(cls, path): return curve @classmethod - def from_values(cls, coords, N=10, NFP=1, basis="rpz", sym=False, name=""): + def from_values(cls, coords, N=10, NFP=1, sym=False, basis="rpz", name=""): """Fit coordinates to FourierRZCurve representation. Parameters @@ -243,10 +244,10 @@ def from_values(cls, coords, N=10, NFP=1, basis="rpz", sym=False, name=""): NFP : int Number of field periods, the curve will have a discrete toroidal symmetry according to NFP. - basis : {"rpz", "xyz"} - basis for input coordinates. Defaults to "rpz" sym : bool Whether to enforce stellarator symmetry. + basis : {"rpz", "xyz"} + basis for input coordinates. Defaults to "rpz" name : str Name for this curve. @@ -291,7 +292,15 @@ def from_values(cls, coords, N=10, NFP=1, basis="rpz", sym=False, name=""): transform = Transform(grid, basis, build_pinv=True) R_n = transform.fit(R) Z_n = transform.fit(Z) - return FourierRZCurve(R_n=R_n, Z_n=Z_n, NFP=NFP, name=name, sym=sym) + return cls( + R_n=R_n, + Z_n=Z_n, + modes_R=basis.modes[:, 2], + modes_Z=basis.modes[:, 2], + NFP=NFP, + sym=sym, + name=name, + ) def _unclose_curve(X, Y, Z): @@ -493,17 +502,18 @@ def from_values(cls, coords, N=10, s=None, basis="xyz", name=""): Parameters ---------- coords: ndarray, shape (num_coords,3) - coordinates to fit a FourierXYZCurve object, with each column + Coordinates to fit a FourierXYZCurve object, with each column corresponding to xyz or rpz depending on the basis argument. N : int Fourier resolution of the new X,Y,Z representation. s : ndarray or "arclength" - arbitrary curve parameter to use for the fitting. + Arbitrary curve parameter to use for the fitting. Should be monotonic, 1D array of same length as coords. if None, defaults linearly spaced in [0,2pi) Alternative, can pass "arclength" to use normalized distance between points. basis : {"rpz", "xyz"} - basis for input coordinates. Defaults to "xyz" + Basis for input coordinates. Defaults to "xyz". + Returns ------- curve : FourierXYZCurve @@ -551,7 +561,7 @@ def from_values(cls, coords, N=10, s=None, basis="xyz", name=""): X_n = transform.fit(X) Y_n = transform.fit(Y) Z_n = transform.fit(Z) - return FourierXYZCurve(X_n=X_n, Y_n=Y_n, Z_n=Z_n, name=name) + return cls(X_n=X_n, Y_n=Y_n, Z_n=Z_n, modes=basis.modes[:, 2], name=name) class FourierPlanarCurve(Curve): @@ -641,7 +651,7 @@ def center(self, new): else: raise ValueError( "center should be a 3 element vector in " - + self._basis + + self.basis + " coordinates, got {}".format(new) ) @@ -658,7 +668,7 @@ def normal(self, new): else: raise ValueError( "normal should be a 3 element vector in " - + self._basis + + self.basis + " coordinates, got {}".format(new) ) @@ -678,6 +688,23 @@ def r_n(self, new): + f"basis with {self.r_basis.num_modes} modes." ) + @property + def basis(self): + """Coordinate system for center and normal vectors.""" + return self._basis + + @basis.setter + def basis(self, new): + assert new.lower() in ["xyz", "rpz"] + if new != self.basis: + if new == "xyz": + self.normal = rpz2xyz_vec(self.normal, phi=self.center[1]) + self.center = rpz2xyz(self.center) + else: + self.center = xyz2rpz(self.center) + self.normal = xyz2rpz_vec(self.normal, phi=self.center[1]) + self._basis = new + def get_coeffs(self, n): """Get Fourier coefficients for given mode number(s).""" n = np.atleast_1d(n).astype(int) @@ -744,10 +771,82 @@ def compute( transforms=transforms, data=data, override_grid=override_grid, - basis_in=self._basis, + basis_in=self.basis, **kwargs, ) + @classmethod + def from_values(cls, coords, N=10, basis="xyz", name=""): + """Fit coordinates to FourierPlanarCurve representation. + + Parameters + ---------- + coords: ndarray, shape (num_coords,3) + Coordinates to fit a FourierPlanarCurve object with each column + corresponding to xyz or rpz depending on the basis argument. + N : int + Fourier resolution of the new r representation. + basis : {"rpz", "xyz"} + Basis for input coordinates. Defaults to "xyz". + name : str + Name for this curve. + + Returns + ------- + curve : FourierPlanarCurve + New representation of the curve parameterized by a Fourier series for r. + + """ + # convert to xyz basis + if basis == "rpz": + coords = rpz2xyz(coords) + coords = np.atleast_2d(coords) + + # center + center = np.mean(coords, axis=0) + coords = coords - center # shift to origin + + # normal + idx = coords.shape[0] // 2 + normal = np.cross(coords[0, :], coords[idx, :]) # 2 arbitrary points on curve + normal = normal / np.linalg.norm(normal) + + # axis and angle of rotation + Z_axis = np.array([0, 0, 1]) + axis = np.cross(Z_axis, normal) + axis = axis / np.linalg.norm(axis) + angle = np.arccos(np.dot(Z_axis, normal)) + rotmat = rotation_matrix(axis, angle) + coords = coords @ rotmat # rotate to X-Y plane + warnif( + np.max(np.abs(coords[:, 2])) > 1e-14, # check that Z=0 for all points + UserWarning, + "Curve values are not planar!", + ) + + # polar radius and angle + r = np.sqrt(coords[:, 0] ** 2 + coords[:, 1] ** 2) + s = np.arctan2(coords[:, 1], coords[:, 0]) + s = np.mod(s + 2 * np.pi, 2 * np.pi) # mod angle to range [0, 2*pi) + idx = np.argsort(s) # sort angle to be monotonically increasing + r = r[idx] + s = s[idx] + + # Fourier transform + basis = FourierSeries(N, NFP=1, sym=False) + grid_fit = LinearGrid(zeta=s, NFP=1) + transform_fit = Transform(grid_fit, basis, build_pinv=True) + r_n = transform_fit.fit(r) + + return cls( + center=center, + normal=normal, + r_n=r_n, + modes=basis.modes[:, 2], + basis="xyz", + name=name, + ) + class SplineXYZCurve(Curve): """Curve parameterized by spline knots in X,Y,Z. @@ -973,7 +1072,7 @@ def compute( ) @classmethod - def from_values(cls, coords, knots=None, method="cubic", name="", basis="xyz"): + def from_values(cls, coords, knots=None, method="cubic", basis="xyz", name=""): """Create SplineXYZCurve from coordinate values. Parameters @@ -984,7 +1083,7 @@ def from_values(cls, coords, knots=None, method="cubic", name="", basis="xyz"): endpoint is included (ie, X[0] == X[-1]), then the final point will be dropped. knots : ndarray - arbitrary curve parameter values to use for spline knots, + Arbitrary curve parameter values to use for spline knots, should be an 1D ndarray of same length as the input. (input length in this case is determined by grid argument, since the input coordinates come from @@ -1000,10 +1099,10 @@ def from_values(cls, coords, knots=None, method="cubic", name="", basis="xyz"): - `'cubic2'`: C2 cubic splines (aka natural splines) - `'catmull-rom'`: C1 cubic centripetal "tension" splines - name : str - name for this curve basis : {"rpz", "xyz"} - basis for input coordinates. Defaults to "xyz" + Basis for input coordinates. Defaults to "xyz". + name : str + Name for this curve. Returns ------- @@ -1013,6 +1112,11 @@ def from_values(cls, coords, knots=None, method="cubic", name="", basis="xyz"): """ if basis == "rpz": coords = rpz2xyz(coords) - return SplineXYZCurve( - coords[:, 0], coords[:, 1], coords[:, 2], knots, method, name + return cls( + X=coords[:, 0], + Y=coords[:, 1], + Z=coords[:, 2], + knots=knots, + method=method, + name=name, ) From fa9865582a9826728d8a024ee2b1301918f95ac5 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Thu, 18 Jul 2024 14:58:42 -0600 Subject: [PATCH 05/19] more robus planar curve fitting --- desc/geometry/curve.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/desc/geometry/curve.py b/desc/geometry/curve.py index f347ec59b2..a0dded1cb1 100644 --- a/desc/geometry/curve.py +++ b/desc/geometry/curve.py @@ -4,7 +4,7 @@ import numpy as np -from desc.backend import jnp, put +from desc.backend import jnp, put, sign from desc.basis import FourierSeries from desc.compute import rpz2xyz, rpz2xyz_vec, xyz2rpz, xyz2rpz_vec from desc.compute.geom_utils import rotation_matrix @@ -807,9 +807,13 @@ def from_values(cls, coords, N=10, basis="xyz", name=""): coords = coords - center # shift to origin # normal - idx = coords.shape[0] // 2 - normal = np.cross(coords[0, :], coords[idx, :]) # 2 arbitrary points on curve - normal = normal / np.linalg.norm(normal) + normal = np.empty((0, 3)) # avg of normals btw all pts in case it is non-planar + for idx in range(1, coords.shape[0], coords.shape[0] // 8): + norm = np.cross(coords[0, :], coords[idx, :]) + norm = norm / np.linalg.norm(norm) + sgn = sign(-int(all(sign(norm) != sign(normal[0, :])))) if idx > 1 else 1 + normal = np.vstack((normal, norm * sgn)) + normal = np.mean(normal, axis=0) # axis and angle of rotation Z_axis = np.array([0, 0, 1]) @@ -821,7 +825,7 @@ def from_values(cls, coords, N=10, basis="xyz", name=""): warnif( np.max(np.abs(coords[:, 2])) > 1e-14, # check that Z=0 for all points UserWarning, - "Curve values are not planar!", + "Curve values are not planar! Using the projection onto a plane.", ) # polar radius and angle From bf30b82fb7dab28df05b1265faf9696587775feb Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Thu, 18 Jul 2024 15:00:16 -0600 Subject: [PATCH 06/19] undo debugging change --- desc/geometry/curve.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/desc/geometry/curve.py b/desc/geometry/curve.py index a0dded1cb1..7f63594334 100644 --- a/desc/geometry/curve.py +++ b/desc/geometry/curve.py @@ -808,7 +808,7 @@ def from_values(cls, coords, N=10, basis="xyz", name=""): # normal normal = np.empty((0, 3)) # avg of normals btw all pts in case it is non-planar - for idx in range(1, coords.shape[0], coords.shape[0] // 8): + for idx in range(1, coords.shape[0]): norm = np.cross(coords[0, :], coords[idx, :]) norm = norm / np.linalg.norm(norm) sgn = sign(-int(all(sign(norm) != sign(normal[0, :])))) if idx > 1 else 1 From 9c849de3cd83bec0dd8aa4a01b1ad0b6969083c2 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Thu, 18 Jul 2024 21:58:13 -0600 Subject: [PATCH 07/19] be explicit about class --- desc/coils.py | 8 ++++---- desc/geometry/curve.py | 10 ++++++---- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/desc/coils.py b/desc/coils.py index 140dc72df7..436e72cc0b 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -448,7 +448,7 @@ def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", sym=False, name= curve = super().from_values( coords=coords, N=N, NFP=NFP, basis=basis, sym=sym, name=name ) - return cls( + return FourierRZCoil( current=current, R_n=curve.R_n, Z_n=curve.Z_n, @@ -549,7 +549,7 @@ def from_values(cls, current, coords, N=10, s=None, basis="xyz", name=""): """ curve = super().from_values(coords=coords, N=N, s=s, basis=basis, name=name) - return cls( + return FourierXYZCoil( current=current, X_n=curve.X_n, Y_n=curve.Y_n, @@ -657,7 +657,7 @@ def from_values(cls, current, coords, N=10, basis="xyz", name=""): """ curve = super().from_values(coords=coords, N=N, basis=basis, name=name) - return cls( + return FourierPlanarCoil( current=current, center=curve.center, normal=curve.normal, @@ -821,7 +821,7 @@ def from_values( curve = super().from_values( coords=coords, knots=knots, method=method, basis=basis, name=name ) - return cls( + return SplineXYZCoil( current=current, X=curve.X, Y=curve.Y, diff --git a/desc/geometry/curve.py b/desc/geometry/curve.py index 7f63594334..477d58e00e 100644 --- a/desc/geometry/curve.py +++ b/desc/geometry/curve.py @@ -292,7 +292,7 @@ def from_values(cls, coords, N=10, NFP=1, sym=False, basis="rpz", name=""): transform = Transform(grid, basis, build_pinv=True) R_n = transform.fit(R) Z_n = transform.fit(Z) - return cls( + return FourierRZCurve( R_n=R_n, Z_n=Z_n, modes_R=basis.modes[:, 2], @@ -561,7 +561,9 @@ def from_values(cls, coords, N=10, s=None, basis="xyz", name=""): X_n = transform.fit(X) Y_n = transform.fit(Y) Z_n = transform.fit(Z) - return cls(X_n=X_n, Y_n=Y_n, Z_n=Z_n, modes=basis.modes[:, 2], name=name) + return FourierXYZCurve( + X_n=X_n, Y_n=Y_n, Z_n=Z_n, modes=basis.modes[:, 2], name=name + ) class FourierPlanarCurve(Curve): @@ -842,7 +844,7 @@ def from_values(cls, coords, N=10, basis="xyz", name=""): transform_fit = Transform(grid_fit, basis, build_pinv=True) r_n = transform_fit.fit(r) - return cls( + return FourierPlanarCurve( center=center, normal=normal, r_n=r_n, @@ -1116,7 +1118,7 @@ def from_values(cls, coords, knots=None, method="cubic", basis="xyz", name=""): """ if basis == "rpz": coords = rpz2xyz(coords) - return cls( + return SplineXYZCurve( X=coords[:, 0], Y=coords[:, 1], Z=coords[:, 2], From 2b6b41b2ff87fc4178298a6be1341912afb405c1 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Sun, 21 Jul 2024 19:12:15 -0600 Subject: [PATCH 08/19] use SVD to compute normal vector --- desc/geometry/curve.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/desc/geometry/curve.py b/desc/geometry/curve.py index 477d58e00e..e659c9f31d 100644 --- a/desc/geometry/curve.py +++ b/desc/geometry/curve.py @@ -4,7 +4,7 @@ import numpy as np -from desc.backend import jnp, put, sign +from desc.backend import jnp, put from desc.basis import FourierSeries from desc.compute import rpz2xyz, rpz2xyz_vec, xyz2rpz, xyz2rpz_vec from desc.compute.geom_utils import rotation_matrix @@ -809,13 +809,8 @@ def from_values(cls, coords, N=10, basis="xyz", name=""): coords = coords - center # shift to origin # normal - normal = np.empty((0, 3)) # avg of normals btw all pts in case it is non-planar - for idx in range(1, coords.shape[0]): - norm = np.cross(coords[0, :], coords[idx, :]) - norm = norm / np.linalg.norm(norm) - sgn = sign(-int(all(sign(norm) != sign(normal[0, :])))) if idx > 1 else 1 - normal = np.vstack((normal, norm * sgn)) - normal = np.mean(normal, axis=0) + U, _, _ = np.linalg.svd(coords.T) + normal = U[:, -1].T # left singular vector of the least singular value # axis and angle of rotation Z_axis = np.array([0, 0, 1]) From efff9d1d6ea4312fbe48edfa7e72b5e4770c8ae1 Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Tue, 23 Jul 2024 11:52:29 -0400 Subject: [PATCH 09/19] add tests, allow normals that are parallel to Zaxis --- CHANGELOG.md | 9 ++++++ desc/coils.py | 67 ++++++++++++++++++++++++++++++++++++++++++ desc/compute/_curve.py | 48 ++++++++++++++++++++---------- desc/geometry/core.py | 41 ++++++++++++++++++++++++-- desc/geometry/curve.py | 14 +++++---- tests/test_coils.py | 26 +++++++++++++++- 6 files changed, 181 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 75d36322a2..2aa9571155 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,15 @@ Changelog ========= +New Features + +- adds from_values method that was present in FourierRZCurve but missing in FourierRZCoil +- also adds new from_values method for FourierPlanarCurve and FourierPlanarCoil + +Minor changes + +- modifies FourierPlanarCurve compute implementation to check if normal is already parallel to the Zaxis, in which case no rotation is needed. + v0.12.0 ------- diff --git a/desc/coils.py b/desc/coils.py index 06f7d57099..c1d3e089b7 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -345,6 +345,73 @@ def to_SplineXYZ(self, knots=None, grid=None, method="cubic", name=""): self.current, coords, knots=knots, method=method, name=name, basis="xyz" ) + def to_FourierRZ(self, N=None, grid=None, NFP=None, sym=False, name=""): + """Convert Coil to FourierRZCoil representation. + + Note that some types of coils may not be representable in this basis. + + Parameters + ---------- + N : int + Fourier resolution of the new R,Z representation. + grid : Grid, int or None + Grid used to evaluate curve coordinates on to fit with FourierRZCoil. + If an integer, uses that many equally spaced points. + NFP : int + Number of field periods, the coil will have a discrete toroidal symmetry + according to NFP. + sym : bool, optional + whether the curve is stellarator-symmetric or not. Default + is False. + name : str + name for this coil + + Returns + ------- + curve : FourierRZCoil + New representation of the coil parameterized by Fourier series for R,Z. + + """ + NFP = 1 or NFP + if grid is None: + grid = LinearGrid(N=2 * N + 1, NFP=NFP) + coords = self.compute("x", grid=grid, basis="xyz")["x"] + return FourierRZCoil.from_values( + self.current, coords, N=N, NFP=NFP, basis="xyz", sym=sym, name=name + ) + + def to_FourierPlanar(self, N=None, grid=None, name=""): + """Convert Coil to FourierPlanarCoil representation. + + Note that some types of coils may not be representable in this basis. + In this case, a least-squares fit will be done to find the + planar coil that best represents the coil. + + Parameters + ---------- + N : int + Fourier resolution of the new FourierPlanarCoil representation. + grid : Grid, int or None + Grid used to evaluate curve coordinates on to fit with FourierPlanarCoil. + If an integer, uses that many equally spaced points. + name : str + name for this coil + + Returns + ------- + coil : FourierPlanarCoil + New representation of the coil parameterized by Fourier series for + minor radius r in a plane specified by a center position and normal + vector. + + """ + if grid is None: + grid = LinearGrid(N=2 * N + 1) + coords = self.compute("x", grid=grid, basis="xyz")["x"] + return FourierPlanarCoil.from_values( + self.current, coords, N=N, basis="xyz", name=name + ) + class FourierRZCoil(_Coil, FourierRZCurve): """Coil parameterized by fourier series for R,Z in terms of toroidal angle phi. diff --git a/desc/compute/_curve.py b/desc/compute/_curve.py index 818da18092..4e4ddae949 100644 --- a/desc/compute/_curve.py +++ b/desc/compute/_curve.py @@ -205,10 +205,14 @@ def _x_FourierPlanarCurve(params, transforms, profiles, data, **kwargs): coords = jnp.array([X, Y, Z]).T # rotate into place Zaxis = jnp.array([0.0, 0.0, 1.0]) # 2D curve in X-Y plane has normal = +Z axis - axis = cross(Zaxis, normal) - angle = jnp.arccos(dot(Zaxis, safenormalize(normal))) - A = rotation_matrix(axis=axis, angle=angle) - coords = jnp.matmul(coords, A.T) + center + # check if we are already in the x-y plane, if so + # no need to rotate into it + unit_normal = safenormalize(normal) + if not jnp.allclose(unit_normal, Zaxis): + axis = cross(Zaxis, normal) + angle = jnp.arccos(dot(Zaxis, unit_normal)) + A = rotation_matrix(axis=axis, angle=angle) + coords = jnp.matmul(coords, A.T) + center coords = jnp.matmul(coords, params["rotmat"].reshape((3, 3)).T) + params["shift"] # convert back to rpz coords = xyz2rpz(coords) @@ -245,10 +249,14 @@ def _x_s_FourierPlanarCurve(params, transforms, profiles, data, **kwargs): coords = jnp.array([dX, dY, dZ]).T # rotate into place Zaxis = jnp.array([0.0, 0.0, 1.0]) # 2D curve in X-Y plane has normal = +Z axis - axis = cross(Zaxis, normal) - angle = jnp.arccos(dot(Zaxis, safenormalize(normal))) - A = rotation_matrix(axis=axis, angle=angle) - coords = jnp.matmul(coords, A.T) + # check if we are already in the x-y plane, if so + # no need to rotate into it + unit_normal = safenormalize(normal) + if not jnp.allclose(unit_normal, Zaxis): + axis = cross(Zaxis, normal) + angle = jnp.arccos(dot(Zaxis, unit_normal)) + A = rotation_matrix(axis=axis, angle=angle) + coords = jnp.matmul(coords, A.T) coords = jnp.matmul(coords, params["rotmat"].reshape((3, 3)).T) # convert back to rpz coords = xyz2rpz_vec(coords, phi=data["phi"]) @@ -290,10 +298,14 @@ def _x_ss_FourierPlanarCurve(params, transforms, profiles, data, **kwargs): coords = jnp.array([d2X, d2Y, d2Z]).T # rotate into place Zaxis = jnp.array([0.0, 0.0, 1.0]) # 2D curve in X-Y plane has normal = +Z axis - axis = cross(Zaxis, normal) - angle = jnp.arccos(dot(Zaxis, safenormalize(normal))) - A = rotation_matrix(axis=axis, angle=angle) - coords = jnp.matmul(coords, A.T) + # check if we are already in the x-y plane, if so + # no need to rotate into it + unit_normal = safenormalize(normal) + if not jnp.allclose(unit_normal, Zaxis): + axis = cross(Zaxis, normal) + angle = jnp.arccos(dot(Zaxis, unit_normal)) + A = rotation_matrix(axis=axis, angle=angle) + coords = jnp.matmul(coords, A.T) coords = jnp.matmul(coords, params["rotmat"].reshape((3, 3)).T) # convert back to rpz coords = xyz2rpz_vec(coords, phi=data["phi"]) @@ -342,10 +354,14 @@ def _x_sss_FourierPlanarCurve(params, transforms, profiles, data, **kwargs): coords = jnp.array([d3X, d3Y, d3Z]).T # rotate into place Zaxis = jnp.array([0.0, 0.0, 1.0]) # 2D curve in X-Y plane has normal = +Z axis - axis = cross(Zaxis, normal) - angle = jnp.arccos(dot(Zaxis, safenormalize(normal))) - A = rotation_matrix(axis=axis, angle=angle) - coords = jnp.matmul(coords, A.T) + # check if we are already in the x-y plane, if so + # no need to rotate into it + unit_normal = safenormalize(normal) + if not jnp.allclose(unit_normal, Zaxis): + axis = cross(Zaxis, normal) + angle = jnp.arccos(dot(Zaxis, unit_normal)) + A = rotation_matrix(axis=axis, angle=angle) + coords = jnp.matmul(coords, A.T) coords = jnp.matmul(coords, params["rotmat"].reshape((3, 3)).T) # convert back to rpz coords = xyz2rpz_vec(coords, phi=data["phi"]) diff --git a/desc/geometry/core.py b/desc/geometry/core.py index fd754dfe83..ef16519d9d 100644 --- a/desc/geometry/core.py +++ b/desc/geometry/core.py @@ -284,7 +284,7 @@ def to_SplineXYZ(self, knots=None, grid=None, method="cubic", name=""): coords, knots=knots, method=method, name=name, basis="xyz" ) - def to_FourierRZ(self, N=None, grid=None, NFP=None, name=""): + def to_FourierRZ(self, N=None, grid=None, NFP=None, sym=False, name=""): """Convert Curve to FourierRZCurve representation. Note that some types of curves may not be representable in this basis. @@ -299,6 +299,9 @@ def to_FourierRZ(self, N=None, grid=None, NFP=None, name=""): NFP : int Number of field periods, the curve will have a discrete toroidal symmetry according to NFP. + sym : bool, optional + whether the curve is stellarator-symmetric or not. Default + is False. name : str name for this curve @@ -314,7 +317,41 @@ def to_FourierRZ(self, N=None, grid=None, NFP=None, name=""): if grid is None: grid = LinearGrid(N=2 * N + 1, NFP=NFP) coords = self.compute("x", grid=grid, basis="xyz")["x"] - return FourierRZCurve.from_values(coords, N=N, NFP=NFP, basis="xyz", name=name) + return FourierRZCurve.from_values( + coords, N=N, NFP=NFP, basis="xyz", name=name, sym=sym + ) + + def to_FourierPlanar(self, N=None, grid=None, name=""): + """Convert Curve to FourierPlanarCurve representation. + + Note that some types of curves may not be representable in this basis. + In this case, a least-squares fit will be done to find the + planar curve that best represents the curve. + + Parameters + ---------- + N : int + Fourier resolution of the new FourierPlanarCurve representation. + grid : Grid, int or None + Grid used to evaluate curve coordinates on to fit with FourierPlanarCurve. + If an integer, uses that many equally spaced points. + name : str + name for this curve + + Returns + ------- + curve : FourierPlanarCurve + New representation of the curve parameterized by Fourier series for + minor radius r in a plane specified by a center position and normal + vector. + + """ + from .curve import FourierPlanarCurve + + if grid is None: + grid = LinearGrid(N=2 * N + 1) + coords = self.compute("x", grid=grid, basis="xyz")["x"] + return FourierPlanarCurve.from_values(coords, N=N, basis="xyz", name=name) class Surface(IOAble, Optimizable, ABC): diff --git a/desc/geometry/curve.py b/desc/geometry/curve.py index e659c9f31d..1c70e05a2d 100644 --- a/desc/geometry/curve.py +++ b/desc/geometry/curve.py @@ -814,11 +814,15 @@ def from_values(cls, coords, N=10, basis="xyz", name=""): # axis and angle of rotation Z_axis = np.array([0, 0, 1]) - axis = np.cross(Z_axis, normal) - axis = axis / np.linalg.norm(axis) - angle = np.arccos(np.dot(Z_axis, normal)) - rotmat = rotation_matrix(axis, angle) - coords = coords @ rotmat # rotate to X-Y plane + # check if we are already in the x-y plane, if so + # no need to rotate into it + if not np.allclose(normal, Z_axis): + axis = np.cross(Z_axis, normal) + axis = axis / np.linalg.norm(axis) + angle = np.arccos(np.dot(Z_axis, normal)) + rotmat = rotation_matrix(axis, angle) + coords = coords @ rotmat # rotate to X-Y plane + warnif( np.max(np.abs(coords[:, 2])) > 1e-14, # check that Z=0 for all points UserWarning, diff --git a/tests/test_coils.py b/tests/test_coils.py index 346e7a90b9..894d705300 100644 --- a/tests/test_coils.py +++ b/tests/test_coils.py @@ -204,13 +204,27 @@ def test_adding_MagneticField_with_Coil_or_CoilSet(self): def test_converting_coil_types(self): """Test conversions between coil representations.""" s = np.linspace(0, 2 * np.pi, 100, endpoint=False) - coil1 = FourierRZCoil(1e6, [0, 10, 2], [-2, 0, 0]) + coil1 = FourierRZCoil(1e6, [0, 10, 0], [0, 0, 0]) coil2 = coil1.to_FourierXYZ(s=s) coil3 = coil1.to_SplineXYZ(knots=s) + coil4 = coil1.to_FourierRZ(N=coil1.N) + coil5 = coil1.to_FourierPlanar(N=10, grid=LinearGrid(N=50)) + grid = LinearGrid(zeta=s) x1 = coil1.compute("x", grid=grid, basis="xyz")["x"] x2 = coil2.compute("x", grid=grid, basis="xyz")["x"] x3 = coil3.compute("x", grid=grid, basis="xyz")["x"] + x4 = coil4.compute("x", grid=grid, basis="xyz")["x"] + # find the arctan angle corresponding to the above pts for the + # planar coil + x5 = coil5.compute( + "x", + grid=LinearGrid( + zeta=np.arctan2(x1[:, 1] - coil5.center[1], x1[:, 0] - coil5.center[0]) + ), + basis="xyz", + )["x"] + B1 = coil1.compute_magnetic_field( np.zeros((1, 3)), source_grid=grid, basis="xyz" ) @@ -220,10 +234,20 @@ def test_converting_coil_types(self): B3 = coil3.compute_magnetic_field( np.zeros((1, 3)), source_grid=grid, basis="xyz" ) + B4 = coil4.compute_magnetic_field( + np.zeros((1, 3)), source_grid=grid, basis="xyz" + ) + B5 = coil5.compute_magnetic_field( + np.zeros((1, 3)), source_grid=grid, basis="xyz" + ) np.testing.assert_allclose(x1, x2, atol=1e-12) np.testing.assert_allclose(x1, x3, atol=1e-12) + np.testing.assert_allclose(x1, x4, atol=1e-12) + np.testing.assert_allclose(x1, x5, atol=1e-6) np.testing.assert_allclose(B1, B2, rtol=1e-8, atol=1e-8) np.testing.assert_allclose(B1, B3, rtol=1e-3, atol=1e-8) + np.testing.assert_allclose(B1, B4, rtol=1e-8, atol=1e-8) + np.testing.assert_allclose(B1, B5, rtol=1e-8, atol=1e-8) class TestCoilSet: From 0037d5097d6f453c6e9f56d83a2eb670b8c97e62 Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Tue, 23 Jul 2024 11:57:49 -0400 Subject: [PATCH 10/19] avoid error if no resolution is passed to to_XXX methods --- desc/coils.py | 4 ++-- desc/geometry/core.py | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/desc/coils.py b/desc/coils.py index c1d3e089b7..60b4cb7d09 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -345,7 +345,7 @@ def to_SplineXYZ(self, knots=None, grid=None, method="cubic", name=""): self.current, coords, knots=knots, method=method, name=name, basis="xyz" ) - def to_FourierRZ(self, N=None, grid=None, NFP=None, sym=False, name=""): + def to_FourierRZ(self, N=10, grid=None, NFP=None, sym=False, name=""): """Convert Coil to FourierRZCoil representation. Note that some types of coils may not be representable in this basis. @@ -380,7 +380,7 @@ def to_FourierRZ(self, N=None, grid=None, NFP=None, sym=False, name=""): self.current, coords, N=N, NFP=NFP, basis="xyz", sym=sym, name=name ) - def to_FourierPlanar(self, N=None, grid=None, name=""): + def to_FourierPlanar(self, N=10, grid=None, name=""): """Convert Coil to FourierPlanarCoil representation. Note that some types of coils may not be representable in this basis. diff --git a/desc/geometry/core.py b/desc/geometry/core.py index ef16519d9d..60b9cec6f6 100644 --- a/desc/geometry/core.py +++ b/desc/geometry/core.py @@ -209,7 +209,7 @@ def __repr__(self): + " (name={})".format(self.name) ) - def to_FourierXYZ(self, N=None, grid=None, s=None, name=""): + def to_FourierXYZ(self, N=10, grid=None, s=None, name=""): """Convert Curve to FourierXYZCurve representation. Parameters @@ -284,7 +284,7 @@ def to_SplineXYZ(self, knots=None, grid=None, method="cubic", name=""): coords, knots=knots, method=method, name=name, basis="xyz" ) - def to_FourierRZ(self, N=None, grid=None, NFP=None, sym=False, name=""): + def to_FourierRZ(self, N=10, grid=None, NFP=None, sym=False, name=""): """Convert Curve to FourierRZCurve representation. Note that some types of curves may not be representable in this basis. @@ -321,7 +321,7 @@ def to_FourierRZ(self, N=None, grid=None, NFP=None, sym=False, name=""): coords, N=N, NFP=NFP, basis="xyz", name=name, sym=sym ) - def to_FourierPlanar(self, N=None, grid=None, name=""): + def to_FourierPlanar(self, N=10, grid=None, name=""): """Convert Curve to FourierPlanarCurve representation. Note that some types of curves may not be representable in this basis. From 29b35aa3404a0bb55b12df35a22f4c0444ff05c6 Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Tue, 23 Jul 2024 12:09:08 -0400 Subject: [PATCH 11/19] undo unneeded changes in _curve compute --- desc/compute/_curve.py | 49 +++++++++++++++--------------------------- desc/geometry/curve.py | 12 ++++------- 2 files changed, 21 insertions(+), 40 deletions(-) diff --git a/desc/compute/_curve.py b/desc/compute/_curve.py index 4e4ddae949..3a36743284 100644 --- a/desc/compute/_curve.py +++ b/desc/compute/_curve.py @@ -205,14 +205,10 @@ def _x_FourierPlanarCurve(params, transforms, profiles, data, **kwargs): coords = jnp.array([X, Y, Z]).T # rotate into place Zaxis = jnp.array([0.0, 0.0, 1.0]) # 2D curve in X-Y plane has normal = +Z axis - # check if we are already in the x-y plane, if so - # no need to rotate into it - unit_normal = safenormalize(normal) - if not jnp.allclose(unit_normal, Zaxis): - axis = cross(Zaxis, normal) - angle = jnp.arccos(dot(Zaxis, unit_normal)) - A = rotation_matrix(axis=axis, angle=angle) - coords = jnp.matmul(coords, A.T) + center + axis = cross(Zaxis, normal) + angle = jnp.arccos(dot(Zaxis, safenormalize(normal))) + A = rotation_matrix(axis=axis, angle=angle) + coords = jnp.matmul(coords, A.T) + center coords = jnp.matmul(coords, params["rotmat"].reshape((3, 3)).T) + params["shift"] # convert back to rpz coords = xyz2rpz(coords) @@ -249,14 +245,10 @@ def _x_s_FourierPlanarCurve(params, transforms, profiles, data, **kwargs): coords = jnp.array([dX, dY, dZ]).T # rotate into place Zaxis = jnp.array([0.0, 0.0, 1.0]) # 2D curve in X-Y plane has normal = +Z axis - # check if we are already in the x-y plane, if so - # no need to rotate into it - unit_normal = safenormalize(normal) - if not jnp.allclose(unit_normal, Zaxis): - axis = cross(Zaxis, normal) - angle = jnp.arccos(dot(Zaxis, unit_normal)) - A = rotation_matrix(axis=axis, angle=angle) - coords = jnp.matmul(coords, A.T) + axis = cross(Zaxis, normal) + angle = jnp.arccos(dot(Zaxis, safenormalize(normal))) + A = rotation_matrix(axis=axis, angle=angle) + coords = jnp.matmul(coords, A.T) coords = jnp.matmul(coords, params["rotmat"].reshape((3, 3)).T) # convert back to rpz coords = xyz2rpz_vec(coords, phi=data["phi"]) @@ -298,14 +290,10 @@ def _x_ss_FourierPlanarCurve(params, transforms, profiles, data, **kwargs): coords = jnp.array([d2X, d2Y, d2Z]).T # rotate into place Zaxis = jnp.array([0.0, 0.0, 1.0]) # 2D curve in X-Y plane has normal = +Z axis - # check if we are already in the x-y plane, if so - # no need to rotate into it - unit_normal = safenormalize(normal) - if not jnp.allclose(unit_normal, Zaxis): - axis = cross(Zaxis, normal) - angle = jnp.arccos(dot(Zaxis, unit_normal)) - A = rotation_matrix(axis=axis, angle=angle) - coords = jnp.matmul(coords, A.T) + axis = cross(Zaxis, normal) + angle = jnp.arccos(dot(Zaxis, safenormalize(normal))) + A = rotation_matrix(axis=axis, angle=angle) + coords = jnp.matmul(coords, A.T) coords = jnp.matmul(coords, params["rotmat"].reshape((3, 3)).T) # convert back to rpz coords = xyz2rpz_vec(coords, phi=data["phi"]) @@ -354,14 +342,11 @@ def _x_sss_FourierPlanarCurve(params, transforms, profiles, data, **kwargs): coords = jnp.array([d3X, d3Y, d3Z]).T # rotate into place Zaxis = jnp.array([0.0, 0.0, 1.0]) # 2D curve in X-Y plane has normal = +Z axis - # check if we are already in the x-y plane, if so - # no need to rotate into it - unit_normal = safenormalize(normal) - if not jnp.allclose(unit_normal, Zaxis): - axis = cross(Zaxis, normal) - angle = jnp.arccos(dot(Zaxis, unit_normal)) - A = rotation_matrix(axis=axis, angle=angle) - coords = jnp.matmul(coords, A.T) + + axis = cross(Zaxis, normal) + angle = jnp.arccos(dot(Zaxis, safenormalize(normal))) + A = rotation_matrix(axis=axis, angle=angle) + coords = jnp.matmul(coords, A.T) coords = jnp.matmul(coords, params["rotmat"].reshape((3, 3)).T) # convert back to rpz coords = xyz2rpz_vec(coords, phi=data["phi"]) diff --git a/desc/geometry/curve.py b/desc/geometry/curve.py index 1c70e05a2d..e2c583ea6f 100644 --- a/desc/geometry/curve.py +++ b/desc/geometry/curve.py @@ -814,14 +814,10 @@ def from_values(cls, coords, N=10, basis="xyz", name=""): # axis and angle of rotation Z_axis = np.array([0, 0, 1]) - # check if we are already in the x-y plane, if so - # no need to rotate into it - if not np.allclose(normal, Z_axis): - axis = np.cross(Z_axis, normal) - axis = axis / np.linalg.norm(axis) - angle = np.arccos(np.dot(Z_axis, normal)) - rotmat = rotation_matrix(axis, angle) - coords = coords @ rotmat # rotate to X-Y plane + axis = np.cross(Z_axis, normal) + angle = np.arccos(np.dot(Z_axis, normal)) + rotmat = rotation_matrix(axis, angle) + coords = coords @ rotmat # rotate to X-Y plane warnif( np.max(np.abs(coords[:, 2])) > 1e-14, # check that Z=0 for all points From 7ba36b70b1ca49a11ec217c6e8b829f8a590828b Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Tue, 23 Jul 2024 12:10:44 -0400 Subject: [PATCH 12/19] remove white space --- desc/compute/_curve.py | 1 - 1 file changed, 1 deletion(-) diff --git a/desc/compute/_curve.py b/desc/compute/_curve.py index 3a36743284..818da18092 100644 --- a/desc/compute/_curve.py +++ b/desc/compute/_curve.py @@ -342,7 +342,6 @@ def _x_sss_FourierPlanarCurve(params, transforms, profiles, data, **kwargs): coords = jnp.array([d3X, d3Y, d3Z]).T # rotate into place Zaxis = jnp.array([0.0, 0.0, 1.0]) # 2D curve in X-Y plane has normal = +Z axis - axis = cross(Zaxis, normal) angle = jnp.arccos(dot(Zaxis, safenormalize(normal))) A = rotation_matrix(axis=axis, angle=angle) From 801832d56435d3ed709f483419a06297e18d7a93 Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Tue, 23 Jul 2024 12:14:52 -0400 Subject: [PATCH 13/19] modify test to test a non-centered coil --- tests/test_coils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_coils.py b/tests/test_coils.py index 894d705300..cc3be5b766 100644 --- a/tests/test_coils.py +++ b/tests/test_coils.py @@ -204,7 +204,7 @@ def test_adding_MagneticField_with_Coil_or_CoilSet(self): def test_converting_coil_types(self): """Test conversions between coil representations.""" s = np.linspace(0, 2 * np.pi, 100, endpoint=False) - coil1 = FourierRZCoil(1e6, [0, 10, 0], [0, 0, 0]) + coil1 = FourierRZCoil(1e6, [0, 10, 1], [0, 0, 0]) coil2 = coil1.to_FourierXYZ(s=s) coil3 = coil1.to_SplineXYZ(knots=s) coil4 = coil1.to_FourierRZ(N=coil1.N) From f35498603d75f56dec7e4ee092865d0486c34579 Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Tue, 23 Jul 2024 12:14:52 -0400 Subject: [PATCH 14/19] modify test to test a non-centered coil --- tests/test_coils.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_coils.py b/tests/test_coils.py index 894d705300..19cc5cfbaa 100644 --- a/tests/test_coils.py +++ b/tests/test_coils.py @@ -204,11 +204,11 @@ def test_adding_MagneticField_with_Coil_or_CoilSet(self): def test_converting_coil_types(self): """Test conversions between coil representations.""" s = np.linspace(0, 2 * np.pi, 100, endpoint=False) - coil1 = FourierRZCoil(1e6, [0, 10, 0], [0, 0, 0]) + coil1 = FourierRZCoil(1e6, [0, 10, 1], [0, 0, 0]) coil2 = coil1.to_FourierXYZ(s=s) coil3 = coil1.to_SplineXYZ(knots=s) coil4 = coil1.to_FourierRZ(N=coil1.N) - coil5 = coil1.to_FourierPlanar(N=10, grid=LinearGrid(N=50)) + coil5 = coil1.to_FourierPlanar(N=10) grid = LinearGrid(zeta=s) x1 = coil1.compute("x", grid=grid, basis="xyz")["x"] @@ -243,7 +243,7 @@ def test_converting_coil_types(self): np.testing.assert_allclose(x1, x2, atol=1e-12) np.testing.assert_allclose(x1, x3, atol=1e-12) np.testing.assert_allclose(x1, x4, atol=1e-12) - np.testing.assert_allclose(x1, x5, atol=1e-6) + np.testing.assert_allclose(x1, x5, atol=1e-12) np.testing.assert_allclose(B1, B2, rtol=1e-8, atol=1e-8) np.testing.assert_allclose(B1, B3, rtol=1e-3, atol=1e-8) np.testing.assert_allclose(B1, B4, rtol=1e-8, atol=1e-8) From 6f6ea23d032b2a49c295e7f0a31369d38fec51ec Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Tue, 23 Jul 2024 20:04:03 -0400 Subject: [PATCH 15/19] update changelog --- CHANGELOG.md | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2aa9571155..cd3e219fc2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,12 +3,8 @@ Changelog New Features -- adds from_values method that was present in FourierRZCurve but missing in FourierRZCoil -- also adds new from_values method for FourierPlanarCurve and FourierPlanarCoil - -Minor changes - -- modifies FourierPlanarCurve compute implementation to check if normal is already parallel to the Zaxis, in which case no rotation is needed. +- adds ``from_values`` method that was present in ``FourierRZCurve`` but missing in ``FourierRZCoil`` +- also adds new ``from_values`` method for ``FourierPlanarCurve`` and ``FourierPlanarCoil`` v0.12.0 From cbb6543733471315beb2ef12044e469d882bb588 Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Tue, 23 Jul 2024 20:10:43 -0400 Subject: [PATCH 16/19] remove NFP from grid of to_FourierRZ calls --- desc/coils.py | 2 +- desc/geometry/core.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/desc/coils.py b/desc/coils.py index 60b4cb7d09..b836899b71 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -374,7 +374,7 @@ def to_FourierRZ(self, N=10, grid=None, NFP=None, sym=False, name=""): """ NFP = 1 or NFP if grid is None: - grid = LinearGrid(N=2 * N + 1, NFP=NFP) + grid = LinearGrid(N=2 * N + 1) coords = self.compute("x", grid=grid, basis="xyz")["x"] return FourierRZCoil.from_values( self.current, coords, N=N, NFP=NFP, basis="xyz", sym=sym, name=name diff --git a/desc/geometry/core.py b/desc/geometry/core.py index 60b9cec6f6..b7715d905a 100644 --- a/desc/geometry/core.py +++ b/desc/geometry/core.py @@ -315,7 +315,7 @@ def to_FourierRZ(self, N=10, grid=None, NFP=None, sym=False, name=""): NFP = 1 or NFP if grid is None: - grid = LinearGrid(N=2 * N + 1, NFP=NFP) + grid = LinearGrid(N=2 * N + 1) coords = self.compute("x", grid=grid, basis="xyz")["x"] return FourierRZCurve.from_values( coords, N=N, NFP=NFP, basis="xyz", name=name, sym=sym From d0eb4e0ba35ca97dfb8ca0d50475b519e8138ca8 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Tue, 23 Jul 2024 20:45:39 -0600 Subject: [PATCH 17/19] add basis option for to_FourierPlanar --- desc/geometry/core.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/desc/geometry/core.py b/desc/geometry/core.py index b7715d905a..b665b5d32b 100644 --- a/desc/geometry/core.py +++ b/desc/geometry/core.py @@ -300,8 +300,7 @@ def to_FourierRZ(self, N=10, grid=None, NFP=None, sym=False, name=""): Number of field periods, the curve will have a discrete toroidal symmetry according to NFP. sym : bool, optional - whether the curve is stellarator-symmetric or not. Default - is False. + Whether the curve is stellarator-symmetric or not. Default is False. name : str name for this curve @@ -321,7 +320,7 @@ def to_FourierRZ(self, N=10, grid=None, NFP=None, sym=False, name=""): coords, N=N, NFP=NFP, basis="xyz", name=name, sym=sym ) - def to_FourierPlanar(self, N=10, grid=None, name=""): + def to_FourierPlanar(self, N=10, grid=None, basis="xyz", name=""): """Convert Curve to FourierPlanarCurve representation. Note that some types of curves may not be representable in this basis. @@ -335,6 +334,8 @@ def to_FourierPlanar(self, N=10, grid=None, name=""): grid : Grid, int or None Grid used to evaluate curve coordinates on to fit with FourierPlanarCurve. If an integer, uses that many equally spaced points. + basis : {'xyz', 'rpz'} + Coordinate system for center and normal vectors. Default = 'xyz'. name : str name for this curve @@ -351,7 +352,7 @@ def to_FourierPlanar(self, N=10, grid=None, name=""): if grid is None: grid = LinearGrid(N=2 * N + 1) coords = self.compute("x", grid=grid, basis="xyz")["x"] - return FourierPlanarCurve.from_values(coords, N=N, basis="xyz", name=name) + return FourierPlanarCurve.from_values(coords, N=N, basis=basis, name=name) class Surface(IOAble, Optimizable, ABC): From 448d23961b528b835b8a7e8e1b4bde988a07e4bc Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Tue, 23 Jul 2024 20:49:50 -0600 Subject: [PATCH 18/19] ocd formatting --- tests/test_coils.py | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/tests/test_coils.py b/tests/test_coils.py index 19cc5cfbaa..a481a76783 100644 --- a/tests/test_coils.py +++ b/tests/test_coils.py @@ -215,15 +215,10 @@ def test_converting_coil_types(self): x2 = coil2.compute("x", grid=grid, basis="xyz")["x"] x3 = coil3.compute("x", grid=grid, basis="xyz")["x"] x4 = coil4.compute("x", grid=grid, basis="xyz")["x"] - # find the arctan angle corresponding to the above pts for the - # planar coil - x5 = coil5.compute( - "x", - grid=LinearGrid( - zeta=np.arctan2(x1[:, 1] - coil5.center[1], x1[:, 0] - coil5.center[0]) - ), - basis="xyz", - )["x"] + grid = LinearGrid( + zeta=np.arctan2(x1[:, 1] - coil5.center[1], x1[:, 0] - coil5.center[0]) + ) # zeta = arctan angle instead of s for the planar coil for the same points + x5 = coil5.compute("x", grid=grid, basis="xyz")["x"] B1 = coil1.compute_magnetic_field( np.zeros((1, 3)), source_grid=grid, basis="xyz" @@ -240,6 +235,7 @@ def test_converting_coil_types(self): B5 = coil5.compute_magnetic_field( np.zeros((1, 3)), source_grid=grid, basis="xyz" ) + np.testing.assert_allclose(x1, x2, atol=1e-12) np.testing.assert_allclose(x1, x3, atol=1e-12) np.testing.assert_allclose(x1, x4, atol=1e-12) From 57688ef54853f6830064ff2282aa2be357787b94 Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Wed, 24 Jul 2024 12:17:10 -0400 Subject: [PATCH 19/19] add basis to coil to_planar method, make sure input coordinates match specified basis --- desc/coils.py | 8 +++++--- desc/geometry/core.py | 2 +- tests/test_coils.py | 4 ++-- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/desc/coils.py b/desc/coils.py index cd4dc0f79c..77f56c2188 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -381,7 +381,7 @@ def to_FourierRZ(self, N=10, grid=None, NFP=None, sym=False, name=""): self.current, coords, N=N, NFP=NFP, basis="xyz", sym=sym, name=name ) - def to_FourierPlanar(self, N=10, grid=None, name=""): + def to_FourierPlanar(self, N=10, grid=None, basis="xyz", name=""): """Convert Coil to FourierPlanarCoil representation. Note that some types of coils may not be representable in this basis. @@ -395,6 +395,8 @@ def to_FourierPlanar(self, N=10, grid=None, name=""): grid : Grid, int or None Grid used to evaluate curve coordinates on to fit with FourierPlanarCoil. If an integer, uses that many equally spaced points. + basis : {'xyz', 'rpz'} + Coordinate system for center and normal vectors. Default = 'xyz'. name : str name for this coil @@ -408,9 +410,9 @@ def to_FourierPlanar(self, N=10, grid=None, name=""): """ if grid is None: grid = LinearGrid(N=2 * N + 1) - coords = self.compute("x", grid=grid, basis="xyz")["x"] + coords = self.compute("x", grid=grid, basis=basis)["x"] return FourierPlanarCoil.from_values( - self.current, coords, N=N, basis="xyz", name=name + self.current, coords, N=N, basis=basis, name=name ) diff --git a/desc/geometry/core.py b/desc/geometry/core.py index b665b5d32b..a389eb4441 100644 --- a/desc/geometry/core.py +++ b/desc/geometry/core.py @@ -351,7 +351,7 @@ def to_FourierPlanar(self, N=10, grid=None, basis="xyz", name=""): if grid is None: grid = LinearGrid(N=2 * N + 1) - coords = self.compute("x", grid=grid, basis="xyz")["x"] + coords = self.compute("x", grid=grid, basis=basis)["x"] return FourierPlanarCurve.from_values(coords, N=N, basis=basis, name=name) diff --git a/tests/test_coils.py b/tests/test_coils.py index 739039c22c..8d8b9c7d4c 100644 --- a/tests/test_coils.py +++ b/tests/test_coils.py @@ -213,7 +213,7 @@ def test_converting_coil_types(self): coil2 = coil1.to_FourierXYZ(s=s) coil3 = coil1.to_SplineXYZ(knots=s) coil4 = coil1.to_FourierRZ(N=coil1.N) - coil5 = coil1.to_FourierPlanar(N=10) + coil5 = coil1.to_FourierPlanar(N=10, basis="rpz") grid = LinearGrid(zeta=s) x1 = coil1.compute("x", grid=grid, basis="xyz")["x"] @@ -248,7 +248,7 @@ def test_converting_coil_types(self): np.testing.assert_allclose(B1, B2, rtol=1e-8, atol=1e-8) np.testing.assert_allclose(B1, B3, rtol=1e-3, atol=1e-8) np.testing.assert_allclose(B1, B4, rtol=1e-8, atol=1e-8) - np.testing.assert_allclose(B1, B5, rtol=1e-8, atol=1e-8) + np.testing.assert_allclose(B1, B5, rtol=1e-6, atol=1e-7) class TestCoilSet: