Let the vacuum potential's Fourier cutoffs exceed the plasma's

src/vmecpp/__init__.py

@@ -260,6 +260,20 @@ class VmecInput(BaseModelWithNumpy):
     the full ntor. < 0 (default) means geometry uses ntor.
     """
 
+    vacuum_mpol: int = 0
+    """Poloidal Fourier cutoff of the vacuum potential in a free-boundary run.
+
+    NESTOR expands the vacuum potential to mpol like the plasma; a value above mpol
+    raises the potential's cutoff alone. 0 (default) means the potential uses mpol.
+    """
+
+    vacuum_ntor: int = 0
+    """Toroidal Fourier cutoff of the vacuum potential in a free-boundary run.
+
+    A value above ntor raises the potential's cutoff alone; nzeta must then be at least
+    2 * vacuum_ntor + 4. 0 (default) means the potential uses ntor.
+    """
+
     ntheta: int = 0
     """Number of poloidal grid points (ntheta >= 0).
 

src/vmecpp/cpp/vmecpp/common/sizes/sizes.cc

@@ -13,9 +13,31 @@
 
 namespace vmecpp {
 
+namespace {
+// The vacuum potential's cutoff is the plasma's unless raised by the input.
+int VacuumCutoff(int vacuum, int plasma) {
+  return vacuum > plasma ? vacuum : plasma;
+}
+
+// The tangential grid has to resolve the larger of the plasma's and the
+// vacuum potential's cutoffs; these are the Nyquist minima computeDerivedSizes
+// applies for the plasma's, evaluated for the larger one.
+int NthetaForVacuum(const VmecINDATA& id) {
+  return std::max(id.ntheta, 2 * VacuumCutoff(id.vacuum_mpol, id.mpol) + 6);
+}
+
+int NzetaForVacuum(const VmecINDATA& id) {
+  const int ntor = VacuumCutoff(id.vacuum_ntor, id.ntor);
+  if (ntor > 0 && id.nzeta < 2 * ntor + 4) {
+    return 2 * ntor + 4;
+  }
+  return id.nzeta;
+}
+}  // namespace
+
 Sizes::Sizes(const VmecINDATA& id)
-    : Sizes(id.lasym, id.nfp, id.mpol, id.ntor, id.ntheta, id.nzeta,
-            id.mpol_geometry, id.ntor_geometry) {}
+    : Sizes(id.lasym, id.nfp, id.mpol, id.ntor, NthetaForVacuum(id),
+            NzetaForVacuum(id), id.mpol_geometry, id.ntor_geometry) {}
 
 Sizes::Sizes(bool lasym, int nfp, int mpol, int ntor, int ntheta, int nzeta,
              int mpol_geometry, int ntor_geometry)

src/vmecpp/cpp/vmecpp/common/vmec_indata/vmec_indata.cc

@@ -174,6 +174,8 @@ VmecINDATA::VmecINDATA() {
   mpol = 6;
   ntor = 0;
   mpol_geometry = -1;
+  vacuum_mpol = 0;
+  vacuum_ntor = 0;
   ntor_geometry = -1;
   ntheta = 0;
   nzeta = 0;
@@ -326,6 +328,8 @@ absl::Status VmecINDATA::WriteTo(H5::H5File& file) const {
   WriteH5Dataset(mpol, "/indata/mpol", file);
   WriteH5Dataset(ntor, "/indata/ntor", file);
   WriteH5Dataset(mpol_geometry, "/indata/mpol_geometry", file);
+  WriteH5Dataset(vacuum_mpol, "/indata/vacuum_mpol", file);
+  WriteH5Dataset(vacuum_ntor, "/indata/vacuum_ntor", file);
   WriteH5Dataset(ntor_geometry, "/indata/ntor_geometry", file);
   WriteH5Dataset(ntheta, "/indata/ntheta", file);
   WriteH5Dataset(nzeta, "/indata/nzeta", file);
@@ -408,6 +412,12 @@ absl::Status VmecINDATA::LoadInto(VmecINDATA& m_indata, H5::H5File& from_file) {
   if (from_file.nameExists("/indata/ntor_geometry")) {
     ReadH5Dataset(m_indata.ntor_geometry, "/indata/ntor_geometry", from_file);
   }
+  if (from_file.nameExists("/indata/vacuum_mpol")) {
+    ReadH5Dataset(m_indata.vacuum_mpol, "/indata/vacuum_mpol", from_file);
+  }
+  if (from_file.nameExists("/indata/vacuum_ntor")) {
+    ReadH5Dataset(m_indata.vacuum_ntor, "/indata/vacuum_ntor", from_file);
+  }
   ReadH5Dataset(m_indata.ntheta, "/indata/ntheta", from_file);
   ReadH5Dataset(m_indata.nzeta, "/indata/nzeta", from_file);
   ReadH5Dataset(m_indata.phiedge, "/indata/phiedge", from_file);
@@ -608,6 +618,22 @@ absl::StatusOr<VmecINDATA> VmecINDATA::FromJson(
     vmec_indata.ntor_geometry = maybe_ntor_geometry->value();
   }
 
+  auto maybe_vacuum_mpol = JsonReadInt(j, "vacuum_mpol");
+  if (!maybe_vacuum_mpol.ok()) {
+    return maybe_vacuum_mpol.status();
+  }
+  if (maybe_vacuum_mpol->has_value()) {
+    vmec_indata.vacuum_mpol = maybe_vacuum_mpol->value();
+  }
+
+  auto maybe_vacuum_ntor = JsonReadInt(j, "vacuum_ntor");
+  if (!maybe_vacuum_ntor.ok()) {
+    return maybe_vacuum_ntor.status();
+  }
+  if (maybe_vacuum_ntor->has_value()) {
+    vmec_indata.vacuum_ntor = maybe_vacuum_ntor->value();
+  }
+
   auto maybe_ntheta = JsonReadInt(j, "ntheta");
   if (!maybe_ntheta.ok()) {
     return maybe_ntheta.status();
@@ -1205,6 +1231,8 @@ absl::StatusOr<std::string> VmecINDATA::ToJson() const {
   output["ntor"] = ntor;
   output["mpol_geometry"] = mpol_geometry;
   output["ntor_geometry"] = ntor_geometry;
+  output["vacuum_mpol"] = vacuum_mpol;
+  output["vacuum_ntor"] = vacuum_ntor;
   output["ntheta"] = ntheta;
   output["nzeta"] = nzeta;
 
@@ -1350,6 +1378,32 @@ absl::Status IsConsistent(const VmecINDATA& vmec_indata,
 
   /* --------------------------------- */
 
+  // vacuum_mpol, vacuum_ntor
+  // * 0 means the vacuum potential uses mpol / ntor; otherwise the cutoff may
+  //   only exceed the plasma's, since the boundary has to fit into it
+  // * the toroidal grid has to resolve the potential's toroidal cutoff
+  if (vmec_indata.vacuum_mpol != 0 &&
+      vmec_indata.vacuum_mpol < vmec_indata.mpol) {
+    return absl::InvalidArgumentError(absl::StrFormat(
+        "input variable 'vacuum_mpol' must be 0 or at least mpol = %d, but "
+        "is %d\n",
+        vmec_indata.mpol, vmec_indata.vacuum_mpol));
+  }
+  if (vmec_indata.vacuum_ntor != 0 &&
+      vmec_indata.vacuum_ntor < vmec_indata.ntor) {
+    return absl::InvalidArgumentError(absl::StrFormat(
+        "input variable 'vacuum_ntor' must be 0 or at least ntor = %d, but "
+        "is %d\n",
+        vmec_indata.ntor, vmec_indata.vacuum_ntor));
+  }
+  if (vmec_indata.vacuum_ntor > vmec_indata.ntor && vmec_indata.nzeta > 0 &&
+      vmec_indata.nzeta < 2 * vmec_indata.vacuum_ntor + 4) {
+    return absl::InvalidArgumentError(absl::StrFormat(
+        "input variable 'nzeta' must be at least 2 * vacuum_ntor + 4 = %d to "
+        "carry the vacuum potential's toroidal cutoff, but is %d\n",
+        2 * vmec_indata.vacuum_ntor + 4, vmec_indata.nzeta));
+  }
+
   const int NS_MIN = 3;
 
   // ns_array

src/vmecpp/cpp/vmecpp/common/vmec_indata/vmec_indata.h

@@ -86,6 +86,15 @@ class VmecINDATA {
   int mpol_geometry;
   int ntor_geometry;
 
+  // Optional larger Fourier cutoffs for the vacuum potential of a
+  // free-boundary run. NESTOR expands the potential to (mpol, ntor) like
+  // the plasma, which limits the vacuum field on helically excursing
+  // boundaries; a value above mpol / ntor raises the potential's cutoff
+  // alone, with nzeta at least 2 * vacuum_ntor + 4. 0 (default) means
+  // "use mpol / ntor".
+  int vacuum_mpol;
+  int vacuum_ntor;
+
   // number of poloidal grid points; if odd: is rounded to next smaller even
   // number
   int ntheta;

src/vmecpp/cpp/vmecpp/common/vmec_indata/vmec_indata_test.cc

@@ -267,6 +267,45 @@ TEST(TestVmecINDATA, CheckSplineProfilesNeedKnots) {
   EXPECT_TRUE(IsConsistent(indata, /*enable_info_messages=*/false).ok());
 }
 
+// The vacuum potential's cutoffs default to the plasma's, may only exceed
+// them, and the toroidal grid has to carry the toroidal one.
+TEST(TestVmecINDATA, VacuumCutoffsStayAboveThePlasmaResolution) {
+  const absl::StatusOr<std::string> json =
+      file_io::ReadFile("vmecpp/test_data/cth_like_free_bdy.json");
+  ASSERT_TRUE(json.ok());
+  absl::StatusOr<VmecINDATA> maybe = VmecINDATA::FromJson(*json);
+  ASSERT_TRUE(maybe.ok());
+  VmecINDATA indata = *maybe;  // mpol 5, ntor 4, nzeta 36
+  EXPECT_EQ(indata.vacuum_mpol, 0);
+  EXPECT_EQ(indata.vacuum_ntor, 0);
+  EXPECT_TRUE(IsConsistent(indata, /*enable_info_messages=*/false).ok());
+
+  indata.vacuum_mpol = 4;
+  EXPECT_EQ(IsConsistent(indata, /*enable_info_messages=*/false).code(),
+            absl::StatusCode::kInvalidArgument);
+  indata.vacuum_mpol = 8;
+  EXPECT_TRUE(IsConsistent(indata, /*enable_info_messages=*/false).ok());
+
+  indata.vacuum_ntor = 3;
+  EXPECT_EQ(IsConsistent(indata, /*enable_info_messages=*/false).code(),
+            absl::StatusCode::kInvalidArgument);
+  // 2 * 16 + 4 = 36 planes fit, 2 * 17 + 4 = 38 do not
+  indata.vacuum_ntor = 16;
+  EXPECT_TRUE(IsConsistent(indata, /*enable_info_messages=*/false).ok());
+  indata.vacuum_ntor = 17;
+  EXPECT_EQ(IsConsistent(indata, /*enable_info_messages=*/false).code(),
+            absl::StatusCode::kInvalidArgument);
+
+  indata.vacuum_ntor = 8;
+  const absl::StatusOr<std::string> round_trip_json = indata.ToJson();
+  ASSERT_TRUE(round_trip_json.ok());
+  const absl::StatusOr<VmecINDATA> round_trip =
+      VmecINDATA::FromJson(*round_trip_json);
+  ASSERT_TRUE(round_trip.ok());
+  EXPECT_EQ(round_trip->vacuum_mpol, 8);
+  EXPECT_EQ(round_trip->vacuum_ntor, 8);
+}
+
 TEST(TestVmecINDATA, ToJson) {
   const absl::StatusOr<std::string> indata_json =
       ReadFile("vmecpp/test_data/cth_like_free_bdy.json");

src/vmecpp/cpp/vmecpp/vmec/handover_storage/handover_storage.cc

@@ -43,15 +43,22 @@ HandoverStorage::HandoverStorage(const Sizes* s) : s_(*s) {
   rAxis.setZero(s_.nZeta);
   zAxis.setZero(s_.nZeta);
 
-  rCC_LCFS.setZero(mnsize);
-  rSS_LCFS.setZero(mnsize);
-  zSC_LCFS.setZero(mnsize);
-  zCS_LCFS.setZero(mnsize);
+  SetVacuumCutoffs(s_.mpol, s_.ntor);
+}
+
+void HandoverStorage::SetVacuumCutoffs(int vacuum_mpol_in, int vacuum_ntor_in) {
+  vacuum_mpol = vacuum_mpol_in;
+  vacuum_ntor = vacuum_ntor_in;
+  const int vacuum_mnsize = vacuum_mpol * (vacuum_ntor + 1);
+  rCC_LCFS.setZero(vacuum_mnsize);
+  rSS_LCFS.setZero(vacuum_mnsize);
+  zSC_LCFS.setZero(vacuum_mnsize);
+  zCS_LCFS.setZero(vacuum_mnsize);
   if (s_.lasym) {
-    rSC_LCFS.setZero(mnsize);
-    rCS_LCFS.setZero(mnsize);
-    zCC_LCFS.setZero(mnsize);
-    zSS_LCFS.setZero(mnsize);
+    rSC_LCFS.setZero(vacuum_mnsize);
+    rCS_LCFS.setZero(vacuum_mnsize);
+    zCC_LCFS.setZero(vacuum_mnsize);
+    zSS_LCFS.setZero(vacuum_mnsize);
   }
 }
 

src/vmecpp/cpp/vmecpp/vmec/handover_storage/handover_storage.h

@@ -162,6 +162,11 @@ class HandoverStorage {
   Eigen::VectorXd rAxis;
   Eigen::VectorXd zAxis;
 
+  // Fourier cutoffs of the vacuum potential; the LCFS rows below are laid
+  // out as n * vacuum_mpol + m for NESTOR, which expands to these.
+  int vacuum_mpol;
+  int vacuum_ntor;
+
   // LCFS geometry for NESTOR
   Eigen::VectorXd rCC_LCFS;
   Eigen::VectorXd rSS_LCFS;
@@ -172,6 +177,10 @@ class HandoverStorage {
   Eigen::VectorXd zCC_LCFS;
   Eigen::VectorXd zSS_LCFS;
 
+  // Size the LCFS rows for a vacuum potential expanded to (vacuum_mpol,
+  // vacuum_ntor); the constructor sizes them for the plasma's cutoffs.
+  void SetVacuumCutoffs(int vacuum_mpol, int vacuum_ntor);
+
   // [nZnT] vacuum magnetic pressure |B_vac^2|/2 at the plasma boundary
   Eigen::VectorXd vacuum_magnetic_pressure;
 

src/vmecpp/cpp/vmecpp/vmec/ideal_mhd_model/ideal_mhd_model.cc

@@ -47,15 +47,16 @@ namespace {
 
 // Set m_h.rCC_LCFS etc. to the corresponding values in the FourierGeometry
 // of the last surface, also transposing m and n dimensions to make the
-// data layout what Nestor expects.
+// data layout what Nestor expects: n * vacuum_mpol + m, where the vacuum
+// potential's cutoffs may exceed the plasma's, leaving the higher modes zero.
 void HandOverBoundaryGeometry(vmecpp::HandoverStorage& m_h,
                               const vmecpp::FourierGeometry& physical_x,
                               const vmecpp::Sizes& sizes, int offset) {
   const int ntorp1 = sizes.ntor + 1;
   for (int m = 0; m < sizes.mpol; ++m) {
     for (int n = 0; n < ntorp1; ++n) {
       const int idx_mn = m * ntorp1 + n;
-      const int idx_nm = n * sizes.mpol + m;
+      const int idx_nm = n * m_h.vacuum_mpol + m;
       m_h.rCC_LCFS[idx_nm] = physical_x.rmncc[offset + idx_mn];
       m_h.zSC_LCFS[idx_nm] = physical_x.zmnsc[offset + idx_mn];
 

src/vmecpp/cpp/vmecpp/vmec/output_quantities/output_quantities.cc

@@ -4769,8 +4769,8 @@ vmecpp::WOutFileContents vmecpp::ComputeWOutFileContents(
   // fixed-boundary run leaves it empty, as Fortran VMEC does.
   const Eigen::VectorXd& vacuum_potential = handover_storage.vacuum_potential;
   if (vacuum_potential.size() > 0) {
-    const int nf = s.ntor;
-    const int mf = s.mpol + 1;
+    const int nf = handover_storage.vacuum_ntor;
+    const int mf = handover_storage.vacuum_mpol + 1;
     const int mnpd = (2 * nf + 1) * (mf + 1);
     wout.potvac = VectorXd::Zero(2 * mnpd);
     if (vacuum_potential.size() <= wout.potvac.size()) {

src/vmecpp/cpp/vmecpp/vmec/pybind11/pybind_vmec.cc

@@ -784,7 +784,9 @@ PYBIND11_MODULE(_vmecpp, m) {
           .def_readwrite("ntheta", &VmecINDATA::ntheta)
           .def_readwrite("nzeta", &VmecINDATA::nzeta)
           .def_readwrite("mpol_geometry", &VmecINDATA::mpol_geometry)
-          .def_readwrite("ntor_geometry", &VmecINDATA::ntor_geometry);
+          .def_readwrite("ntor_geometry", &VmecINDATA::ntor_geometry)
+          .def_readwrite("vacuum_mpol", &VmecINDATA::vacuum_mpol)
+          .def_readwrite("vacuum_ntor", &VmecINDATA::vacuum_ntor);
 
   // multi-grid steps
   DefEigenProperty(pyindata, "ns_array", &VmecINDATA::ns_array);

src/vmecpp/cpp/vmecpp/vmec/vmec/vmec.cc

@@ -176,10 +176,23 @@ absl::StatusOr<std::unique_ptr<Vmec>> Vmec::FromIndata(
 }
 
 // initialize based on input file contents
+namespace {
+// Fourier cutoffs of the vacuum potential: the plasma's unless raised.
+int VacuumMpol(const vmecpp::VmecINDATA& indata) {
+  return std::max(indata.vacuum_mpol, indata.mpol);
+}
+int VacuumNtor(const vmecpp::VmecINDATA& indata) {
+  return std::max(indata.vacuum_ntor, indata.ntor);
+}
+}  // namespace
+
 Vmec::Vmec(const VmecINDATA& indata, std::optional<int> max_threads,
            OutputMode verbose, InterruptCallback interrupt_callback)
     : indata_(indata),
       s_(indata_),
+      vacuum_s_(indata_.lasym, indata_.nfp, VacuumMpol(indata_),
+                VacuumNtor(indata_), s_.ntheta, s_.nZeta, VacuumMpol(indata_),
+                VacuumNtor(indata_)),
       t_(&s_),
       b_(&s_, &t_, kSignOfJacobian),
       h_(&s_),
@@ -199,11 +212,11 @@ Vmec::Vmec(const VmecINDATA& indata, std::optional<int> max_threads,
   fc_.haveToFlipTheta = b_.setupFromIndata(indata_, verbose_);
 
   if (fc_.lfreeb) {
-    // tangential Fourier resolution
+    // tangential Fourier resolution of the vacuum potential
     // 0 : ntor
-    int nf = s_.ntor;
+    int nf = vacuum_s_.ntor;
     // 0 : (mpol + 1)
-    int mf = s_.mpol + 1;
+    int mf = vacuum_s_.mpol + 1;
     int mnpd = (2 * nf + 1) * (mf + 1);
     // For lasym = true the scalar potential carries both sin(mu-nv) and
     // cos(mu-nv) coefficients, doubling the Nestor linear system to
@@ -212,6 +225,7 @@ Vmec::Vmec(const VmecINDATA& indata, std::optional<int> max_threads,
     matrixShare.setZero(mnpd_dim * mnpd_dim);
     bvecShare.setZero(mnpd_dim);
 
+    h_.SetVacuumCutoffs(vacuum_s_.mpol, vacuum_s_.ntor);
     h_.vacuum_magnetic_pressure.setZero(s_.nZnT);
     h_.initial_plasma_pressure_at_boundary.setZero(s_.nZnT);
     h_.initial_vacuum_pressure_at_boundary.setZero(s_.nZnT);
@@ -505,7 +519,7 @@ void Vmec::SetupVacuumSolvers() {
 
     if (indata_.free_boundary_method == FreeBoundaryMethod::NESTOR) {
       fb_vac_[vac_thread_id] = std::make_unique<Nestor>(
-          &s_, tp_vac_[vac_thread_id].get(), &mgrid_,
+          &vacuum_s_, tp_vac_[vac_thread_id].get(), &mgrid_,
           std::span<double>(matrixShare.data(), matrixShare.size()),
           std::span<double>(bvecShare.data(), bvecShare.size()),
           std::span<double>(h_.vacuum_magnetic_pressure.data(),
@@ -518,7 +532,7 @@ void Vmec::SetupVacuumSolvers() {
                             vacuum_reduce_slots_.size()));
     } else if (indata_.free_boundary_method == FreeBoundaryMethod::ONLY_COILS) {
       fb_vac_[vac_thread_id] = std::make_unique<OnlyCoils>(
-          &s_, tp_vac_[vac_thread_id].get(), &mgrid_,
+          &vacuum_s_, tp_vac_[vac_thread_id].get(), &mgrid_,
           std::span<double>(h_.vacuum_magnetic_pressure.data(),
                             h_.vacuum_magnetic_pressure.size()),
           std::span<double>(h_.vacuum_b_r.data(), h_.vacuum_b_r.size()),

src/vmecpp/cpp/vmecpp/vmec/vmec/vmec.h

@@ -195,6 +195,8 @@ class Vmec {
 
   VmecINDATA indata_;
   Sizes s_;
+  // Fourier cutoffs of the vacuum potential on the plasma's tangential grid
+  Sizes vacuum_s_;
   FourierBasisFastPoloidal t_;
   Boundaries b_;
   VmecConstants constants_;