Skip to content

Commit

Permalink
Let linspace_XXX take in check_intersection (#1401)
Browse files Browse the repository at this point in the history
Resolves #1400
  • Loading branch information
dpanici authored Nov 19, 2024
2 parents 0b52a24 + 0b8a1d7 commit ff9d414
Show file tree
Hide file tree
Showing 2 changed files with 78 additions and 31 deletions.
31 changes: 25 additions & 6 deletions desc/coils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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)
Expand All @@ -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.
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
-------
Expand Down Expand Up @@ -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):
Expand Down
78 changes: 53 additions & 25 deletions tests/test_objective_funs.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand All @@ -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()
Expand All @@ -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()
Expand All @@ -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()
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
)
Expand All @@ -957,18 +975,24 @@ 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))

# mixture of toroidal field coilset, vertical field coilset, and extra coil
# 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(
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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()
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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())
Expand Down

0 comments on commit ff9d414

Please sign in to comment.