From 0b8a1d76f7d19b39a5fcd399a2b2ad1489e6e1ad Mon Sep 17 00:00:00 2001 From: Dario Panici Date: Mon, 18 Nov 2024 16:27:23 -0500 Subject: [PATCH] let linspace_XXX take in check_intersection --- desc/coils.py | 31 +++++++++++--- tests/test_objective_funs.py | 78 ++++++++++++++++++++++++------------ 2 files changed, 78 insertions(+), 31 deletions(-) diff --git a/desc/coils.py b/desc/coils.py index ea0da348a5..2de3d3d825 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -1633,7 +1633,14 @@ def compute_magnetic_vector_potential( @classmethod def linspaced_angular( - cls, coil, current=None, axis=[0, 0, 1], angle=2 * np.pi, n=10, endpoint=False + cls, + coil, + current=None, + axis=[0, 0, 1], + angle=2 * np.pi, + n=10, + endpoint=False, + check_intersection=True, ): """Create a CoilSet by repeating a coil at equal spacing around the torus. @@ -1651,6 +1658,8 @@ def linspaced_angular( Number of copies of original coil. endpoint : bool Whether to include a coil at final rotation angle. Default = False. + check_intersection : bool + whether to check the resulting coilsets for intersecting coils. """ assert isinstance(coil, _Coil) and not isinstance(coil, CoilSet) @@ -1664,11 +1673,17 @@ def linspaced_angular( coili.rotate(axis=axis, angle=phi[i]) coili.current = currents[i] coils.append(coili) - return cls(*coils) + return cls(*coils, check_intersection=check_intersection) @classmethod def linspaced_linear( - cls, coil, current=None, displacement=[2, 0, 0], n=4, endpoint=False + cls, + coil, + current=None, + displacement=[2, 0, 0], + n=4, + endpoint=False, + check_intersection=True, ): """Create a CoilSet by repeating a coil at equal spacing in a straight line. @@ -1685,6 +1700,8 @@ def linspaced_linear( Number of copies of original coil. endpoint : bool Whether to include a coil at final displacement location. Default = False. + check_intersection : bool + whether to check the resulting coilsets for intersecting coils. """ assert isinstance(coil, _Coil) and not isinstance(coil, CoilSet) @@ -1699,10 +1716,10 @@ def linspaced_linear( coili.translate(a[i] * displacement) coili.current = currents[i] coils.append(coili) - return cls(*coils) + return cls(*coils, check_intersection=check_intersection) @classmethod - def from_symmetry(cls, coils, NFP=1, sym=False): + def from_symmetry(cls, coils, NFP=1, sym=False, check_intersection=True): """Create a coil group by reflection and symmetry. Given coils over one field period, repeat coils NFP times between @@ -1721,6 +1738,8 @@ def from_symmetry(cls, coils, NFP=1, sym=False): sym : bool (optional) Whether to enforce stellarator symmetry. If True, the coils will be duplicated 2*NFP times. Default = False. + check_intersection : bool + whether to check the resulting coilsets for intersecting coils. Returns ------- @@ -1763,7 +1782,7 @@ def from_symmetry(cls, coils, NFP=1, sym=False): rotated_coils.rotate(axis=[0, 0, 1], angle=2 * jnp.pi * k / NFP) coilset += rotated_coils - return cls(*coilset) + return cls(*coilset, check_intersection=check_intersection) @classmethod def from_makegrid_coilfile(cls, coil_file, method="cubic", check_intersection=True): diff --git a/tests/test_objective_funs.py b/tests/test_objective_funs.py index a429358407..40c204466b 100644 --- a/tests/test_objective_funs.py +++ b/tests/test_objective_funs.py @@ -526,7 +526,7 @@ def test(eq): def test_boundary_error_biestsc(self): """Test calculation of boundary error using BIEST w/ sheet current.""" coil = FourierXYZCoil(5e5) - coilset = CoilSet.linspaced_angular(coil, n=100) + coilset = CoilSet.linspaced_angular(coil, n=100, check_intersection=False) coil_grid = LinearGrid(N=20) eq = Equilibrium(L=3, M=3, N=3, Psi=np.pi) eq.surface = FourierCurrentPotentialField.from_surface( @@ -548,7 +548,7 @@ def test_boundary_error_biestsc(self): def test_boundary_error_biest(self): """Test calculation of boundary error using BIEST.""" coil = FourierXYZCoil(5e5) - coilset = CoilSet.linspaced_angular(coil, n=100) + coilset = CoilSet.linspaced_angular(coil, n=100, check_intersection=False) coil_grid = LinearGrid(N=20) eq = Equilibrium(L=3, M=3, N=3, Psi=np.pi) eq.solve() @@ -565,7 +565,7 @@ def test_boundary_error_biest(self): def test_boundary_error_vacuum(self): """Test calculation of vacuum boundary error.""" coil = FourierXYZCoil(5e5) - coilset = CoilSet.linspaced_angular(coil, n=100) + coilset = CoilSet.linspaced_angular(coil, n=100, check_intersection=False) coil_grid = LinearGrid(N=20) eq = Equilibrium(L=3, M=3, N=3, Psi=np.pi) eq.solve() @@ -582,7 +582,7 @@ def test_boundary_error_vacuum(self): def test_boundary_error_nestor(self): """Test calculation of boundary error using NESTOR.""" coil = FourierXYZCoil(5e5) - coilset = CoilSet.linspaced_angular(coil, n=100) + coilset = CoilSet.linspaced_angular(coil, n=100, check_intersection=False) coil_grid = LinearGrid(N=20) eq = Equilibrium(L=3, M=3, N=3, Psi=np.pi) eq.solve() @@ -833,8 +833,12 @@ def test(coil, grid=None): assert len(f) == obj.dim_f coil = FourierPlanarCoil(r_n=2, basis="rpz") - coils = CoilSet.linspaced_linear(coil, n=3, displacement=[0, 3, 0]) - mixed_coils = MixedCoilSet.linspaced_linear(coil, n=2, displacement=[0, 7, 0]) + coils = CoilSet.linspaced_linear( + coil, n=3, displacement=[0, 3, 0], check_intersection=False + ) + mixed_coils = MixedCoilSet.linspaced_linear( + coil, n=2, displacement=[0, 7, 0], check_intersection=False + ) nested_coils = MixedCoilSet(coils, mixed_coils, check_intersection=False) grid = None # default grid @@ -856,8 +860,12 @@ def test(coil, grid=None): assert f.shape == (obj.dim_f,) coil = FourierPlanarCoil(r_n=2, basis="rpz") - coils = CoilSet.linspaced_linear(coil, n=3, displacement=[0, 3, 0]) - mixed_coils = MixedCoilSet.linspaced_linear(coil, n=2, displacement=[0, 7, 0]) + coils = CoilSet.linspaced_linear( + coil, n=3, displacement=[0, 3, 0], check_intersection=False + ) + mixed_coils = MixedCoilSet.linspaced_linear( + coil, n=2, displacement=[0, 7, 0], check_intersection=False + ) nested_coils = MixedCoilSet(coils, mixed_coils, check_intersection=False) grid = LinearGrid(N=5) # single grid @@ -879,8 +887,12 @@ def test(coil, grid=None): assert len(f) == obj.dim_f coil = FourierPlanarCoil(r_n=2, basis="rpz") - coils = CoilSet.linspaced_linear(coil, n=3, displacement=[0, 3, 0]) - mixed_coils = MixedCoilSet.linspaced_linear(coil, n=2, displacement=[0, 7, 0]) + coils = CoilSet.linspaced_linear( + coil, n=3, displacement=[0, 3, 0], check_intersection=False + ) + mixed_coils = MixedCoilSet.linspaced_linear( + coil, n=2, displacement=[0, 7, 0], check_intersection=False + ) nested_coils = MixedCoilSet(coils, mixed_coils, check_intersection=False) grid = [LinearGrid(N=5)] * 5 # single list of grids @@ -902,8 +914,12 @@ def test(coil, grid=None): assert len(f) == obj.dim_f coil = FourierPlanarCoil(r_n=2, basis="rpz") - coils = CoilSet.linspaced_linear(coil, n=3, displacement=[0, 3, 0]) - mixed_coils = MixedCoilSet.linspaced_linear(coil, n=2, displacement=[0, 7, 0]) + coils = CoilSet.linspaced_linear( + coil, n=3, displacement=[0, 3, 0], check_intersection=False + ) + mixed_coils = MixedCoilSet.linspaced_linear( + coil, n=2, displacement=[0, 7, 0], check_intersection=False + ) nested_coils = MixedCoilSet(coils, mixed_coils, check_intersection=False) grid = [[LinearGrid(N=5)] * 3, [LinearGrid(N=5)] * 2] # nested list of grids @@ -937,7 +953,9 @@ def test(coils, mindist, grid=None, expect_intersect=False, tol=None): n = 3 disp = 5 coil = FourierPlanarCoil(r_n=1, normal=[0, 0, 1]) - coils_linear = CoilSet.linspaced_linear(coil, n=n, displacement=[0, 0, disp]) + coils_linear = CoilSet.linspaced_linear( + coil, n=n, displacement=[0, 0, disp], check_intersection=False + ) test(coils_linear, disp / n) # planar toroidal coils, without symmetry @@ -946,7 +964,7 @@ def test(coils, mindist, grid=None, expect_intersect=False, tol=None): center = 3 r = 1 coil = FourierPlanarCoil(center=[center, 0, 0], normal=[0, 1, 0], r_n=r) - coils_angular = CoilSet.linspaced_angular(coil, n=4) + coils_angular = CoilSet.linspaced_angular(coil, n=4, check_intersection=False) test( coils_angular, np.sqrt(2) * (center - r), grid=LinearGrid(zeta=4), tol=1e-5 ) @@ -957,7 +975,9 @@ def test(coils, mindist, grid=None, expect_intersect=False, tol=None): center = 3 r = 1 coil = FourierPlanarCoil(center=[center, 0, 0], normal=[0, 1, 0], r_n=r) - coils = CoilSet.linspaced_angular(coil, angle=np.pi / 2, n=5, endpoint=True) + coils = CoilSet.linspaced_angular( + coil, angle=np.pi / 2, n=5, endpoint=True, check_intersection=False + ) coils_sym = CoilSet(coils[1::2], NFP=2, sym=True) test(coils_sym, 2 * (center - r) * np.sin(np.pi / 8), grid=LinearGrid(zeta=4)) @@ -965,10 +985,14 @@ def test(coils, mindist, grid=None, expect_intersect=False, tol=None): # TF coils instersect with the middle VF coil # extra coil is 5 m from middle VF coil tf_coil = FourierPlanarCoil(center=[2, 0, 0], normal=[0, 1, 0], r_n=1) - tf_coilset = CoilSet.linspaced_angular(tf_coil, n=4) + tf_coilset = CoilSet.linspaced_angular(tf_coil, n=4, check_intersection=False) vf_coil = FourierRZCoil(R_n=3, Z_n=-1) vf_coilset = CoilSet.linspaced_linear( - vf_coil, displacement=[0, 0, 2], n=3, endpoint=True + vf_coil, + displacement=[0, 0, 2], + n=3, + endpoint=True, + check_intersection=False, ) xyz_coil = FourierXYZCoil(X_n=[0, 6, 1], Y_n=[0, 0, 0], Z_n=[-1, 0, 0]) coils_mixed = MixedCoilSet( @@ -1029,7 +1053,7 @@ def test( ) eq = Equilibrium(surface=surf, NFP=1, M=2, N=0, sym=True) coil = FourierPlanarCoil(center=[R0, 0, 0], normal=[0, 1, 0], r_n=[a + offset]) - coils = CoilSet.linspaced_angular(coil, n=8) + coils = CoilSet.linspaced_angular(coil, n=8, check_intersection=False) test( eq, coils, @@ -1059,7 +1083,9 @@ def test( ) eq = Equilibrium(surface=surf, NFP=1, M=2, N=0, sym=True) coil = FourierPlanarCoil(center=[R0, 0, 0], normal=[0, 1, 0], r_n=[a + offset]) - coils = CoilSet.linspaced_angular(coil, angle=np.pi / 2, n=5, endpoint=True) + coils = CoilSet.linspaced_angular( + coil, angle=np.pi / 2, n=5, endpoint=True, check_intersection=False + ) coils = CoilSet(coils[1::2], NFP=2, sym=True) test( eq, @@ -1090,7 +1116,9 @@ def test( ) eq = Equilibrium(surface=surf, NFP=1, M=2, N=0, sym=True) coil = FourierPlanarCoil(center=[R0, 0, 0], normal=[0, 1, 0], r_n=[a + offset]) - coils = CoilSet.linspaced_angular(coil, angle=np.pi / 2, n=5, endpoint=True) + coils = CoilSet.linspaced_angular( + coil, angle=np.pi / 2, n=5, endpoint=True, check_intersection=False + ) coils = CoilSet(coils[1::2], NFP=2, sym=True) test( eq, @@ -1204,10 +1232,10 @@ def test_coil_linking_number(self): coil = FourierPlanarCoil(center=[10, 1, 0]) # regular modular coilset from symmetry, so that there are 10 coils, half going # one way and half going the other way - coilset = CoilSet.from_symmetry(coil, NFP=5, sym=True) + coilset = CoilSet.from_symmetry(coil, NFP=5, sym=True, check_intersection=False) coil2 = FourierRZCoil() # add a coil along the axis that links all the other coils - coilset2 = MixedCoilSet(coilset, coil2) + coilset2 = MixedCoilSet(coilset, coil2, check_intersection=False) obj = CoilSetLinkingNumber(coilset2) obj.build() @@ -1788,7 +1816,7 @@ def test(obj, eq, surface, d, print_init=False): def test_boundary_error_print(capsys): """Test that the boundary error objectives print correctly.""" coil = FourierXYZCoil(5e5) - coilset = CoilSet.linspaced_angular(coil, n=100) + coilset = CoilSet.linspaced_angular(coil, n=100, check_intersection=False) coil_grid = LinearGrid(N=20) eq = Equilibrium(L=3, M=3, N=3, Psi=np.pi) @@ -2667,7 +2695,7 @@ def test_compute_scalar_resolution_others(self, objective): def test_compute_scalar_resolution_coils(self, objective): """Coil objectives.""" coil = FourierXYZCoil() - coilset = CoilSet.linspaced_angular(coil) + coilset = CoilSet.linspaced_angular(coil, check_intersection=False) f = np.zeros_like(self.res_array, dtype=float) for i, res in enumerate(self.res_array): obj = ObjectiveFunction( @@ -2904,7 +2932,7 @@ def test_objective_no_nangrad(self, objective): def test_objective_no_nangrad_coils(self, objective): """Coil objectives.""" coil = FourierXYZCoil() - coilset = CoilSet.linspaced_angular(coil, n=3) + coilset = CoilSet.linspaced_angular(coil, n=3, check_intersection=False) obj = ObjectiveFunction(objective(coilset), use_jit=False) obj.build(verbose=0) g = obj.grad(obj.x())