Remove vaccum objective option from Equilibrium.solve

desc/continuation.py

@@ -464,7 +464,7 @@ def solve_continuation_automatic(  # noqa: C901
     ----------
     eq : Equilibrium
         Unsolved Equilibrium with the final desired boundary, profiles, resolution.
-    objective : {"force", "energy", "vacuum"}
+    objective : {"force", "energy"}
         function to solve for equilibrium solution
     optimizer : str or Optimizer (optional)
         optimizer to use
@@ -637,7 +637,7 @@ def solve_continuation(  # noqa: C901
     ----------
     eqfam : EquilibriaFamily or list of Equilibria
         Equilibria to solve for at each step.
-    objective : {"force", "energy", "vacuum"}
+    objective : {"force", "energy"}
         function to solve for equilibrium solution
     optimizer : str or Optimizer (optional)
         optimizer to use

desc/equilibrium/equilibrium.py

@@ -1809,7 +1809,7 @@ def solve(
 
         Parameters
         ----------
-        objective : {"force", "forces", "energy", "vacuum"}
+        objective : {"force", "forces", "energy"}
             Objective function to solve. Default = force balance on unified grid.
         constraints : Tuple
             set of constraints to enforce. Default = fixed boundary/profiles

desc/input_reader.py

@@ -189,6 +189,7 @@
         iota_flag = False
         pres_flag = False
         curr_flag = False
+        vac_flag = False
         inputs["output_path"] = self.output_path
 
         if self.args is not None and self.args.quiet:
@@ -360,7 +361,11 @@
             # solver methods
             match = re.search(r"objective", argument, re.IGNORECASE)
             if match:
-                inputs["objective"] = words[0].lower()
+                method = words[0].lower()
+                if method == "vacuum":
+                    method = "force"
+                    vac_flag = True
+                inputs["objective"] = method
                 flag = True
             match = re.search(r"optimizer", argument, re.IGNORECASE)
             if match:
@@ -558,20 +563,20 @@
         if curr_flag and iota_flag:
             raise OSError(colored("Cannot specify both iota and current.", "red"))
 
+        if vac_flag and (pres_flag or iota_flag or curr_flag):
+            warnings.warn(
+                "Vacuum objective assumes 0 pressure and 0 current, "
+                + "ignoring provided pressure, iota, and current profiles"
+            )
+            _ = inputs.pop("iota", None)
+            _ = inputs.pop("current", None)
+            _ = inputs.pop("pressure", None)
+
         # remove unused profile
         if iota_flag:
-            if inputs["objective"] != "vacuum":
-                del inputs["current"]
-            else:  # if vacuum objective from input file, use zero current
-                del inputs["iota"]
+            _ = inputs.pop("current", None)
         else:
-            del inputs["iota"]
-
-        if inputs["objective"] == "vacuum" and (pres_flag or iota_flag or curr_flag):
-            warnings.warn(
-                "Vacuum objective does not use any profiles, "
-                + "ignoring pressure, iota, and current"
-            )
+            _ = inputs.pop("iota", None)
 
         # sort axis array
         inputs["axis"] = inputs["axis"][inputs["axis"][:, 0].argsort()]
@@ -644,7 +649,7 @@
                 else:
                     inputs_ii[key] = inputs[key]
             # apply pressure ratio
-            if inputs_ii["pres_ratio"] is not None:
+            if "pressure" in inputs_ii and inputs_ii["pres_ratio"] is not None:
                 inputs_ii["pressure"][:, 1] *= inputs_ii["pres_ratio"]
             # apply current ratio
             if "current" in inputs_ii and inputs_ii["curr_ratio"] is not None:
@@ -824,6 +829,7 @@
         from desc.grid import LinearGrid
         from desc.io.equilibrium_io import load
         from desc.profiles import PowerSeriesProfile
+        from desc.utils import copy_coeffs
 
         f = open(outfile, "w+")
 
@@ -902,11 +908,15 @@
         pres_profile.change_resolution(L=L_profile)
         profile.change_resolution(L=L_profile)
 
+        prof_modes = np.zeros((L_profile, 3))
+        prof_modes[:, 0] = np.arange(L_profile)
+        p1 = copy_coeffs(pres_profile.params, pres_profile.basis.modes, prof_modes)
+        p2 = copy_coeffs(profile.params, profile.basis.modes, prof_modes)
         f.write("\n# pressure and rotational transform/current profiles\n")
         for l in range(L_profile):
             f.write(
                 "l: {:3d}  p = {:15.8E}  {} = {:15.8E}\n".format(
-                    int(l), pres_profile.params[l], char, profile.params[l]
+                    int(l), p1[l], char, p2[l]
                 )
             )
 

desc/objectives/getters.py

@@ -1,12 +1,6 @@
 """Utilities for getting standard groups of objectives and constraints."""
 
-from ._equilibrium import (
-    CurrentDensity,
-    Energy,
-    ForceBalance,
-    HelicalForceBalance,
-    RadialForceBalance,
-)
+from ._equilibrium import Energy, ForceBalance, HelicalForceBalance, RadialForceBalance
 from .linear_objectives import (
     AxisRSelfConsistency,
     AxisZSelfConsistency,
@@ -51,10 +45,10 @@ def get_equilibrium_objective(eq, mode="force", normalize=True):
     ----------
     eq : Equilibrium
         Equilibrium that will be optimized to satisfy the Objective.
-    mode : one of {"force", "forces", "energy", "vacuum"}
+    mode : one of {"force", "forces", "energy"}
         which objective to return. "force" computes force residuals on unified grid.
         "forces" uses two different grids for radial and helical forces. "energy" is
-        for minimizing MHD energy. "vacuum" directly minimizes current density.
+        for minimizing MHD energy.
     normalize : bool
         Whether to normalize units of objective.
 
@@ -74,10 +68,6 @@ def get_equilibrium_objective(eq, mode="force", normalize=True):
             RadialForceBalance(eq=eq, normalize=normalize, normalize_target=normalize),
             HelicalForceBalance(eq=eq, normalize=normalize, normalize_target=normalize),
         )
-    elif mode == "vacuum":
-        objectives = CurrentDensity(
-            eq=eq, normalize=normalize, normalize_target=normalize
-        )
     else:
         raise ValueError("got an unknown equilibrium objective type '{}'".format(mode))
     return ObjectiveFunction(objectives)

tests/conftest.py

@@ -212,7 +212,7 @@ def HELIOTRON_vac(tmpdir_factory):
     print("cwd=", cwd)
 
     args = ["-o", str(desc_h5_path), input_filename, "-vv"]
-    with pytest.warns(UserWarning, match="Vacuum objective does not use any profiles"):
+    with pytest.warns(UserWarning, match="Vacuum objective assumes 0 pressure"):
         main(args)
 
     HELIOTRON_vacuum_out = {

tests/test_input_output.py

@@ -134,13 +134,11 @@ def test_desc_output_to_input(tmpdir_factory):
     arr1 = arr1[arr1[:, 1].argsort()]
     arr1mneg = arr1[arr1[:, 1] < 0]
     arr1mpos = arr1[arr1[:, 1] >= 0]
-    pres1 = ir1.parse_inputs()[-1]["pressure"]
 
     desc_input_truth = "./tests/inputs/LandremanPaul2022_QA_reactorScale_lowRes"
     with pytest.warns(UserWarning):
         ir2 = InputReader(cl_args=[str(desc_input_truth)])
         arr2 = ir2.parse_inputs()[-1]["surface"]
-        pres2 = ir2.parse_inputs()[-1]["pressure"]
     arr2 = arr2[arr2[:, 1].argsort()]
     arr2mneg = arr2[arr2[:, 1] < 0]
     arr2mpos = arr2[arr2[:, 1] >= 0]
@@ -162,9 +160,6 @@ def test_desc_output_to_input(tmpdir_factory):
         atol=1e-8,
     )
 
-    if np.linalg.norm(pres1[:, 1]) > 0:
-        np.testing.assert_allclose(pres1(pres1[:, 1] > 0), pres2(pres2[:, 1] > 0))
-
     outfile_path = "./tests/inputs/iotest_HELIOTRON.h5"
     tmpdir = tmpdir_factory.mktemp("desc_inputs")
     tmp_path = tmpdir.join("input_iotest_HELIOTRON")
@@ -175,13 +170,16 @@ def test_desc_output_to_input(tmpdir_factory):
     ir1.desc_output_to_input(str(tmp_path), str(tmpout_path))
     ir1 = InputReader(cl_args=[str(tmp_path)])
     arr1 = ir1.parse_inputs()[-1]["surface"]
+    pres1 = ir1.parse_inputs()[-1]["pressure"]
+
     arr1 = arr1[arr1[:, 1].argsort()]
     arr1mneg = arr1[arr1[:, 1] < 0]
     arr1mpos = arr1[arr1[:, 1] >= 0]
 
     desc_input_truth = "./tests/inputs/iotest_HELIOTRON"
     ir2 = InputReader(cl_args=[str(desc_input_truth)])
     arr2 = ir2.parse_inputs()[-1]["surface"]
+    pres2 = ir2.parse_inputs()[-1]["pressure"]
     arr2 = arr2[arr2[:, 1].argsort()]
     arr2mneg = arr2[arr2[:, 1] < 0]
     arr2mpos = arr2[arr2[:, 1] >= 0]
@@ -204,8 +202,7 @@ def test_desc_output_to_input(tmpdir_factory):
         atol=1e-8,
     )
 
-    if np.linalg.norm(pres1[:, 1]) > 0:
-        np.testing.assert_allclose(pres1(pres1[:, 1] > 0), pres2(pres2[:, 1] > 0))
+    np.testing.assert_allclose(pres1[pres1[:, 1] > 0], pres2[pres2[:, 1] > 0])
 
 
 @pytest.mark.unit
@@ -315,9 +312,11 @@ def test_vacuum_objective_with_iota_yields_current(self):
         # load an input file with vacuum obj but also an iota profile specified
         with pytest.warns(UserWarning):
             ir = InputReader(input_path)
-        # ensure that a current profile instead of an iota profile is used
+        # ensure that no profiles used
         assert "iota" not in ir.inputs[0].keys()
-        assert "current" in ir.inputs[0].keys()
+        assert "current" not in ir.inputs[0].keys()
+        assert "pressure" not in ir.inputs[0].keys()
+        assert ir.inputs[0]["objective"] == "force"
 
     @pytest.mark.unit
     def test_node_pattern_warning(self):