diff --git a/src/simsopt/field/coilset.py b/src/simsopt/field/coilset.py index 363cdffb4..d2ed74f8c 100644 --- a/src/simsopt/field/coilset.py +++ b/src/simsopt/field/coilset.py @@ -64,7 +64,7 @@ def __init__(self, base_coils=None, coils=None, surface=None): else: if coils is not None: raise ValueError("If base_coils is None, coils must be None as well") - base_curves = self._circlecurves_around_surface(self._surface, nfp=self._surface.nfp, coils_per_period=10) + base_curves = self._circlecurves_around_surface(self._surface, coils_per_period=10) base_currents = [Current(1e5) for _ in base_curves] # fix the currents on the default coilset for current in base_currents: @@ -101,7 +101,7 @@ def for_surface(cls, surf, coil_current=1e5, coils_per_period=5, nfp=None, curre """ if nfp is None: nfp = surf.nfp - base_curves = CoilSet._circlecurves_around_surface(surf, coils_per_period=coils_per_period, nfp=nfp, **kwargs) + base_curves = CoilSet._circlecurves_around_surface(surf, coils_per_period=coils_per_period, **kwargs) base_currents = [Current(coil_current) for _ in base_curves] if current_constraint == "fix_one": base_currents[0].fix_all() @@ -132,49 +132,6 @@ def to_makegrid_file(self, filename): [coil.current for coil in self.coils], nfp=1) - - @classmethod - def for_spec_equil(cls, spec, coils_per_period=5, current_constraint="fix_all", **kwargs): - """ - Create a CoilSet for a given SPEC equilibrium. The coils are created using - :obj:`create_equally_spaced_curves` with the given parameters. - - Args: - spec: The SPEC object for which to create the coils - coils_per_period: the number of coils per field period - nfp: The number of field periods. - current_constraint: "fix_one" or "fix_all" or "free_all" - - Keyword Args (passed to the create_equally_spaced_curves function) - coils_per_period: The number of coils per field period - order: The order of the Fourier expansion - R0: major radius of a torus on which the initial coils are placed - R1: The radius of the coils - factor: if R0 or R1 are None, they are factor times - the first sine/cosine coefficient (factor > 1) - use_stellsym: Whether to use stellarator symmetry - """ - nfp = spec.nfp - use_stellsym = spec.stellsym - total_current = spec.poloidal_current_amperes - total_coil_number = coils_per_period * nfp * (1 + int(use_stellsym)) # only coils for half-period - if spec.freebound: - surface = spec.computational_boundary - else: - surface = spec.boundary - base_curves = CoilSet._circlecurves_around_surface(surface, coils_per_period=coils_per_period, nfp=nfp, use_stellsym=use_stellsym, **kwargs) - base_currents = [Current(total_current/total_coil_number) for _ in base_curves] - if current_constraint == "fix_one": - base_currents[0].fix_all() - elif current_constraint == "fix_all": - [base_current.fix_all() for base_current in base_currents] - elif current_constraint == "free_all": - pass - else: - raise ValueError("current_constraint must be 'fix_one', 'fix_all' or 'free_all'") - base_coils = [Coil(curv, curr) for (curv, curr) in zip(base_curves, base_currents)] - coils = coils_via_symmetries(base_curves, base_currents, nfp, stellsym=True) - return cls(base_coils, coils, surface) def reduce(self, target_function, nsv='nonzero'): """ @@ -228,16 +185,13 @@ def base_coils(self, base_coils): self.append_parent(coilparent) @staticmethod - def _circlecurves_around_surface(surf, nfp=None, coils_per_period=4, order=6, R0=None, R1=None, use_stellsym=None, factor=2.): + def _circlecurves_around_surface(surf, coils_per_period=4, order=6, R0=None, R1=None, use_stellsym=None, factor=2.): """ return a set of base curves for a surface using the surfaces properties where possible """ + nfp = surf.nfp if use_stellsym is None: use_stellsym = surf.stellsym - if nfp is None: - nfp = surf.nfp - elif nfp != surf.nfp: - raise ValueError("nfp must equal surf.nfp") if R0 is None: R0 = surf.to_RZFourier().get_rc(0, 0) if R1 is None: @@ -460,6 +414,8 @@ def __init__(self, coilset=None, nsv=None, s_diag=None, u_matrix=None, vh_matrix raise ValueError("If any of [s_diag, u_matrix, vh_matrix] are None, all must be None") if s_diag is None: + if nsv > coilset.dof_size: + raise ValueError("nsv must be equal to or smaller than the coilset's number of DOFs if initializing with None") s_diag = np.ones(nsv) u_matrix = np.eye(nsv) vh_matrix = np.eye(nsv) @@ -500,8 +456,6 @@ def from_function(cls, coilset, target_function, nsv='nonzero'): u_matrix, s_diag, vh_matrix = np.linalg.svd(jaccers) if nsv == 'nonzero': nsv = len(s_diag) - else: - assert isinstance(nsv, int), "nsv must be an integer or 'nonzero'" return cls(coilset, nsv, s_diag, u_matrix, vh_matrix, target_function) def recalculate_reduced_basis(self, target_function=None): @@ -634,53 +588,56 @@ def plot_singular_vector(self, n, eps=1e-4, show_delta_B=True, engine='mayavi', show_delta_B: whether to show the change in the magnetic field engine: the plotting engine to use. Defaults to 'mayavi' """ - from mayavi import mlab - if n > self.nsv: - raise ValueError("n must be smaller than the number of singular values") - singular_vector = np.array(self.rsv[n]) - - plotsurf = self.surface.copy(range=SurfaceRZFourier.RANGE_FULL_TORUS) - if show_delta_B: - bs = self.bs - initial_points = self.bs.get_points_cart_ref() - bs.set_points(plotsurf.gamma().reshape((-1, 3))) - - current_x = np.copy(self.x) - startpositions = [np.copy(coil.curve.gamma()) for coil in self.coilset.coils] - plot(self.coilset.coils, close=True, engine='mayavi',tube_radius=0.02, color=(0, 0, 0), show=False, **kwargs) - if show_delta_B: - startB = np.copy(np.sum(bs.B().reshape((plotsurf.quadpoints_phi.size, plotsurf.quadpoints_theta.size, 3)) * plotsurf.unitnormal()*-1, axis=2)) - startB = np.concatenate((startB, startB[:1, :]), axis=0) - startB = np.concatenate((startB, startB[:, :1]), axis=1) - - # Perturb the coils by the singular vector - self.coilset.x = self.coilset.x + singular_vector*eps - newpositions = [np.copy(coil.curve.gamma()) for coil in self.coilset.coils] - if show_delta_B: - changedB = np.copy(np.sum(bs.B().reshape((plotsurf.quadpoints_phi.size, plotsurf.quadpoints_theta.size, 3)) * plotsurf.unitnormal()*-1, axis=2)) - # close the plot - changedB = np.concatenate((changedB, changedB[:1, :]), axis=0) - changedB = np.concatenate((changedB, changedB[:, :1]), axis=1) - # plot the displacement vectors - for newcoilpos, startcoilpos in zip(newpositions, startpositions): - diffs = (0.05/eps) * (startcoilpos - newcoilpos) - x = startcoilpos[:, 0] - y = startcoilpos[:, 1] - z = startcoilpos[:, 2] - # enlarge the difference vectors for better visibility - dx = diffs[:, 0] - dy = diffs[:, 1] - dz = diffs[:, 2] - mlab.quiver3d(x, y, z, dx, dy, dz, line_width=4, **kwargs) - - if show_delta_B: - plot([plotsurf,], engine='mayavi', wireframe=False, close=True, scalars=changedB-startB, colormap='Reds', show=False, **kwargs) + if engine == 'mayavi': + from mayavi import mlab + if n > self.nsv: + raise ValueError("n must be smaller than the number of singular values") + singular_vector = np.array(self.rsv[n]) + + plotsurf = self.surface.copy(range=SurfaceRZFourier.RANGE_FULL_TORUS) + if show_delta_B: + bs = self.bs + initial_points = self.bs.get_points_cart_ref() + bs.set_points(plotsurf.gamma().reshape((-1, 3))) + + current_x = np.copy(self.x) + startpositions = [np.copy(coil.curve.gamma()) for coil in self.coilset.coils] + plot(self.coilset.coils, close=True, engine='mayavi',tube_radius=0.02, color=(0, 0, 0), show=False, **kwargs) + if show_delta_B: + startB = np.copy(np.sum(bs.B().reshape((plotsurf.quadpoints_phi.size, plotsurf.quadpoints_theta.size, 3)) * plotsurf.unitnormal()*-1, axis=2)) + startB = np.concatenate((startB, startB[:1, :]), axis=0) + startB = np.concatenate((startB, startB[:, :1]), axis=1) + + # Perturb the coils by the singular vector + self.coilset.x = self.coilset.x + singular_vector*eps + newpositions = [np.copy(coil.curve.gamma()) for coil in self.coilset.coils] + if show_delta_B: + changedB = np.copy(np.sum(bs.B().reshape((plotsurf.quadpoints_phi.size, plotsurf.quadpoints_theta.size, 3)) * plotsurf.unitnormal()*-1, axis=2)) + # close the plot + changedB = np.concatenate((changedB, changedB[:1, :]), axis=0) + changedB = np.concatenate((changedB, changedB[:, :1]), axis=1) + # plot the displacement vectors + for newcoilpos, startcoilpos in zip(newpositions, startpositions): + diffs = (0.05/eps) * (startcoilpos - newcoilpos) + x = startcoilpos[:, 0] + y = startcoilpos[:, 1] + z = startcoilpos[:, 2] + # enlarge the difference vectors for better visibility + dx = diffs[:, 0] + dy = diffs[:, 1] + dz = diffs[:, 2] + mlab.quiver3d(x, y, z, dx, dy, dz, line_width=4, **kwargs) + + if show_delta_B: + plot([plotsurf,], engine='mayavi', wireframe=False, close=True, scalars=changedB-startB, colormap='Reds', show=False, **kwargs) + else: + plot([plotsurf,], engine='mayavi', wireframe=False, close=True, colormap='Reds', show=False, **kwargs) + # plot the original coils again + self.x = current_x + plot(self.coilset.coils, close=True, tube_radius=0.02, engine='mayavi', color=(1, 1, 1), show=show, **kwargs) + # set bs set points back + if show_delta_B: + bs.set_points(initial_points) + return changedB, startB else: - plot([plotsurf,], engine='mayavi', wireframe=False, close=True, colormap='Reds', show=False, **kwargs) - # plot the original coils again - self.x = current_x - plot(self.coilset.coils, close=True, tube_radius=0.02, engine='mayavi', color=(1, 1, 1), show=show, **kwargs) - # set bs set points back - if show_delta_B: - bs.set_points(initial_points) - return changedB, startB \ No newline at end of file + raise ValueError("plotting a ReducedCoilSet is only implemented in Mayavi. you can access and plot the surface and coils attributes directly.") \ No newline at end of file diff --git a/src/simsopt/field/normal_field.py b/src/simsopt/field/normal_field.py index b4874b8c7..57ca43e1e 100644 --- a/src/simsopt/field/normal_field.py +++ b/src/simsopt/field/normal_field.py @@ -449,8 +449,6 @@ def get_vns_vnc_asarray(self, mpol=None, ntor=None): vns = self.get_vns_asarray(mpol, ntor) vnc = self.get_vnc_asarray(mpol, ntor) - if vnc is None: - vnc = np.zeros_like(vns) return vns, vnc def set_vns_asarray(self, vns, mpol=None, ntor=None): @@ -622,7 +620,7 @@ def target_function(coilset): cnf = CoilNormalField(coilset) # dummy CoilNormalField to evaluate the vnc and vnc output = cnf.vns.ravel()[coilset.surface.ntor+1:] #remove leading zeros if not coilset.surface.stellsym: - np.append(output, cnf.vnc.ravel()[coilset.surface.ntor:]) #remove leading zeros + output = np.append(output, cnf.vnc.ravel()[coilset.surface.ntor:]) #remove leading zeros return np.ravel(output) reduced_coilset = thiscoilset.reduce(target_function, nsv=nsv) @@ -672,7 +670,7 @@ def set_vns_vnc_asarray(self, *args, **kwargs): raise AttributeError('you cannot set fourier components, the coils do this!') def change_resolution(self, *args, **kwargs): - raise NotImplementedError('CoilNormalField.change_resolution() not implemented') + raise ValueError('CoilNormalField has no resolution, change parameters in its coilset') def fixed_range(self, *args, **kwargs): raise ValueError('no sense in fixing anything in a CoilNormalField') diff --git a/src/simsopt/geo/surfacerzfourier.py b/src/simsopt/geo/surfacerzfourier.py index 32429e93b..2b9a301dd 100644 --- a/src/simsopt/geo/surfacerzfourier.py +++ b/src/simsopt/geo/surfacerzfourier.py @@ -460,19 +460,21 @@ def copy(self, **kwargs): else: kwargs["quadpoints_phi"] = self.quadpoints_phi kwargs["quadpoints_theta"] = self.quadpoints_theta - elif quadpoints_theta is None: - if ntheta is not otherntheta or grid_range is not None: - kwargs["quadpoints_theta"] = Surface.get_theta_quadpoints(ntheta, range=grid_range) + else: + if quadpoints_theta is None: + if ntheta is not otherntheta or grid_range is not None: + kwargs["quadpoints_theta"] = Surface.get_theta_quadpoints(ntheta) + else: + kwargs["quadpoints_theta"] = self.quadpoints_theta else: - kwargs["quadpoints_theta"] = self.quadpoints_theta - elif quadpoints_phi is None: - if nphi is not othernphi or grid_range is not None: - kwargs["quadpoints_phi"] = Surface.get_phi_quadpoints(nfp, nphi, range=grid_range) + kwargs["quadpoints_theta"] = quadpoints_theta + if quadpoints_phi is None: + if nphi is not othernphi or grid_range is not None: + kwargs["quadpoints_phi"] = Surface.get_phi_quadpoints(nphi, range=grid_range, nfp=nfp) + else: + kwargs["quadpoints_phi"] = self.quadpoints_phi else: - kwargs["quadpoints_phi"] = self.quadpoints_phi - else: - kwargs["quadpoints_phi"] = quadpoints_phi - kwargs["quadpoints_theta"] = quadpoints_theta + kwargs["quadpoints_phi"] = quadpoints_phi # create new surface in old resolution surf = SurfaceRZFourier(mpol=self.mpol, ntor=self.ntor, nfp=nfp, stellsym=stellsym, **kwargs) diff --git a/tests/field/test_coilset.py b/tests/field/test_coilset.py index fd9d332c3..4032c3d61 100644 --- a/tests/field/test_coilset.py +++ b/tests/field/test_coilset.py @@ -19,6 +19,27 @@ def test_default_properties(self): coilset = CoilSet() self.assertEqual(len(coilset.base_coils), 10) self.assertEqual(len(coilset.coils), 20) + + def test_wrong_init(self): + curves, currents, _ = get_ncsx_data() + coils = [Coil(curve, current) for curve, current in zip(curves, currents)] + with self.assertRaises(ValueError): + CoilSet(coils = coils) + + def test_for_surface_classmethod(self): + s = SurfaceRZFourier(nfp=2, mpol=3, ntor=3) + coilset1 = CoilSet.for_surface(s, current_constraint='fix_all') + for coil in coilset1.base_coils: + self.assertEqual(coil.current.dof_size, 0) + coilset2 = CoilSet.for_surface(s, current_constraint='fix_one') + self.assertEqual(coilset2.base_coils[0].current.dof_size, 0) + for coil in coilset2.base_coils[1:]: + self.assertEqual(coil.current.dof_size, 1) + coilset3 = CoilSet.for_surface(s, current_constraint='free_all') + for coil in coilset3.base_coils: + self.assertEqual(coil.current.dof_size, 1) + + def test_to_from_mgrid(self): order = 25 @@ -50,15 +71,26 @@ def test_surface_setter_different_nfp(self): def test_surface_setter_nonstellsym(self): # Test the surface setter method - new_surface = SurfaceRZFourier(nfp=1, stellsym=False) - self.coilset.surface = new_surface - self.assertEqual(self.coilset.surface.deduced_range, SurfaceRZFourier.RANGE_FULL_TORUS) # for nfp==1 full torus is field period. + first_surface = SurfaceRZFourier(nfp=2, stellsym=True) + coilset = CoilSet(surface=first_surface) + self.assertEqual(coilset.surface.deduced_range, SurfaceRZFourier.RANGE_HALF_PERIOD) + second_surface = SurfaceRZFourier(nfp=2, stellsym=False) + coilset.surface = second_surface + self.assertEqual(coilset.surface.deduced_range, SurfaceRZFourier.RANGE_FIELD_PERIOD) def test_surface_setter_stellsym(self): # Test the surface setter method new_surface = SurfaceRZFourier(nfp=1, stellsym=True) self.coilset.surface = new_surface self.assertEqual(self.coilset.surface.deduced_range, SurfaceRZFourier.RANGE_HALF_PERIOD) + + def test_surface_setter_field_period(self): + s = SurfaceRZFourier(nfp=2, stellsym=False).copy(range='half period') + coilset = CoilSet(surface=s) + self.assertEqual(coilset.surface.deduced_range, SurfaceRZFourier.RANGE_FIELD_PERIOD) + s2 = SurfaceRZFourier(nfp=2, stellsym=True).copy(range='field period') + coilset = CoilSet(surface=s2) + self.assertEqual(coilset.surface.deduced_range, SurfaceRZFourier.RANGE_HALF_PERIOD) def test_base_coils(self): # Test the base_coils property @@ -118,6 +150,11 @@ def test_meansquared_curvature_penalty(self): # Test the meansquared_curvature_penalty function penalty = self.coilset.meansquared_curvature_penalty() self.assertIsNotNone(penalty.J()) + + def test_meansquared_curvature_threshold_penalty(self): + # Test the meansquared_curvature_penalty function + penalty = self.coilset.meansquared_curvature_threshold(CURVATURE_THRESHOLD=0.1) + self.assertIsNotNone(penalty.J()) def test_arc_length_variation_penalty(self): # Test the arc_length_variation_penalty function @@ -141,6 +178,10 @@ def test_coilset_to_vtk(self): self.coilset.to_vtk("test") self.assertTrue(os.path.exists("test_coils.vtu")) self.assertTrue(os.path.exists("test_surface.vts")) + with ScratchDir("."): + self.coilset.to_vtk("test2", add_biotsavart=False) + self.assertTrue(os.path.exists("test2_coils.vtu")) + self.assertTrue(os.path.exists("test2_surface.vts")) def test_save_load(self): with ScratchDir("."): @@ -174,12 +215,18 @@ def test_partially_empty_init(self): reduced_coilset = ReducedCoilSet() self.assertIsNotNone(reduced_coilset) with self.assertRaises(ValueError): - reduced_coilset.recalculate_reduced_basis() + reduced_coilset.recalculate_reduced_basis() #target function not set; fail reduced_coilset2 = ReducedCoilSet(self.unreduced_coilset) self.assertIsNotNone(reduced_coilset2) self.assertEqual(len(reduced_coilset2.x), len(self.unreduced_coilset.x)) reduced_coilset3 = ReducedCoilSet(self.unreduced_coilset, nsv=10) self.assertEqual(len(reduced_coilset3.x), 10) + with self.assertRaises(ValueError): + ReducedCoilSet(u_matrix=np.random.random((20,20))) + with self.assertRaises(ValueError): + ReducedCoilSet(self.unreduced_coilset, nsv=1000) + with self.assertRaises(TypeError): + ReducedCoilSet(self.unreduced_coilset, nsv=np.pi) def test_nsv_setter(self): self.coilset.nsv = 10 @@ -187,6 +234,11 @@ def test_nsv_setter(self): self.assertEqual(len(self.coilset.x), 10) self.coilset.recalculate_reduced_basis() self.assertEqual(len(self.coilset.x), 10) + with self.assertRaises(ValueError): + self.coilset.nsv = 1000 + with self.assertRaises(TypeError): + self.coilset.nsv = np.pi + self.coilset.nsv = 'nonzero' def test_surface_setter(self): with self.assertRaises(ValueError): @@ -220,6 +272,20 @@ def test_taylor_test_of_value_decomposition(self): new_function_value = self.test_target_function(self.coilset.coilset) function_diff = new_function_value - initial_function_value np.testing.assert_allclose(lsv, function_diff/(epsilon*singular_value), atol=1e-4) + + def test_wrong_setters(self): + with self.assertRaises(ValueError): + self.coilset.rsv = np.random.random(10) + with self.assertRaises(ValueError): + self.coilset.lsv = np.random.random(10) + with self.assertRaises(ValueError): + self.coilset.singular_values = np.random.random(10) + with self.assertRaises(ValueError): + self.coilset.coilset = self.unreduced_coilset + with self.assertRaises(ValueError): + self.coilset.coils = self.unreduced_coilset.coils + with self.assertRaises(ValueError): + self.set_dofs(np.random.random(self.coilset.nsv+1)) diff --git a/tests/field/test_normal_field.py b/tests/field/test_normal_field.py index dc81878af..3e55dd965 100644 --- a/tests/field/test_normal_field.py +++ b/tests/field/test_normal_field.py @@ -435,7 +435,7 @@ def test_reduce_coilset(self): epsilon = 1e-5 initial_vns = np.copy(cnf.vns) dofs = np.zeros(cnf.dof_size) - dofs[0] = 1e-5 + dofs[0] = epsilon cnf.x = dofs vnsdiff = initial_vns - cnf.vns vnsdiff_unraveled = vnsdiff.ravel()[ @@ -465,6 +465,8 @@ def test_inherited_methods_handled_correctly(self): cnf.get_index_in_dofs() with self.assertRaises(ValueError): cnf.fixed_range() + with self.assertRaises(ValueError): + cnf.change_resolution() def test_real_space_field(self): @@ -479,7 +481,108 @@ def test_real_space_field(self): thetasize, phisize = cnf.surface.quadpoints_theta.size, cnf.surface.quadpoints_phi.size directly_evaluated = np.copy(np.sum(cnf.coilset.bs.B().reshape(phisize, thetasize, 3) * cnf.surface.unitnormal(), axis=2)) #evaluate the field on the surface np.testing.assert_allclose(real_space_field, directly_evaluated, atol=1e-3) + + def test_cache_vns(self): + cnf = CoilNormalField() + self.assertIsNone(cnf._vns) + tmp = cnf.vns + np.testing.assert_equal(cnf._vns, tmp) + self.assertIsNotNone(cnf._vnc) + + def test_cache_vnc(self): + cnf = CoilNormalField() + self.assertIsNone(cnf._vnc) + tmp = cnf.vnc + np.testing.assert_equal(cnf._vnc, tmp) + self.assertIsNotNone(cnf._vns) + + def test_nonstellsym_reduce(self): + """ + nonstellaratorsymmetric fields have both vns and vnc components on the field + and a different function is used to reduce the coilset. + This test ensures that the reduction works for non-stellaratorsymmetric fields + """ + surface = SurfaceRZFourier(nfp=2, stellsym=False, mpol=6, ntor=6) + base_curves = CoilSet._circlecurves_around_surface( + surface, coils_per_period=2, order=4 + ) + base_coils = [Coil(curve, Current(1e5)) for curve in base_curves] + coilset = CoilSet(base_coils=base_coils, surface=surface) + cnf = CoilNormalField(coilset) + cnf.reduce_coilset() + + initial_vns = np.copy(cnf.vns) + initial_vnc = np.copy(cnf.vnc) + + epsilon = 1e-5 + dofs = np.zeros(cnf.dof_size); dofs[0] = epsilon; cnf.x = dofs + + vnsdiff = initial_vns - cnf.vns + vncdiff = initial_vnc - cnf.vnc + + diff_unraveled = np.concatenate( + (vnsdiff.ravel()[ + cnf.coilset.surface.ntor + 1: + ], + vncdiff.ravel()[ + cnf.coilset.surface.ntor: + ] + ) + ) + np.testing.assert_allclose( + cnf.coilset.lsv[0], + -1*diff_unraveled / (epsilon * cnf.coilset.singular_values[0]), + atol=1e-4, + ) + def test_vns_vnc_asarray(self): + cnf = CoilNormalField() + vns = cnf.get_vns_asarray() + vnc = cnf.get_vnc_asarray() + self.assertIsNotNone(vns) + self.assertIsNotNone(vnc) + + def test_wrong_index(self): + """ + test if accessing a wrong m or n raises an error + """ + surface = SurfaceRZFourier(nfp=3, stellsym=False, mpol=3, ntor=2) + base_curves = CoilSet._circlecurves_around_surface( + surface, coils_per_period=2, order=4 + ) + base_coils = [Coil(curve, Current(1e5)) for curve in base_curves] + coilset = CoilSet(base_coils=base_coils, surface=surface) + cnf = CoilNormalField(coilset) + with self.assertRaises(ValueError): + cnf.get_index_in_array(m=4, n=1) + with self.assertRaises(ValueError): + cnf.get_index_in_array(m=1, n=3) + with self.assertRaises(ValueError): + cnf.get_index_in_array(m=0, n=-1) + with self.assertRaises(ValueError): + cnf.get_index_in_array(m=3, n=-3) + with self.assertRaises(ValueError): + cnf.get_index_in_array(m=4, n=-3) + + def test_double_reduction(self): + """ + test if reducing a coilset twice works + """ + surface = SurfaceRZFourier(nfp=2, stellsym=False, mpol=6, ntor=6) + base_curves = CoilSet._circlecurves_around_surface( + surface, coils_per_period=2, order=4 + ) + base_coils = [Coil(curve, Current(1e5)) for curve in base_curves] + coilset = CoilSet(base_coils=base_coils, surface=surface) + cnf = CoilNormalField(coilset) + cnf.reduce_coilset() + num1 = int(cnf.dof_names[0][14]) # grab 'ReducedCoilSet*N* from first dof name + cnf.reduce_coilset() + num2 = int(cnf.dof_names[0][14]) # grab 'ReducedCoilSet*N* from first dof name + self.assertEqual(num2-num1, 1) # Coilset sucessfully replaced + + + diff --git a/tests/geo/test_surface_rzfourier.py b/tests/geo/test_surface_rzfourier.py index fce6fec00..1dc66467f 100755 --- a/tests/geo/test_surface_rzfourier.py +++ b/tests/geo/test_surface_rzfourier.py @@ -771,6 +771,24 @@ def test_fourier_transform_scalar(self): # Check that the result is the same as the original field: np.testing.assert_allclose(field/2*np.pi**2, field2) + def test_copy_method(self): + s = SurfaceRZFourier(mpol=4, ntor=5, nfp=3) + s2 = s.copy(quadpoints_phi=Surface.get_phi_quadpoints(nphi=100,range='field period')) + self.assertEqual(len(s2.quadpoints_phi), 100) + s3 = s.copy(quadpoints_theta=Surface.get_theta_quadpoints(ntheta=50)) + self.assertEqual(len(s3.quadpoints_theta), 50) + s4 = s.copy(ntheta = 42) + self.assertEqual(len(s4.quadpoints_theta), 42) + s5 = s.copy(nphi = 17) + self.assertEqual(len(s5.quadpoints_phi), 17) + s6 = s.copy(range='field period') + self.assertEqual(s6.deduced_range, Surface.RANGE_FIELD_PERIOD) + s7 = s.copy(nfp=10) + self.assertEqual(s7.nfp, 10) + s8 = s.copy(mpol=5, ntor=6) + self.assertEqual(s8.mpol, 5) + self.assertEqual(s8.ntor, 6) + class SurfaceRZPseudospectralTests(unittest.TestCase): def test_names(self): """