diff --git a/desc/objectives/__init__.py b/desc/objectives/__init__.py index 7c6c95a68d..33ff6e22fd 100644 --- a/desc/objectives/__init__.py +++ b/desc/objectives/__init__.py @@ -1,7 +1,7 @@ """Classes defining objectives for equilibrium and optimization.""" from ._bootstrap import BootstrapRedlConsistency -from ._coils import CoilCurvature, CoilLength, CoilTorsion +from ._coils import CoilCurvature, CoilLength, CoilTorsion, QuadraticFlux from ._equilibrium import ( CurrentDensity, Energy, diff --git a/desc/objectives/_coils.py b/desc/objectives/_coils.py index 38cb32ad29..ec4f1daf9d 100644 --- a/desc/objectives/_coils.py +++ b/desc/objectives/_coils.py @@ -10,8 +10,10 @@ tree_structure, tree_unflatten, ) -from desc.compute import get_transforms +from desc.compute import compute as compute_fun +from desc.compute import get_profiles, get_transforms from desc.grid import LinearGrid, _Grid +from desc.singularities import compute_B_plasma from desc.utils import Timer, errorif from .normalization import compute_scaling_factors @@ -46,7 +48,7 @@ class _CoilObjective(_Objective): loss_function : {None, 'mean', 'min', 'max'}, optional Loss function to apply to the objective values once computed. This loss function is called on the raw compute value, before any shifting, scaling, or - normalization. Operates over all coils, not each individial coil. + normalization. Operates over all coils, not each individual coil. deriv_mode : {"auto", "fwd", "rev"} Specify how to compute jacobian matrix, either forward mode or reverse mode AD. "auto" selects forward or reverse mode based on the size of the input and output @@ -270,7 +272,7 @@ class CoilLength(_CoilObjective): loss_function : {None, 'mean', 'min', 'max'}, optional Loss function to apply to the objective values once computed. This loss function is called on the raw compute value, before any shifting, scaling, or - normalization. Operates over all coils, not each individial coil. + normalization. Operates over all coils, not each individual coil. deriv_mode : {"auto", "fwd", "rev"} Specify how to compute jacobian matrix, either forward mode or reverse mode AD. "auto" selects forward or reverse mode based on the size of the input and output @@ -400,7 +402,7 @@ class CoilCurvature(_CoilObjective): loss_function : {None, 'mean', 'min', 'max'}, optional Loss function to apply to the objective values once computed. This loss function is called on the raw compute value, before any shifting, scaling, or - normalization. Operates over all coils, not each individial coil. + normalization. Operates over all coils, not each individual coil. deriv_mode : {"auto", "fwd", "rev"} Specify how to compute jacobian matrix, either forward mode or reverse mode AD. "auto" selects forward or reverse mode based on the size of the input and output @@ -515,7 +517,7 @@ class CoilTorsion(_CoilObjective): loss_function : {None, 'mean', 'min', 'max'}, optional Loss function to apply to the objective values once computed. This loss function is called on the raw compute value, before any shifting, scaling, or - normalization. Operates over all coils, not each individial coil. + normalization. Operates over all coils, not each individual coil. deriv_mode : {"auto", "fwd", "rev"} Specify how to compute jacobian matrix, either forward mode or reverse mode AD. "auto" selects forward or reverse mode based on the size of the input and output @@ -597,3 +599,203 @@ def compute(self, params, constants=None): data = tree_flatten(data, is_leaf=lambda x: isinstance(x, dict))[0] out = jnp.concatenate([dat["torsion"] for dat in data]) return out + + +class QuadraticFlux(_Objective): + """Target B*n = 0 on LCFS. + + Uses virtual casing to find plasma component of B and penalizes + (B_coil + B_plasma)*n. The equilibrium is kept fixed while the + field is unfixed. + + Parameters + ---------- + eq : Equilibrium + Equilibrium upon whose surface the normal field error will be minimized. + The equilibrium is kept fixed during the optimization with this objective. + field : MagneticField + External field produced by coils or other source, which will be optimized to + minimize the normal field error on the provided equilibrium's surface. + target : float, ndarray, optional + Target value(s) of the objective. Only used if bounds is None. + len(target) must be equal to Objective.dim_f. + Default target is zero. + bounds : tuple, optional + Lower and upper bounds on the objective. Overrides target. + len(bounds[0]) and len(bounds[1]) must be equal to Objective.dim_f + weight : float, ndarray, optional + Weighting to apply to the Objective, relative to other Objectives. + len(weight) must be equal to Objective.dim_f + normalize : bool + Whether to compute the error in physical units or non-dimensionalize. + normalize_target : bool + Whether target and bounds should be normalized before comparing to computed + values. If `normalize` is `True` and the target is in physical units, + this should also be set to True. + source_grid : Grid, optional + Collocation grid containing the nodes for plasma source terms. + Default grid is detailed in the docs for ``compute_B_plasma`` + eval_grid : Grid, optional + Collocation grid containing the nodes on the plasma surface at which the + magnetic field is being calculated and where to evaluate Bn errors. + Default grid is: LinearGrid(rho=np.array([1.0]), M=eq.M_grid, N=eq.N_grid, + NFP=int(eq.NFP), sym=False) + field_grid : Grid, optional + Grid used to discretize field (e.g. grid for the magnetic field source from + coils). Default grid is determined by the specific MagneticField object, see + the docs of that object's ``compute_magnetic_field`` method for more detail. + vacuum : bool + If true, B_plasma (the contribution to the normal field on the boundary from the + plasma currents) is set to zero. + name : str + Name of the objective function. + + """ + + _scalar = False + _linear = False + _print_value_fmt = "Boundary normal field error: {:10.3e} " + _units = "(T m^2)" + _coordinates = "rtz" + + def __init__( + self, + eq, + field, + target=None, + bounds=None, + weight=1, + normalize=True, + normalize_target=True, + source_grid=None, + eval_grid=None, + field_grid=None, + vacuum=False, + name="Quadratic flux", + ): + if target is None and bounds is None: + target = 0 + self._source_grid = source_grid + self._eval_grid = eval_grid + self._eq = eq + self._field = field + self._field_grid = field_grid + self._vacuum = vacuum + things = [field] + super().__init__( + things=things, + target=target, + bounds=bounds, + weight=weight, + normalize=normalize, + normalize_target=normalize_target, + name=name, + ) + + def build(self, use_jit=True, verbose=1): + """Build constant arrays. + + Parameters + ---------- + use_jit : bool, optional + Whether to just-in-time compile the objective and derivatives. + verbose : int, optional + Level of output. + + """ + eq = self._eq + + if self._eval_grid is None: + eval_grid = LinearGrid( + rho=np.array([1.0]), + M=eq.M_grid, + N=eq.N_grid, + NFP=int(eq.NFP), + sym=False, + ) + self._eval_grid = eval_grid + else: + eval_grid = self._eval_grid + + self._data_keys = ["R", "Z", "n_rho", "phi", "|e_theta x e_zeta|"] + + timer = Timer() + if verbose > 0: + print("Precomputing transforms") + timer.start("Precomputing transforms") + + self._dim_f = eval_grid.num_nodes + + w = eval_grid.weights + w *= jnp.sqrt(eval_grid.num_nodes) + + eval_profiles = get_profiles(self._data_keys, obj=eq, grid=eval_grid) + eval_transforms = get_transforms(self._data_keys, obj=eq, grid=eval_grid) + eval_data = compute_fun( + "desc.equilibrium.equilibrium.Equilibrium", + self._data_keys, + params=eq.params_dict, + transforms=eval_transforms, + profiles=eval_profiles, + ) + + # pre-compute B_plasma because we are assuming eq is fixed + if self._vacuum: + Bplasma = jnp.zeros(eval_grid.num_nodes) + + else: + Bplasma = compute_B_plasma( + eq, eval_grid, self._source_grid, normal_only=True + ) + + self._constants = { + "field": self._field, + "field_grid": self._field_grid, + "quad_weights": w, + "eval_data": eval_data, + "B_plasma": Bplasma, + } + + timer.stop("Precomputing transforms") + if verbose > 1: + timer.disp("Precomputing transforms") + + if self._normalize: + scales = compute_scaling_factors(eq) + self._normalization = scales["B"] * scales["R0"] * scales["a"] + + super().build(use_jit=use_jit, verbose=verbose) + + def compute(self, field_params, constants=None): + """Compute boundary force error. + + Parameters + ---------- + field_params : dict + Dictionary of the external field's degrees of freedom. + constants : dict + Dictionary of constant data, eg transforms, profiles etc. Defaults to + self.constants + + Returns + ------- + f : ndarray + Bnorm from B_ext and B_plasma + + """ + if constants is None: + constants = self.constants + + # B_plasma from equilibrium precomputed + eval_data = constants["eval_data"] + B_plasma = constants["B_plasma"] + + x = jnp.array([eval_data["R"], eval_data["phi"], eval_data["Z"]]).T + + # B_ext is not pre-computed because field is not fixed + B_ext = constants["field"].compute_magnetic_field( + x, source_grid=constants["field_grid"], basis="rpz", params=field_params + ) + B_ext = jnp.sum(B_ext * eval_data["n_rho"], axis=-1) + f = (B_ext + B_plasma) * eval_data["|e_theta x e_zeta|"] + return f diff --git a/desc/objectives/getters.py b/desc/objectives/getters.py index fc157d07a8..d9151c6cef 100644 --- a/desc/objectives/getters.py +++ b/desc/objectives/getters.py @@ -2,6 +2,8 @@ import numpy as np +from desc.utils import is_any_instance + from ._equilibrium import Energy, ForceBalance, HelicalForceBalance, RadialForceBalance from .linear_objectives import ( AxisRSelfConsistency, @@ -220,10 +222,6 @@ def get_NAE_constraints( def maybe_add_self_consistency(thing, constraints): """Add self consistency constraints if needed.""" - - def _is_any_instance(things, cls): - return any([isinstance(t, cls) for t in things]) - # Equilibrium if ( hasattr(thing, "Ra_n") @@ -232,22 +230,22 @@ def _is_any_instance(things, cls): and hasattr(thing, "Zb_lmn") and hasattr(thing, "L_lmn") ): - if not _is_any_instance(constraints, BoundaryRSelfConsistency): + if not is_any_instance(constraints, BoundaryRSelfConsistency): constraints += (BoundaryRSelfConsistency(eq=thing),) - if not _is_any_instance(constraints, BoundaryZSelfConsistency): + if not is_any_instance(constraints, BoundaryZSelfConsistency): constraints += (BoundaryZSelfConsistency(eq=thing),) - if not _is_any_instance(constraints, FixLambdaGauge): + if not is_any_instance(constraints, FixLambdaGauge): constraints += (FixLambdaGauge(eq=thing),) - if not _is_any_instance(constraints, AxisRSelfConsistency): + if not is_any_instance(constraints, AxisRSelfConsistency): constraints += (AxisRSelfConsistency(eq=thing),) - if not _is_any_instance(constraints, AxisZSelfConsistency): + if not is_any_instance(constraints, AxisZSelfConsistency): constraints += (AxisZSelfConsistency(eq=thing),) # Curve elif hasattr(thing, "shift") and hasattr(thing, "rotmat"): - if not _is_any_instance(constraints, FixCurveShift): + if not is_any_instance(constraints, FixCurveShift): constraints += (FixCurveShift(curve=thing),) - if not _is_any_instance(constraints, FixCurveRotation): + if not is_any_instance(constraints, FixCurveRotation): constraints += (FixCurveRotation(curve=thing),) return constraints diff --git a/desc/optimize/optimizer.py b/desc/optimize/optimizer.py index 63c57bf7f7..08a4a8a696 100644 --- a/desc/optimize/optimizer.py +++ b/desc/optimize/optimizer.py @@ -13,7 +13,15 @@ maybe_add_self_consistency, ) from desc.objectives.utils import combine_args -from desc.utils import Timer, flatten_list, get_instance, unique_list, warnif +from desc.utils import ( + Timer, + errorif, + flatten_list, + get_instance, + is_any_instance, + unique_list, + warnif, +) from ._constraint_wrappers import LinearConstraintProjection, ProximalProjection @@ -147,6 +155,11 @@ def optimize( # noqa: C901 - FIXME: simplify this """ if not isinstance(constraints, (tuple, list)): constraints = (constraints,) + errorif( + not isinstance(objective, ObjectiveFunction), + TypeError, + "objective should be of type ObjectiveFunction.", + ) # get unique things things, indices = unique_list(flatten_list(things, flatten_tuple=True)) @@ -161,10 +174,19 @@ def optimize( # noqa: C901 - FIXME: simplify this # need local import to avoid circular dependencies from desc.equilibrium import Equilibrium + from desc.objectives import QuadraticFlux # eq may be None eq = get_instance(things, Equilibrium) if eq is not None: + # check if stage 2 objectives are here: + errorif( + is_any_instance(constraints, QuadraticFlux) + or is_any_instance(objective.objectives, QuadraticFlux), + ValueError, + "QuadraticFlux objective assumes Equilibrium is fixed but Equilibrium " + + "is in things to optimize.", + ) # save these for later eq_params_init = eq.params_dict.copy() diff --git a/desc/singularities.py b/desc/singularities.py index 1289792daf..2c6f0f0bbe 100644 --- a/desc/singularities.py +++ b/desc/singularities.py @@ -834,16 +834,7 @@ def compute_B_plasma(eq, eval_grid, source_grid=None, normal_only=False): sym=False, ) - data_keys = [ - "K_vc", - "B", - "R", - "phi", - "Z", - "e^rho", - "n_rho", - "|e_theta x e_zeta|", - ] + data_keys = ["K_vc", "B", "R", "phi", "Z", "e^rho", "n_rho", "|e_theta x e_zeta|"] eval_data = eq.compute(data_keys, grid=eval_grid) source_data = eq.compute(data_keys, grid=source_grid) @@ -860,13 +851,9 @@ def compute_B_plasma(eq, eval_grid, source_grid=None, normal_only=False): interpolator = DFTInterpolator(eval_grid, source_grid, s, q) if hasattr(eq.surface, "Phi_mn"): source_data["K_vc"] += eq.surface.compute("K", grid=source_grid)["K"] - Bplasma = virtual_casing_biot_savart( - eval_data, - source_data, - interpolator, - ) + Bplasma = virtual_casing_biot_savart(eval_data, source_data, interpolator) # need extra factor of B/2 bc we're evaluating on plasma surface Bplasma = Bplasma + eval_data["B"] / 2 if normal_only: - Bplasma = jnp.sum(Bplasma * source_data["n_rho"], axis=1) + Bplasma = jnp.sum(Bplasma * eval_data["n_rho"], axis=1) return Bplasma diff --git a/desc/utils.py b/desc/utils.py index 57a58b9c20..710532d432 100644 --- a/desc/utils.py +++ b/desc/utils.py @@ -603,3 +603,8 @@ def unique_list(thelist): unique.append(x) inds.append(unique.index(x)) return unique, inds + + +def is_any_instance(things, cls): + """Check if any of things is an instance of cls.""" + return any([isinstance(t, cls) for t in things]) diff --git a/tests/inputs/vacuum_circular_tokamak.h5 b/tests/inputs/vacuum_circular_tokamak.h5 new file mode 100644 index 0000000000..0b592fce08 Binary files /dev/null and b/tests/inputs/vacuum_circular_tokamak.h5 differ diff --git a/tests/inputs/vacuum_nonaxisym.h5 b/tests/inputs/vacuum_nonaxisym.h5 new file mode 100644 index 0000000000..be1d81dca4 Binary files /dev/null and b/tests/inputs/vacuum_nonaxisym.h5 differ diff --git a/tests/test_examples.py b/tests/test_examples.py index 3f93c58921..aa6e8c2503 100644 --- a/tests/test_examples.py +++ b/tests/test_examples.py @@ -17,7 +17,11 @@ from desc.geometry import FourierRZToroidalSurface from desc.grid import LinearGrid from desc.io import load -from desc.magnetic_fields import OmnigenousField, SplineMagneticField +from desc.magnetic_fields import ( + OmnigenousField, + SplineMagneticField, + ToroidalMagneticField, +) from desc.objectives import ( AspectRatio, BoundaryError, @@ -44,6 +48,7 @@ Omnigenity, PlasmaVesselDistance, PrincipalCurvature, + QuadraticFlux, QuasisymmetryBoozer, QuasisymmetryTwoTerm, VacuumBoundaryError, @@ -1230,3 +1235,44 @@ def test_single_coil_optimization(): np.testing.assert_allclose( coil.compute("torsion", grid=grid)["torsion"], target, atol=1e-5 ) + + +@pytest.mark.unit +def test_quadratic_flux_optimization_with_analytic_field(): + """Test analytic field optimization to reduce quadratic flux. + + Checks that B goes to zero for non-axisymmetric eq and axisymmetric field. + """ + eq = get("precise_QA") + field = ToroidalMagneticField(1, 1) + eval_grid = LinearGrid( + rho=np.array([1.0]), + M=eq.M_grid, + N=eq.N_grid, + NFP=eq.NFP, + sym=False, + ) + + optimizer = Optimizer("lsq-exact") + + constraints = (FixParameter(field, ["R0"]),) + quadflux_obj = QuadraticFlux( + eq=eq, + field=field, + eval_grid=eval_grid, + vacuum=True, + ) + objective = ObjectiveFunction(quadflux_obj) + things, __ = optimizer.optimize( + field, + objective=objective, + constraints=constraints, + ftol=1e-14, + gtol=1e-14, + copy=True, + verbose=3, + ) + + # optimizer should zero out field since that's the easiest way + # to get to Bnorm = 0 + np.testing.assert_allclose(things[0].B0, 0, atol=1e-12) diff --git a/tests/test_objective_funs.py b/tests/test_objective_funs.py index d8cf475014..ac0a0a3512 100644 --- a/tests/test_objective_funs.py +++ b/tests/test_objective_funs.py @@ -20,10 +20,12 @@ from desc.examples import get from desc.geometry import FourierRZToroidalSurface from desc.grid import ConcentricGrid, LinearGrid, QuadratureGrid +from desc.io import load from desc.magnetic_fields import ( FourierCurrentPotentialField, OmnigenousField, SplineMagneticField, + ToroidalMagneticField, ) from desc.objectives import ( AspectRatio, @@ -49,6 +51,7 @@ PlasmaVesselDistance, Pressure, PrincipalCurvature, + QuadraticFlux, QuasisymmetryBoozer, QuasisymmetryTripleProduct, QuasisymmetryTwoTerm, @@ -816,6 +819,64 @@ def test(coil, grid=None): test(mixed_coils, grid=[LinearGrid(N=5)] * len(mixed_coils.coils)) test(nested_coils, grid=nested_grids) + @pytest.mark.unit + def test_quadratic_flux(self): + """Test calculation of quadratic flux on the boundary.""" + t_field = ToroidalMagneticField(1, 1) + + # test that torus (axisymmetric) Bnorm is exactly 0 + eq = load("./tests/inputs/vacuum_circular_tokamak.h5") + obj = QuadraticFlux(eq, t_field) + obj.build(eq, verbose=2) + f = obj.compute(field_params=t_field.params_dict) + np.testing.assert_allclose(f, 0, rtol=1e-14, atol=1e-14) + + # test non-axisymmetric surface + eq = desc.examples.get("precise_QA", "all")[0] + eq.change_resolution(4, 4, 4, 8, 8, 8) + eval_grid = LinearGrid( + rho=np.array([1.0]), + M=eq.M_grid, + N=eq.N_grid, + NFP=eq.NFP, + sym=False, + ) + source_grid = LinearGrid( + rho=np.array([1.0]), + M=eq.M_grid, + N=eq.N_grid, + NFP=eq.NFP, + sym=False, + ) + + obj = QuadraticFlux(eq, t_field, eval_grid=eval_grid, source_grid=source_grid) + Bnorm = t_field.compute_Bnormal( + eq.surface, eval_grid=eval_grid, source_grid=source_grid + )[0] + obj.build(eq) + dA = eq.compute("|e_theta x e_zeta|", grid=eval_grid)["|e_theta x e_zeta|"] + f = obj.compute(field_params=t_field.params_dict) + + np.testing.assert_allclose(f, Bnorm * dA, atol=2e-4, rtol=1e-2) + + # equilibrium that has B_plasma == 0 + eq = load("./tests/inputs/vacuum_nonaxisym.h5") + + eval_grid = LinearGrid( + rho=np.array([1.0]), + M=eq.M_grid, + N=eq.N_grid, + NFP=eq.NFP, + sym=False, + ) + obj = QuadraticFlux(eq, t_field, vacuum=True, eval_grid=eval_grid) + Bnorm = t_field.compute_Bnormal(eq.surface, eval_grid=eval_grid)[0] + obj.build(eq) + f = obj.compute(field_params=t_field.params_dict) + dA = eq.compute("|e_theta x e_zeta|", grid=eval_grid)["|e_theta x e_zeta|"] + # check that they're the same since we set B_plasma = 0 + np.testing.assert_allclose(f, Bnorm * dA, atol=1e-14) + @pytest.mark.regression def test_derivative_modes(): @@ -1626,6 +1687,7 @@ class TestComputeScalarResolution: CoilLength, CoilTorsion, CoilCurvature, + QuadraticFlux, # need to avoid blowup near the axis MercierStability, # don't test these since they depend on what user wants @@ -1700,7 +1762,7 @@ def test_compute_scalar_resolution_boundary_error(self): eq.change_resolution( L_grid=int(eq.L * res), M_grid=int(eq.M * res), N_grid=int(eq.N * res) ) - obj = ObjectiveFunction(BoundaryError(self.eq, ext_field), use_jit=False) + obj = ObjectiveFunction(BoundaryError(eq, ext_field), use_jit=False) obj.build(verbose=0) f[i] = obj.compute_scalar(obj.x()) np.testing.assert_allclose(f, f[-1], rtol=5e-2) @@ -1726,14 +1788,38 @@ def test_compute_scalar_resolution_vacuum_boundary_error(self): eq.change_resolution( L_grid=int(eq.L * res), M_grid=int(eq.M * res), N_grid=int(eq.N * res) ) - obj = ObjectiveFunction( - VacuumBoundaryError(self.eq, ext_field), use_jit=False - ) + obj = ObjectiveFunction(VacuumBoundaryError(eq, ext_field), use_jit=False) with pytest.warns(UserWarning): obj.build(verbose=0) f[i] = obj.compute_scalar(obj.x()) np.testing.assert_allclose(f, f[-1], rtol=5e-2) + @pytest.mark.regression + def test_compute_scalar_resolution_quadratic_flux(self): + """VacuumBoundaryError.""" + ext_field = SplineMagneticField.from_mgrid(r"tests/inputs/mgrid_solovev.nc") + + pres = PowerSeriesProfile([1.25e-1, 0, -1.25e-1]) + iota = PowerSeriesProfile([-4.9e-1, 0, 3.0e-1]) + surf = FourierRZToroidalSurface( + R_lmn=[4.0, 1.0], + modes_R=[[0, 0], [1, 0]], + Z_lmn=[-1.0], + modes_Z=[[-1, 0]], + NFP=1, + ) + eq = Equilibrium(M=6, N=0, Psi=1.0, surface=surf, pressure=pres, iota=iota) + + f = np.zeros_like(self.res_array, dtype=float) + for i, res in enumerate(self.res_array): + eq.change_resolution( + L_grid=int(eq.L * res), M_grid=int(eq.M * res), N_grid=int(eq.N * res) + ) + obj = ObjectiveFunction(QuadraticFlux(eq, ext_field), use_jit=False) + obj.build(verbose=0) + f[i] = obj.compute_scalar(obj.x()) + np.testing.assert_allclose(f, f[-1], rtol=5e-2) + @pytest.mark.regression def test_compute_scalar_resolution_generic_scalar(self): """Generic objective with scalar qty.""" @@ -1903,6 +1989,7 @@ class TestObjectiveNaNGrad: CoilLength, CoilCurvature, CoilTorsion, + QuadraticFlux, # we don't test these since they depend too much on what exactly the user wants GenericObjective, LinearObjectiveFromUser, @@ -1991,6 +2078,28 @@ def test_objective_no_nangrad_vacuum_boundary_error(self): g = obj.grad(obj.x(eq, ext_field)) assert not np.any(np.isnan(g)), "vacuum boundary error" + @pytest.mark.unit + def test_objective_no_nangrad_quadratic_flux(self): + """QuadraticFlux.""" + ext_field = SplineMagneticField.from_mgrid(r"tests/inputs/mgrid_solovev.nc") + + pres = PowerSeriesProfile([1.25e-1, 0, -1.25e-1]) + iota = PowerSeriesProfile([-4.9e-1, 0, 3.0e-1]) + surf = FourierRZToroidalSurface( + R_lmn=[4.0, 1.0], + modes_R=[[0, 0], [1, 0]], + Z_lmn=[-1.0], + modes_Z=[[-1, 0]], + NFP=1, + ) + + eq = Equilibrium(M=6, N=0, Psi=1.0, surface=surf, pressure=pres, iota=iota) + + obj = ObjectiveFunction(QuadraticFlux(eq, ext_field), use_jit=False) + obj.build() + g = obj.grad(obj.x(ext_field)) + assert not np.any(np.isnan(g)), "quadratic flux" + @pytest.mark.unit @pytest.mark.parametrize( "objective", sorted(other_objectives, key=lambda x: str(x.__name__))