Use cubic interpolation for vacuum field tables

src/vmecpp/cpp/vmecpp/free_boundary/mgrid_provider/mgrid_provider.cc

@@ -387,36 +387,47 @@ absl::Status MGridProvider::interpolate(int ztMin, int ztMax, int nZeta,
     double r = std::max(minR, std::min(rLCFS[kl], maxR));
     double z = std::max(minZ, std::min(zLCFS[kl], maxZ));
 
-    // DETERMINE INTEGER INDICES (IR,JZ) FOR LOWER LEFT R, Z CORNER GRID POINT
-    int ir = static_cast<int>(floor((r - minR) / deltaR));
-    int jz = static_cast<int>(floor((z - minZ) / deltaZ));
-    int ir1 = std::min(numR - 1, ir + 1);
-    int jz1 = std::min(numZ - 1, jz + 1);
-
-    // COMPUTE RI, ZJ AND PR, QZ AT GRID POINT (IR , JZ)
-    double ri = minR + ir * deltaR;
-    double zj = minZ + jz * deltaZ;
-    double pr = (r - ri) / deltaR;
-    double qz = (z - zj) / deltaZ;
-
-    // COMPUTE WEIGHTS WIJ FOR 4 CORNER GRID POINTS
-    double w22 = pr * qz;                //    p *   q
-    double w21 = pr - w22;               //    p *(1-q) = p - p*q
-    double w12 = qz - w22;               // (1-p)*   q  = q - p*q
-    double w11 = 1.0 + w22 - (pr + qz);  // (1-p)*(1-q) = 1 + p*q - (p + q)
-
-    // COMPUTE B FIELD AT R, PHI, Z BY INTERPOLATION
-    int kj_i_ = (k * numZ + jz) * numR + ir;
-    int kj1i_ = (k * numZ + jz1) * numR + ir;
-    int kj_i1 = (k * numZ + jz) * numR + ir1;
-    int kj1i1 = (k * numZ + jz1) * numR + ir1;
-
-    m_interpBr[kl - ztMin] =
-        w11 * bR[kj_i_] + w12 * bR[kj1i_] + w21 * bR[kj_i1] + w22 * bR[kj1i1];
-    m_interpBp[kl - ztMin] =
-        w11 * bP[kj_i_] + w12 * bP[kj1i_] + w21 * bP[kj_i1] + w22 * bP[kj1i1];
-    m_interpBz[kl - ztMin] =
-        w11 * bZ[kj_i_] + w12 * bZ[kj1i_] + w21 * bZ[kj_i1] + w22 * bZ[kj1i1];
+    // Use a centred four-node stencil, shifted inward at the grid edges.
+    // Smaller tables use the polynomial supported by their available nodes.
+    const int r_nodes = std::min(4, numR);
+    const int z_nodes = std::min(4, numZ);
+    const double r_index = (r - minR) / deltaR;
+    const double z_index = (z - minZ) / deltaZ;
+    const int r_start =
+        std::clamp(static_cast<int>(floor(r_index)) - 1, 0, numR - r_nodes);
+    const int z_start =
+        std::clamp(static_cast<int>(floor(z_index)) - 1, 0, numZ - z_nodes);
+    const auto weights = [](double u, int nodes) -> std::array<double, 4> {
+      if (nodes == 4) {
+        return {-(u - 1) * (u - 2) * (u - 3) / 6, u * (u - 2) * (u - 3) / 2,
+                -u * (u - 1) * (u - 3) / 2, u * (u - 1) * (u - 2) / 6};
+      }
+      std::array<double, 4> result{};
+      for (int i = 0; i < nodes; ++i) {
+        result[i] = 1.0;
+        for (int j = 0; j < nodes; ++j) {
+          if (i != j) result[i] *= (u - j) / (i - j);
+        }
+      }
+      return result;
+    };
+    const auto r_weights = weights(r_index - r_start, r_nodes);
+    const auto z_weights = weights(z_index - z_start, z_nodes);
+    double br = 0.0;
+    double bp = 0.0;
+    double bz = 0.0;
+    for (int j = 0; j < z_nodes; ++j) {
+      for (int i = 0; i < r_nodes; ++i) {
+        const int index = (k * numZ + z_start + j) * numR + r_start + i;
+        const double weight = r_weights[i] * z_weights[j];
+        br += weight * bR[index];
+        bp += weight * bP[index];
+        bz += weight * bZ[index];
+      }
+    }
+    m_interpBr[kl - ztMin] = br;
+    m_interpBp[kl - ztMin] = bp;
+    m_interpBz[kl - ztMin] = bz;
   }  // kl
 
   absl::Status status = absl::OkStatus();

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

