diff --git a/CHANGELOG.md b/CHANGELOG.md index 75d36322a2..cd3e219fc2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,11 @@ 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`` + v0.12.0 ------- diff --git a/desc/coils.py b/desc/coils.py index efea3b0ee7..49c56d5968 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -346,6 +346,75 @@ 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=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. + + 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) + 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=10, grid=None, basis="xyz", 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. + basis : {'xyz', 'rpz'} + Coordinate system for center and normal vectors. Default = 'xyz'. + 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=basis)["x"] + return FourierPlanarCoil.from_values( + self.current, coords, N=N, basis=basis, name=name + ) + class FourierRZCoil(_Coil, FourierRZCurve): """Coil parameterized by fourier series for R,Z in terms of toroidal angle phi. @@ -353,7 +422,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 @@ -414,6 +483,50 @@ 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", sym=False, name=""): + """Fit coordinates to FourierRZCoil representation. + + Parameters + ---------- + current : float + 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. + 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" + sym : bool + Whether to enforce stellarator symmetry. + name : str + name for this coil + + + Returns + ------- + coil : FourierRZCoil + New representation of the coil parameterized by Fourier series for R,Z. + + """ + curve = super().from_values( + coords=coords, N=N, NFP=NFP, basis=basis, sym=sym, name=name + ) + return FourierRZCoil( + current=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], + NFP=NFP, + sym=curve.sym, + name=name, + ) + class FourierXYZCoil(_Coil, FourierXYZCurve): """Coil parameterized by fourier series for X,Y,Z in terms of arbitrary angle s. @@ -421,7 +534,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 @@ -484,7 +597,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 @@ -503,12 +616,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) - return cls( - current, + curve = super().from_values(coords=coords, N=N, s=s, basis=basis, name=name) + return FourierXYZCoil( + 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, ) @@ -586,6 +700,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 FourierPlanarCoil( + 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. @@ -593,7 +742,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. @@ -705,7 +854,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. @@ -737,9 +886,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) - return cls( - current, + curve = super().from_values( + coords=coords, knots=knots, method=method, basis=basis, name=name + ) + return SplineXYZCoil( + current=current, X=curve.X, Y=curve.Y, Z=curve.Z, diff --git a/desc/geometry/core.py b/desc/geometry/core.py index fd754dfe83..a389eb4441 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, 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. @@ -299,6 +299,8 @@ 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 @@ -312,9 +314,45 @@ def to_FourierRZ(self, N=None, grid=None, NFP=None, 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) + return FourierRZCurve.from_values( + coords, N=N, NFP=NFP, basis="xyz", name=name, sym=sym + ) + + 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. + 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. + basis : {'xyz', 'rpz'} + Coordinate system for center and normal vectors. Default = 'xyz'. + 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=basis)["x"] + return FourierPlanarCurve.from_values(coords, N=N, basis=basis, name=name) class Surface(IOAble, Optimizable, ABC): diff --git a/desc/geometry/curve.py b/desc/geometry/curve.py index a62ed8fdab..e2c583ea6f 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", name="", sym=False): + def from_values(cls, coords, N=10, NFP=1, sym=False, basis="rpz", name=""): """Fit coordinates to FourierRZCurve representation. Parameters @@ -243,8 +244,12 @@ def from_values(cls, coords, N=10, NFP=1, basis="rpz", name="", sym=False): NFP : int Number of field periods, the curve will have a discrete toroidal symmetry according to NFP. + sym : bool + Whether to enforce stellarator symmetry. basis : {"rpz", "xyz"} basis for input coordinates. Defaults to "rpz" + name : str + Name for this curve. Returns ------- @@ -287,7 +292,15 @@ def from_values(cls, coords, N=10, NFP=1, basis="rpz", name="", sym=False): 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 FourierRZCurve( + 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): @@ -489,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 @@ -547,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 FourierXYZCurve(X_n=X_n, Y_n=Y_n, Z_n=Z_n, 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): @@ -637,7 +653,7 @@ def center(self, new): else: raise ValueError( "center should be a 3 element vector in " - + self._basis + + self.basis + " coordinates, got {}".format(new) ) @@ -654,7 +670,7 @@ def normal(self, new): else: raise ValueError( "normal should be a 3 element vector in " - + self._basis + + self.basis + " coordinates, got {}".format(new) ) @@ -674,6 +690,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) @@ -740,10 +773,81 @@ 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 + 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]) + 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 + UserWarning, + "Curve values are not planar! Using the projection onto a plane.", + ) + + # 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 FourierPlanarCurve( + 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. @@ -969,7 +1073,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 @@ -980,7 +1084,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 @@ -996,10 +1100,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 ------- @@ -1010,5 +1114,10 @@ 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 + X=coords[:, 0], + Y=coords[:, 1], + Z=coords[:, 2], + knots=knots, + method=method, + name=name, ) diff --git a/tests/test_coils.py b/tests/test_coils.py index defa278c32..8d8b9c7d4c 100644 --- a/tests/test_coils.py +++ b/tests/test_coils.py @@ -209,13 +209,22 @@ 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, 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, basis="rpz") + 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"] + 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" ) @@ -225,10 +234,21 @@ 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-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) + np.testing.assert_allclose(B1, B5, rtol=1e-6, atol=1e-7) class TestCoilSet: