From 4cbe3ff2cc83d8eadb76961c9c2929431a9e55f5 Mon Sep 17 00:00:00 2001 From: Issra Ali Date: Fri, 10 Jul 2026 13:36:22 -0400 Subject: [PATCH 1/9] Added support for eq_fixed option in BoundaryError --- desc/objectives/_free_boundary.py | 101 +++++++++++++++++++++++------- 1 file changed, 78 insertions(+), 23 deletions(-) diff --git a/desc/objectives/_free_boundary.py b/desc/objectives/_free_boundary.py index b62603289d..5d8c0489ff 100644 --- a/desc/objectives/_free_boundary.py +++ b/desc/objectives/_free_boundary.py @@ -478,6 +478,7 @@ def __init__( eval_grid=None, field_grid=None, field_fixed=False, + eq_fixed=False, name="Boundary error", jac_chunk_size=None, *, @@ -487,12 +488,20 @@ def __init__( ): if target is None and bounds is None: target = 0 + self._eq = eq self._source_grid = source_grid self._eval_grid = eval_grid self._st, self._sz = s if isinstance(s, (tuple, list)) else (s, s) self._q = q self._field = [field] if not isinstance(field, list) else field self._field_grid = field_grid + errorif( + field_fixed and eq_fixed, + ValueError, + "At least one of eq_fixed or field_fixed must be false.", + ) + self._field_fixed = field_fixed + self._eq_fixed = eq_fixed self._bs_chunk_size = bs_chunk_size B_plasma_chunk_size = parse_argname_change( B_plasma_chunk_size, kwargs, "loop", "B_plasma_chunk_size" @@ -501,7 +510,9 @@ def __init__( B_plasma_chunk_size = None self._B_plasma_chunk_size = B_plasma_chunk_size self._sheet_current = hasattr(eq.surface, "Phi_mn") - things = [eq] + things = [] + if not eq_fixed: + things.append(self._eq) if not field_fixed: things.append(self._field) super().__init__( @@ -628,6 +639,9 @@ def build(self, use_jit=True, verbose=1): neq = 3 if self._sheet_current else 2 # number of equations we're using + if self._eq_fixed: + eq = self._eq + self._constants = { "eval_transforms": eval_transforms, "eval_profiles": eval_profiles, @@ -651,6 +665,39 @@ def build(self, use_jit=True, verbose=1): ) ) + if self._eq_fixed: + source_data = compute_fun( + "desc.equilibrium.equilibrium.Equilibrium", + self._eq_data_keys, + params=self._eq.params_dict, + transforms=source_transforms, + profiles=source_profiles, + ) + eval_data = ( + source_data + if self._use_same_grid + else compute_fun( + "desc.equilibrium.equilibrium.Equilibrium", + self._eq_data_keys, + params=self._eq.params_dict, + transforms=eval_transforms, + profiles=eval_profiles, + ) + ) + + Bplasma = virtual_casing_biot_savart( + eval_data, + source_data, + interpolator, + chunk_size=self._B_plasma_chunk_size, + ) + # need extra factor of B/2 bc we're evaluating on plasma surface + Bplasma = Bplasma + eval_data["B"] / 2 + + self._constants["source_data"]=source_data + self._constants["eval_data"]=eval_data + self._constants["Bplasma"]=Bplasma + timer.stop("Precomputing transforms") if verbose > 1: timer.disp("Precomputing transforms") @@ -698,24 +745,40 @@ def compute(self, eq_params, *field_params, constants=None): if field_params == (): # common case for field_fixed=True field_params = None constants = self._get_deprecated_constants(constants) - source_data = compute_fun( - "desc.equilibrium.equilibrium.Equilibrium", - self._eq_data_keys, - params=eq_params, - transforms=constants["source_transforms"], - profiles=constants["source_profiles"], - ) - eval_data = ( - source_data - if self._use_same_grid - else compute_fun( + + if self._eq_fixed: + source_data = constants["source_data"] + eval_data = constants["eval_data"] + Bplasma = constants["source_data"] + else: + source_data = compute_fun( "desc.equilibrium.equilibrium.Equilibrium", self._eq_data_keys, params=eq_params, - transforms=constants["eval_transforms"], - profiles=constants["eval_profiles"], + transforms=constants["source_transforms"], + profiles=constants["source_profiles"], ) - ) + eval_data = ( + source_data + if self._use_same_grid + else compute_fun( + "desc.equilibrium.equilibrium.Equilibrium", + self._eq_data_keys, + params=eq_params, + transforms=constants["eval_transforms"], + profiles=constants["eval_profiles"], + ) + ) + + Bplasma = virtual_casing_biot_savart( + eval_data, + source_data, + constants["interpolator"], + chunk_size=self._B_plasma_chunk_size, + ) + # need extra factor of B/2 bc we're evaluating on plasma surface + Bplasma = Bplasma + eval_data["B"] / 2 + if self._sheet_current: p = "desc.magnetic_fields._current_potential.FourierCurrentPotentialField" sheet_params = { @@ -745,14 +808,6 @@ def compute(self, eq_params, *field_params, constants=None): ) source_data["K_vc"] += sheet_source_data["K"] - Bplasma = virtual_casing_biot_savart( - eval_data, - source_data, - constants["interpolator"], - chunk_size=self._B_plasma_chunk_size, - ) - # need extra factor of B/2 bc we're evaluating on plasma surface - Bplasma = Bplasma + eval_data["B"] / 2 x = jnp.array([eval_data["R"], eval_data["phi"], eval_data["Z"]]).T # can always pass in field params. If they're None, it just uses the # defaults for the given field. From fd5eb5f6ee9f725cbc170512f0c0f98499b49a8d Mon Sep 17 00:00:00 2001 From: Issra Ali Date: Fri, 10 Jul 2026 14:59:12 -0400 Subject: [PATCH 2/9] Fixed a bug and moved some of the sheet_current stuff to be precomputed --- desc/objectives/_free_boundary.py | 87 ++++++++++++++++++++++--------- 1 file changed, 62 insertions(+), 25 deletions(-) diff --git a/desc/objectives/_free_boundary.py b/desc/objectives/_free_boundary.py index 5d8c0489ff..23c1bdc2b8 100644 --- a/desc/objectives/_free_boundary.py +++ b/desc/objectives/_free_boundary.py @@ -698,6 +698,38 @@ def build(self, use_jit=True, verbose=1): self._constants["eval_data"]=eval_data self._constants["Bplasma"]=Bplasma + # sheet current stuff + if self._sheet_current: + p = "desc.magnetic_fields._current_potential.FourierCurrentPotentialField" + sheet_params = { + "R_lmn": self._eq.params_dict["Rb_lmn"], + "Z_lmn": self._eq.params_dict["Zb_lmn"], + "I": self._eq.params_dict["I"], + "G": self._eq.params_dict["G"], + "Phi_mn": self._eq.params_dict["Phi_mn"], + } + sheet_source_data = compute_fun( + p, + self._sheet_data_keys, + params=sheet_params, + transforms=self._constants["sheet_source_transforms"], + profiles={}, + ) + sheet_eval_data = ( + sheet_source_data + if self._use_same_grid + else compute_fun( + p, + self._sheet_data_keys, + params=sheet_params, + transforms=self._constants["sheet_eval_transforms"], + profiles={}, + ) + ) + + self._constants["sheet_eval_data"]=sheet_eval_data + source_data["K_vc"] += sheet_source_data["K"] + timer.stop("Precomputing transforms") if verbose > 1: timer.disp("Precomputing transforms") @@ -727,7 +759,7 @@ def compute(self, eq_params, *field_params, constants=None): Parameters ---------- eq_params : dict - Dictionary of equilibrium degrees of freedom, eg Equilibrium.params_dict + Dictionary of equilibrium degrees of freedom, eg Equilibrium.params_dict, if equilibrium is not fixed field_params : dict Dictionary of field parameters, if field is not fixed. constants : dict @@ -744,12 +776,17 @@ def compute(self, eq_params, *field_params, constants=None): """ if field_params == (): # common case for field_fixed=True field_params = None + if eq_params == (): + eq_params = None + constants = self._get_deprecated_constants(constants) if self._eq_fixed: source_data = constants["source_data"] eval_data = constants["eval_data"] - Bplasma = constants["source_data"] + Bplasma = constants["Bplasma"] + if self._sheet_current: + sheet_eval_data = constants["sheet_eval_data"] else: source_data = compute_fun( "desc.equilibrium.equilibrium.Equilibrium", @@ -779,34 +816,34 @@ def compute(self, eq_params, *field_params, constants=None): # need extra factor of B/2 bc we're evaluating on plasma surface Bplasma = Bplasma + eval_data["B"] / 2 - if self._sheet_current: - p = "desc.magnetic_fields._current_potential.FourierCurrentPotentialField" - sheet_params = { - "R_lmn": eq_params["Rb_lmn"], - "Z_lmn": eq_params["Zb_lmn"], - "I": eq_params["I"], - "G": eq_params["G"], - "Phi_mn": eq_params["Phi_mn"], - } - sheet_source_data = compute_fun( - p, - self._sheet_data_keys, - params=sheet_params, - transforms=constants["sheet_source_transforms"], - profiles={}, - ) - sheet_eval_data = ( - sheet_source_data - if self._use_same_grid - else compute_fun( + if self._sheet_current: + p = "desc.magnetic_fields._current_potential.FourierCurrentPotentialField" + sheet_params = { + "R_lmn": eq_params["Rb_lmn"], + "Z_lmn": eq_params["Zb_lmn"], + "I": eq_params["I"], + "G": eq_params["G"], + "Phi_mn": eq_params["Phi_mn"], + } + sheet_source_data = compute_fun( p, self._sheet_data_keys, params=sheet_params, - transforms=constants["sheet_eval_transforms"], + transforms=constants["sheet_source_transforms"], profiles={}, ) - ) - source_data["K_vc"] += sheet_source_data["K"] + sheet_eval_data = ( + sheet_source_data + if self._use_same_grid + else compute_fun( + p, + self._sheet_data_keys, + params=sheet_params, + transforms=constants["sheet_eval_transforms"], + profiles={}, + ) + ) + source_data["K_vc"] += sheet_source_data["K"] x = jnp.array([eval_data["R"], eval_data["phi"], eval_data["Z"]]).T # can always pass in field params. If they're None, it just uses the From 0da0a74a697649e5ac43ecaebf2f001146b5fac7 Mon Sep 17 00:00:00 2001 From: Issra Ali Date: Fri, 10 Jul 2026 16:29:04 -0400 Subject: [PATCH 3/9] Added eq_fixed as a static attribute, updated changelog, made edits conformant with DESC formatting --- CHANGELOG.md | 1 + desc/objectives/_free_boundary.py | 28 ++++++++++++++++++---------- 2 files changed, 19 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d04b9a0491..6226411c4b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,6 +25,7 @@ New Features - Adds an optional attribute `ion_density` to the `Equilibrium` class, to allow the ion density profile to be set independently of the electron density and effective atomic number. Also adds compute functions for ``"ni_rr"`` and ``"Zeff_rr"``. - Modernizes dependencies to use [``nvidia-ml-py``](https://pypi.org/project/nvidia-ml-py/) in place of [``nvgpu``](https://github.com/rossumai/nvgpu). If you are updating an existing software environment uninstall ``pynvml`` first and then reinstall the dependencies to correctly get ``nvidia-ml-py``. +- Option for ``eq_fixed`` added to ``BoundaryError`` to remove equilibrium as a degree of freedom from boundary error objective. Bug Fixes diff --git a/desc/objectives/_free_boundary.py b/desc/objectives/_free_boundary.py index 23c1bdc2b8..89d8ed5a34 100644 --- a/desc/objectives/_free_boundary.py +++ b/desc/objectives/_free_boundary.py @@ -447,6 +447,7 @@ class BoundaryError(_Objective): "_B_plasma_chunk_size", "_bs_chunk_size", "_eq_data_keys", + "_eq_fixed", "_field_fixed", "_q", "_sheet_current", @@ -499,7 +500,7 @@ def __init__( field_fixed and eq_fixed, ValueError, "At least one of eq_fixed or field_fixed must be false.", - ) + ) self._field_fixed = field_fixed self._eq_fixed = eq_fixed self._bs_chunk_size = bs_chunk_size @@ -694,13 +695,16 @@ def build(self, use_jit=True, verbose=1): # need extra factor of B/2 bc we're evaluating on plasma surface Bplasma = Bplasma + eval_data["B"] / 2 - self._constants["source_data"]=source_data - self._constants["eval_data"]=eval_data - self._constants["Bplasma"]=Bplasma + self._constants["source_data"] = source_data + self._constants["eval_data"] = eval_data + self._constants["Bplasma"] = Bplasma - # sheet current stuff + # sheet current stuff if self._sheet_current: - p = "desc.magnetic_fields._current_potential.FourierCurrentPotentialField" + p = ( + "desc.magnetic_fields._current_potential." + "FourierCurrentPotentialField" + ) sheet_params = { "R_lmn": self._eq.params_dict["Rb_lmn"], "Z_lmn": self._eq.params_dict["Zb_lmn"], @@ -727,7 +731,7 @@ def build(self, use_jit=True, verbose=1): ) ) - self._constants["sheet_eval_data"]=sheet_eval_data + self._constants["sheet_eval_data"] = sheet_eval_data source_data["K_vc"] += sheet_source_data["K"] timer.stop("Precomputing transforms") @@ -759,7 +763,8 @@ def compute(self, eq_params, *field_params, constants=None): Parameters ---------- eq_params : dict - Dictionary of equilibrium degrees of freedom, eg Equilibrium.params_dict, if equilibrium is not fixed + Dictionary of equilibrium degrees of freedom, eg Equilibrium.params_dict, + if equilibrium is not fixed field_params : dict Dictionary of field parameters, if field is not fixed. constants : dict @@ -780,7 +785,7 @@ def compute(self, eq_params, *field_params, constants=None): eq_params = None constants = self._get_deprecated_constants(constants) - + if self._eq_fixed: source_data = constants["source_data"] eval_data = constants["eval_data"] @@ -817,7 +822,10 @@ def compute(self, eq_params, *field_params, constants=None): Bplasma = Bplasma + eval_data["B"] / 2 if self._sheet_current: - p = "desc.magnetic_fields._current_potential.FourierCurrentPotentialField" + p = ( + "desc.magnetic_fields._current_potential." + "FourierCurrentPotentialField" + ) sheet_params = { "R_lmn": eq_params["Rb_lmn"], "Z_lmn": eq_params["Zb_lmn"], From 0b0788a5a1aa3c96012829e55cfe41686c1dbb98 Mon Sep 17 00:00:00 2001 From: Issra Ali Date: Sat, 11 Jul 2026 15:49:26 -0400 Subject: [PATCH 4/9] Fixed signature of compute(), wrote a unit test --- desc/objectives/_free_boundary.py | 27 +++++++++++++++--------- tests/test_objective_funs.py | 34 +++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+), 10 deletions(-) diff --git a/desc/objectives/_free_boundary.py b/desc/objectives/_free_boundary.py index 89d8ed5a34..3efdeb556a 100644 --- a/desc/objectives/_free_boundary.py +++ b/desc/objectives/_free_boundary.py @@ -542,7 +542,7 @@ def build(self, use_jit=True, verbose=1): """ from desc.magnetic_fields import SumMagneticField - eq = self.things[0] + eq = self._eq if self._source_grid is None: # for axisymmetry we still need to know about toroidal effects, so its @@ -757,16 +757,19 @@ def build(self, use_jit=True, verbose=1): super().build(use_jit=use_jit, verbose=verbose) - def compute(self, eq_params, *field_params, constants=None): + def compute(self, params_1, params_2=None, constants=None): """Compute boundary force error. Parameters ---------- - eq_params : dict - Dictionary of equilibrium degrees of freedom, eg Equilibrium.params_dict, - if equilibrium is not fixed - field_params : dict - Dictionary of field parameters, if field is not fixed. + params_1 : dict + Dictionary of equilibrium or surface degrees of freedom, + eg ``Equilibrium.params_dict`` + Only required if ``self._eq_fixed = False``. + params_2 : dict + Dictionary of field degrees of freedom, eg ``Field.params_dict`` if + self._coils_fixed is False, else is the equilibrium or surface degrees of + freedom constants : dict Dictionary of constant data, eg transforms, profiles etc. Defaults to self.constants. (Deprecated) @@ -779,10 +782,14 @@ def compute(self, eq_params, *field_params, constants=None): √g||μ₀𝐊 − 𝐧 × [𝐁]|| in T*m^2 """ - if field_params == (): # common case for field_fixed=True + if self._eq_fixed: + field_params = params_1 + elif self._field_fixed: + eq_params = params_1 field_params = None - if eq_params == (): - eq_params = None + else: + eq_params = params_1 + field_params = params_2 constants = self._get_deprecated_constants(constants) diff --git a/tests/test_objective_funs.py b/tests/test_objective_funs.py index 5403aef6ba..19fba0c0ac 100644 --- a/tests/test_objective_funs.py +++ b/tests/test_objective_funs.py @@ -610,6 +610,40 @@ def test_boundary_error_biest(self): # next n should be B^2 errors np.testing.assert_allclose(f[n : 2 * n], 0, atol=5e-2) + @pytest.mark.unit + def test_boundary_error_things_fixed(self): + """Test that BoundaryError runs correctly for eq_fixed/field_fixed combos.""" + eq = Equilibrium(iota=PowerSeriesProfile(0)) + eq.solve() + field = ToroidalMagneticField(B0=1, R0=1) + + def test(eq_fixed=False, field_fixed=False): + obj = BoundaryError( + eq, + field, + eq_fixed=eq_fixed, + field_fixed=field_fixed, + ) + obj.build() + if eq_fixed: + f = obj.compute_unscaled(*obj.xs(field)) + elif field_fixed: + f = obj.compute_unscaled(*obj.xs(eq)) + else: + f = obj.compute_unscaled(*obj.xs(eq, field)) + n = len(f) // 2 + # first n should be B*n errors + np.testing.assert_allclose(f[:n], 0, atol=1e-4) + # next n should be B^2 errors + np.testing.assert_allclose(f[n : 2 * n], 0, atol=5e-2) + + test(eq_fixed=False, field_fixed=False) + test(eq_fixed=True) + test(field_fixed=True) + + with pytest.raises(ValueError, match="At least one"): + BoundaryError(eq, field, eq_fixed=True, field_fixed=True) + @pytest.mark.unit def test_boundary_error_vacuum(self): """Test calculation of vacuum boundary error.""" From aa36d9ecdb63844303f8327cab88ef46762b8885 Mon Sep 17 00:00:00 2001 From: Issra Ali Date: Sat, 11 Jul 2026 21:33:04 -0400 Subject: [PATCH 5/9] Fixed docstring of compute() --- desc/objectives/_free_boundary.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/desc/objectives/_free_boundary.py b/desc/objectives/_free_boundary.py index 3efdeb556a..ec0485148c 100644 --- a/desc/objectives/_free_boundary.py +++ b/desc/objectives/_free_boundary.py @@ -763,13 +763,13 @@ def compute(self, params_1, params_2=None, constants=None): Parameters ---------- params_1 : dict - Dictionary of equilibrium or surface degrees of freedom, - eg ``Equilibrium.params_dict`` - Only required if ``self._eq_fixed = False``. + Dictionary of equilibrium degrees of freedom, eg + ``Equilibrium.params_dict`` if self._eq_fixed is False, else is the field + degrees of freedom params_2 : dict - Dictionary of field degrees of freedom, eg ``Field.params_dict`` if - self._coils_fixed is False, else is the equilibrium or surface degrees of - freedom + Dictionary of field degrees of freedom, + eg ``CoilSet.params_dict`` + Only required if ``self._field_fixed = False``. constants : dict Dictionary of constant data, eg transforms, profiles etc. Defaults to self.constants. (Deprecated) From 87c46027db6653b5d48ff54d604a55113edc5992 Mon Sep 17 00:00:00 2001 From: Issra Ali Date: Tue, 14 Jul 2026 11:29:05 -0400 Subject: [PATCH 6/9] Fixed signature of compute() to accept any number of params, fixed unit test accordingly. --- desc/objectives/_free_boundary.py | 12 ++++++------ tests/test_objective_funs.py | 13 ++++++++----- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/desc/objectives/_free_boundary.py b/desc/objectives/_free_boundary.py index ec0485148c..5f218cb6f8 100644 --- a/desc/objectives/_free_boundary.py +++ b/desc/objectives/_free_boundary.py @@ -757,7 +757,7 @@ def build(self, use_jit=True, verbose=1): super().build(use_jit=use_jit, verbose=verbose) - def compute(self, params_1, params_2=None, constants=None): + def compute(self, *params, constants=None): """Compute boundary force error. Parameters @@ -783,14 +783,14 @@ def compute(self, params_1, params_2=None, constants=None): """ if self._eq_fixed: - field_params = params_1 + field_params = params + eq_params = None elif self._field_fixed: - eq_params = params_1 + eq_params = params[0] field_params = None else: - eq_params = params_1 - field_params = params_2 - + eq_params = params[0] + field_params = params[1:] constants = self._get_deprecated_constants(constants) if self._eq_fixed: diff --git a/tests/test_objective_funs.py b/tests/test_objective_funs.py index 19fba0c0ac..ebbeaac764 100644 --- a/tests/test_objective_funs.py +++ b/tests/test_objective_funs.py @@ -613,9 +613,11 @@ def test_boundary_error_biest(self): @pytest.mark.unit def test_boundary_error_things_fixed(self): """Test that BoundaryError runs correctly for eq_fixed/field_fixed combos.""" - eq = Equilibrium(iota=PowerSeriesProfile(0)) + eq = Equilibrium(L=3, M=3, N=3, Psi=np.pi) eq.solve() - field = ToroidalMagneticField(B0=1, R0=1) + coil = FourierXYZCoil(5e5) + coilset = CoilSet.linspaced_angular(coil, n=100, check_intersection=False) + field=[coilset, ToroidalMagneticField(B0=0, R0=1)] def test(eq_fixed=False, field_fixed=False): obj = BoundaryError( @@ -626,11 +628,12 @@ def test(eq_fixed=False, field_fixed=False): ) obj.build() if eq_fixed: - f = obj.compute_unscaled(*obj.xs(field)) + f = obj.compute(*[ff.params_dict for ff in field]) elif field_fixed: - f = obj.compute_unscaled(*obj.xs(eq)) + f = obj.compute(eq.params_dict) else: - f = obj.compute_unscaled(*obj.xs(eq, field)) + f = obj.compute(eq.params_dict, *[ff.params_dict for ff in field]) + f=obj.compute_scaled_error(*obj.xs()) n = len(f) // 2 # first n should be B*n errors np.testing.assert_allclose(f[:n], 0, atol=1e-4) From dc4cd3ed6529860b6d488c7136884dd7c67423bd Mon Sep 17 00:00:00 2001 From: Issra Ali Date: Tue, 14 Jul 2026 12:55:46 -0400 Subject: [PATCH 7/9] Fixed some errors in the test, added a test for sheet current case --- tests/test_objective_funs.py | 47 ++++++++++++++++++++++++++++++++---- 1 file changed, 42 insertions(+), 5 deletions(-) diff --git a/tests/test_objective_funs.py b/tests/test_objective_funs.py index ebbeaac764..dbe6593724 100644 --- a/tests/test_objective_funs.py +++ b/tests/test_objective_funs.py @@ -617,7 +617,7 @@ def test_boundary_error_things_fixed(self): eq.solve() coil = FourierXYZCoil(5e5) coilset = CoilSet.linspaced_angular(coil, n=100, check_intersection=False) - field=[coilset, ToroidalMagneticField(B0=0, R0=1)] + field = [coilset, ToroidalMagneticField(B0=0, R0=1)] def test(eq_fixed=False, field_fixed=False): obj = BoundaryError( @@ -628,12 +628,11 @@ def test(eq_fixed=False, field_fixed=False): ) obj.build() if eq_fixed: - f = obj.compute(*[ff.params_dict for ff in field]) + f = obj.compute_scaled_error(*[ff.params_dict for ff in field]) elif field_fixed: - f = obj.compute(eq.params_dict) + f = obj.compute_scaled_error(eq.params_dict) else: - f = obj.compute(eq.params_dict, *[ff.params_dict for ff in field]) - f=obj.compute_scaled_error(*obj.xs()) + f = obj.compute_scaled_error(eq.params_dict, *[ff.params_dict for ff in field]) n = len(f) // 2 # first n should be B*n errors np.testing.assert_allclose(f[:n], 0, atol=1e-4) @@ -647,6 +646,44 @@ def test(eq_fixed=False, field_fixed=False): with pytest.raises(ValueError, match="At least one"): BoundaryError(eq, field, eq_fixed=True, field_fixed=True) + @pytest.mark.unit + def test_boundary_error_things_fixed_sheet_current(self): + """Test that BoundaryError runs correctly for eq_fixed/field_fixed combos, + with sheet currents.""" + coil = FourierXYZCoil(5e5) + coilset = CoilSet.linspaced_angular(coil, n=100, check_intersection=False) + field = [coilset, ToroidalMagneticField(B0=0, R0=1)] + eq = Equilibrium(L=3, M=3, N=3, Psi=np.pi) + eq.surface = FourierCurrentPotentialField.from_surface( + eq.surface, M_Phi=eq.M, N_Phi=eq.N + ) + eq.solve() + + def test(eq_fixed=False, field_fixed=False): + obj = BoundaryError( + eq, + field, + eq_fixed=eq_fixed, + field_fixed=field_fixed, + ) + obj.build() + if eq_fixed: + f = obj.compute_scaled_error(*[ff.params_dict for ff in field]) + elif field_fixed: + f = obj.compute_scaled_error(eq.params_dict) + else: + f = obj.compute_scaled_error(eq.params_dict, *[ff.params_dict for ff in field]) + n = len(f) // 3 + # first n should be B*n errors + np.testing.assert_allclose(f[:n], 0, atol=1e-4) + # next n should be B^2 errors + np.testing.assert_allclose(f[n : 2 * n], 0, atol=5e-2) + # last n should be K errors + np.testing.assert_allclose(f[2 * n :], 0, atol=3e-2) + + with pytest.raises(ValueError, match="At least one"): + BoundaryError(eq, field, eq_fixed=True, field_fixed=True) + @pytest.mark.unit def test_boundary_error_vacuum(self): """Test calculation of vacuum boundary error.""" From 63c850ab921e0f014bfbb5624db12cb214c0021f Mon Sep 17 00:00:00 2001 From: Issra Ali Date: Tue, 14 Jul 2026 13:27:09 -0400 Subject: [PATCH 8/9] Made formatting changes --- tests/test_objective_funs.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/tests/test_objective_funs.py b/tests/test_objective_funs.py index dbe6593724..e99ab3d85d 100644 --- a/tests/test_objective_funs.py +++ b/tests/test_objective_funs.py @@ -632,7 +632,9 @@ def test(eq_fixed=False, field_fixed=False): elif field_fixed: f = obj.compute_scaled_error(eq.params_dict) else: - f = obj.compute_scaled_error(eq.params_dict, *[ff.params_dict for ff in field]) + f = obj.compute_scaled_error( + eq.params_dict, *[ff.params_dict for ff in field] + ) n = len(f) // 2 # first n should be B*n errors np.testing.assert_allclose(f[:n], 0, atol=1e-4) @@ -648,8 +650,7 @@ def test(eq_fixed=False, field_fixed=False): @pytest.mark.unit def test_boundary_error_things_fixed_sheet_current(self): - """Test that BoundaryError runs correctly for eq_fixed/field_fixed combos, - with sheet currents.""" + """Test BoundaryError for eq_fixed/field_fixed combos, with sheet currents.""" coil = FourierXYZCoil(5e5) coilset = CoilSet.linspaced_angular(coil, n=100, check_intersection=False) field = [coilset, ToroidalMagneticField(B0=0, R0=1)] @@ -672,7 +673,9 @@ def test(eq_fixed=False, field_fixed=False): elif field_fixed: f = obj.compute_scaled_error(eq.params_dict) else: - f = obj.compute_scaled_error(eq.params_dict, *[ff.params_dict for ff in field]) + f = obj.compute_scaled_error( + eq.params_dict, *[ff.params_dict for ff in field] + ) n = len(f) // 3 # first n should be B*n errors np.testing.assert_allclose(f[:n], 0, atol=1e-4) From b016feca52c1c61e64f97a204698e027f35c57e8 Mon Sep 17 00:00:00 2001 From: Issra Ali Date: Tue, 14 Jul 2026 14:55:18 -0400 Subject: [PATCH 9/9] Improved compute() docstring, cleaned up tests, fixed changelog --- CHANGELOG.md | 6 +++++- desc/objectives/_free_boundary.py | 17 ++++++----------- tests/test_objective_funs.py | 18 ++---------------- 3 files changed, 13 insertions(+), 28 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 905d489068..5c03516b4d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,11 @@ Changelog ========= +New Features + +- Added option for ``eq_fixed`` added to ``BoundaryError`` to remove equilibrium as a degree of freedom from boundary error objective. + + Performance Improvements - Speeds up the ``"qr"`` trust-region subproblem and Newton-step solves in the least-squares optimizers by reusing the Jacobian QR factorization across the Levenberg-Marquardt parameter sweep. On ``jax >= 0.10.0`` this uses ``qr_multiply`` to additionally avoid forming ``Q`` explicitly; on older versions a fallback preserves the same results. @@ -27,7 +32,6 @@ New Features - Adds an optional attribute `ion_density` to the `Equilibrium` class, to allow the ion density profile to be set independently of the electron density and effective atomic number. Also adds compute functions for ``"ni_rr"`` and ``"Zeff_rr"``. - Modernizes dependencies to use [``nvidia-ml-py``](https://pypi.org/project/nvidia-ml-py/) in place of [``nvgpu``](https://github.com/rossumai/nvgpu). If you are updating an existing software environment uninstall ``pynvml`` first and then reinstall the dependencies to correctly get ``nvidia-ml-py``. -- Option for ``eq_fixed`` added to ``BoundaryError`` to remove equilibrium as a degree of freedom from boundary error objective. Bug Fixes diff --git a/desc/objectives/_free_boundary.py b/desc/objectives/_free_boundary.py index 5f218cb6f8..e99d323c9c 100644 --- a/desc/objectives/_free_boundary.py +++ b/desc/objectives/_free_boundary.py @@ -640,9 +640,6 @@ def build(self, use_jit=True, verbose=1): neq = 3 if self._sheet_current else 2 # number of equations we're using - if self._eq_fixed: - eq = self._eq - self._constants = { "eval_transforms": eval_transforms, "eval_profiles": eval_profiles, @@ -762,14 +759,12 @@ def compute(self, *params, constants=None): Parameters ---------- - params_1 : dict - Dictionary of equilibrium degrees of freedom, eg - ``Equilibrium.params_dict`` if self._eq_fixed is False, else is the field - degrees of freedom - params_2 : dict - Dictionary of field degrees of freedom, - eg ``CoilSet.params_dict`` - Only required if ``self._field_fixed = False``. + params : dict + 1 or more dictionaries of params assigned in order of self.things, + respecting which degrees of freedom are fixed (via `eq_fixed` + and `field_fixed`). If eq_fixed, `params` are field_params, if + neither are fixed, `params[0]` takes equilibrium params, `params[1:]` + are field params)./ constants : dict Dictionary of constant data, eg transforms, profiles etc. Defaults to self.constants. (Deprecated) diff --git a/tests/test_objective_funs.py b/tests/test_objective_funs.py index e99ab3d85d..d4571bac9d 100644 --- a/tests/test_objective_funs.py +++ b/tests/test_objective_funs.py @@ -627,14 +627,7 @@ def test(eq_fixed=False, field_fixed=False): field_fixed=field_fixed, ) obj.build() - if eq_fixed: - f = obj.compute_scaled_error(*[ff.params_dict for ff in field]) - elif field_fixed: - f = obj.compute_scaled_error(eq.params_dict) - else: - f = obj.compute_scaled_error( - eq.params_dict, *[ff.params_dict for ff in field] - ) + f = obj.compute_scaled_error(*obj.xs()) n = len(f) // 2 # first n should be B*n errors np.testing.assert_allclose(f[:n], 0, atol=1e-4) @@ -668,14 +661,7 @@ def test(eq_fixed=False, field_fixed=False): field_fixed=field_fixed, ) obj.build() - if eq_fixed: - f = obj.compute_scaled_error(*[ff.params_dict for ff in field]) - elif field_fixed: - f = obj.compute_scaled_error(eq.params_dict) - else: - f = obj.compute_scaled_error( - eq.params_dict, *[ff.params_dict for ff in field] - ) + f = obj.compute_scaled_error(*obj.xs()) n = len(f) // 3 # first n should be B*n errors np.testing.assert_allclose(f[:n], 0, atol=1e-4)