@@ -813,11 +813,9 @@ TEST(TestVmec, MultiGridFreeBoundary) {
   const auto output = vmecpp::run(*indata, std::nullopt, 1);
   ASSERT_TRUE(output.ok());
 
-  // Regression guard for issue #330/#640 and other changes to the multigrid
-  // convergence path. 344 with the historical unbalanced stage entry; 321
-  // since the vacuum state is seeded across multigrid transitions (the
-  // second stage enters force-balanced instead of kicking the boundary).
-  EXPECT_EQ(output->wout.niter, 321);
+  // Regression guard for the seeded vacuum state across multigrid transitions
+  // with cubic vacuum field interpolation.
+  EXPECT_EQ(output->wout.niter, 322);
 }  // MultiGridFreeBoundary
 
 // The free-boundary threed1 section covers the poloidal range the run is solved

src/vmecpp/cpp/vmecpp_large_cpp_tests/free_boundary/mgrid_provider/mgrid_provider_test.cc

@@ -4,6 +4,7 @@
 // SPDX-License-Identifier: MIT
 #include "vmecpp/free_boundary/mgrid_provider/mgrid_provider.h"
 
+#include <algorithm>
 #include <cmath>
 #include <fstream>
 #include <string>
@@ -460,4 +461,82 @@ TEST(MGridProviderValidation, LoadFileReadsCoilGroupNames) {
   }
 }
 
+TEST(MGridPolynomialInterpolation,
+     ReproducesTensorPolynomialsOnAvailableStencil) {
+  for (const int num_r : {2, 3, 4, 7}) {
+    for (const int num_z : {2, 3, 4, 8}) {
+      for (const int degree : {1, 2, 3}) {
+        makegrid::MagneticFieldResponseTable table;
+        auto& parameters = table.parameters;
+        parameters.normalize_by_currents = false;
+        parameters.number_of_field_periods = 3;
+        parameters.r_grid_minimum = 1.0;
+        parameters.r_grid_maximum = 3.0;
+        parameters.z_grid_minimum = -0.7;
+        parameters.z_grid_maximum = 0.9;
+        parameters.number_of_r_grid_points = num_r;
+        parameters.number_of_z_grid_points = num_z;
+        parameters.number_of_phi_grid_points = 5;
+        const int num_cells = num_r * num_z * 5;
+        table.b_r.resize(2, num_cells);
+        table.b_p.resize(2, num_cells);
+        table.b_z.resize(2, num_cells);
+        const auto polynomial = [degree](double r, double z, int k) {
+          return std::pow(r, degree) + 2 * std::pow(z, degree) +
+                 std::pow(r * z, degree) + 0.3 * r * z + k;
+        };
+        for (int k = 0; k < 5; ++k) {
+          for (int j = 0; j < num_z; ++j) {
+            for (int i = 0; i < num_r; ++i) {
+              const double value = polynomial(1.0 + 2.0 * i / (num_r - 1),
+                                              -0.7 + 1.6 * j / (num_z - 1), k);
+              const int index = (k * num_z + j) * num_r + i;
+              table.b_r(0, index) = value;
+              table.b_p(0, index) = 2 * value;
+              table.b_z(0, index) = -value;
+              table.b_r(1, index) = 3 * value;
+              table.b_p(1, index) = 6 * value;
+              table.b_z(1, index) = -3 * value;
+            }
+          }
+        }
+        Eigen::VectorXd currents(2);
+        currents << 2.0, -0.25;
+        SCOPED_TRACE(
+            absl::StrFormat("nr=%d nz=%d polynomial=%d", num_r, num_z, degree));
+        MGridProvider provider;
+        ASSERT_TRUE(provider.LoadFields(table, currents).ok());
+        Eigen::VectorXd r(205), z(205), br(205), bp(205), bz(205);
+        for (int i = 0; i < 205; ++i) {
+          r[i] = 1.0 + 2.0 * (((i * 37) % 205) + 0.5) / 205.0;
+          z[i] = -0.7 + 1.6 * (((i * 71) % 205) + 0.5) / 205.0;
+        }
+        r[0] = 1.0;
+        z[0] = -0.7;
+        r[1] = 3.0;
+        z[1] = 0.9;
+        r[2] = 3.0;
+        z[2] = -0.7;
+        r[3] = 1.0;
+        z[3] = 0.9;
+        ASSERT_TRUE(
+            provider.interpolate(0, 205, 5, 205, r, z, br, bp, bz).ok());
+        double error = 0.0;
+        for (int i = 0; i < 205; ++i) {
+          const double expected = 1.25 * polynomial(r[i], z[i], i % 5);
+          error = std::max({error, std::abs(br[i] - expected),
+                            std::abs(bp[i] - 2 * expected),
+                            std::abs(bz[i] + expected)});
+        }
+        const bool exact = num_r > degree && num_z > degree;
+        if (exact) {
+          EXPECT_LT(error, 1e-11);
+        } else {
+          EXPECT_GT(error, 1e-4);
+        }
+      }
+    }
+  }
+}
+
 }  // namespace vmecpp