Fix the sign of the cross terms in NESTOR's analytic add-back

benchmarks/test_benchmarks.py

@@ -154,7 +154,7 @@ def test_bench_free_boundary(benchmark, free_boundary_input, response_table):
         rounds=3,
         warmup_rounds=1,
     )
-    assert result.wout.volume == pytest.approx(0.3075, rel=1e-3)
+    assert result.wout.volume == pytest.approx(0.3070, rel=1e-3)
 
 
 # ---------------------------------------------------------------------------

src/vmecpp/cpp/vmecpp/free_boundary/nestor/BUILD.bazel

@@ -18,3 +18,17 @@ cc_library(
         "//vmecpp/free_boundary/free_boundary_base:free_boundary_base",
     ],
 )
+
+cc_test(
+    name = "nestor_test",
+    srcs = ["nestor_test.cc"],
+    deps = [
+        ":nestor",
+        "//vmecpp/common/fourier_basis_fast_toroidal",
+        "//vmecpp/common/sizes:sizes",
+        "//vmecpp/free_boundary/mgrid_provider:mgrid_provider",
+        "//vmecpp/free_boundary/tangential_partitioning:tangential_partitioning",
+        "@googletest//:gtest_main",
+    ],
+    size = "small",
+)

src/vmecpp/cpp/vmecpp/free_boundary/nestor/nestor_test.cc

@@ -0,0 +1,158 @@
+// SPDX-FileCopyrightText: 2024-present Proxima Fusion GmbH
+// <info@proximafusion.com>
+//
+// SPDX-License-Identifier: MIT
+#include "vmecpp/free_boundary/nestor/nestor.h"
+
+#include <cmath>
+#include <vector>
+
+#include "gtest/gtest.h"
+#include "vmecpp/common/fourier_basis_fast_toroidal/fourier_basis_fast_toroidal.h"
+#include "vmecpp/common/sizes/sizes.h"
+#include "vmecpp/free_boundary/mgrid_provider/mgrid_provider.h"
+#include "vmecpp/free_boundary/tangential_partitioning/tangential_partitioning.h"
+
+namespace vmecpp {
+
+// Exterior Neumann problem with a known solution: the external field is that
+// of a line of vertical point dipoles along the circular magnetic axis of a
+// helically deformed torus, i.e. an interior source with a single-valued
+// potential that is harmonic in the whole exterior. The vacuum field that
+// NESTOR reconstructs must vanish on the surface, so |B|^2/2 there is the
+// solver's error, measured against the external field's own |B_ext|^2/2.
+// The field is supplied exactly on the surface (no mgrid interpolation).
+TEST(TestNestor, InteriorDipoleLineGivesZeroFieldOnHelicalTorus) {
+  const bool lasym = false;
+  const int nfp = 2;
+  const int mpol = 8;
+  const int ntor = 4;
+  const int ntheta = 0;
+  const int nzeta = 24;
+  Sizes s(lasym, nfp, mpol, ntor, ntheta, nzeta);
+  FourierBasisFastToroidal fb(&s);
+
+  // circular torus R0 = 1, a = 0.3 with a helical (m = 1, n = 1) deformation
+  const double R0 = 1.0;
+  const double a = 0.3;
+  const double eps = 0.05;
+  std::vector<double> rmnc(s.mnmax, 0.0);
+  std::vector<double> zmns(s.mnmax, 0.0);
+  std::vector<int> xm(s.mnmax);
+  std::vector<int> xn(s.mnmax);
+  {
+    int mn = 0;
+    for (int n = 0; n < s.ntor + 1; ++n) {
+      xm[mn] = 0;
+      xn[mn] = n;
+      ++mn;
+    }
+    for (int m = 1; m < s.mpol; ++m) {
+      for (int n = -s.ntor; n < s.ntor + 1; ++n) {
+        xm[mn] = m;
+        xn[mn] = n;
+        ++mn;
+      }
+    }
+    ASSERT_EQ(mn, s.mnmax);
+    for (mn = 0; mn < s.mnmax; ++mn) {
+      if (xm[mn] == 0 && xn[mn] == 0) rmnc[mn] = R0;
+      if (xm[mn] == 1 && xn[mn] == 0) {
+        rmnc[mn] = a;
+        zmns[mn] = a;
+      }
+      if (xm[mn] == 1 && xn[mn] == 1) {
+        rmnc[mn] = eps;
+        zmns[mn] = eps;
+      }
+    }
+  }
+  std::vector<double> rCC(s.mnsize), rSS(s.mnsize), zSC(s.mnsize),
+      zCS(s.mnsize);
+  std::vector<double> rSC, rCS, zCC, zSS;
+  fb.cos_to_cc_ss(rmnc, rCC, rSS, s.ntor, s.mpol);
+  fb.sin_to_sc_cs(zmns, zSC, zCS, s.ntor, s.mpol);
+
+  // surface points on NESTOR's grid: kl = l * nZeta + k, theta_l = 2 pi l /
+  // nThetaEven, zeta_k = 2 pi k / (nfp nZeta)
+  Eigen::VectorXd br = Eigen::VectorXd::Zero(s.nZnT);
+  Eigen::VectorXd bp = Eigen::VectorXd::Zero(s.nZnT);
+  Eigen::VectorXd bz = Eigen::VectorXd::Zero(s.nZnT);
+  const int num_dipoles = 96 * nfp;
+  const double moment = 1.0 / 96.0;  // total vertical moment 1 A m^2 per period
+  double ref_sum = 0.0;
+  for (int kl = 0; kl < s.nZnT; ++kl) {
+    const int l = kl / s.nZeta;
+    const int k = kl % s.nZeta;
+    const double theta = 2.0 * M_PI * l / s.nThetaEven;
+    const double zeta = 2.0 * M_PI * k / (nfp * s.nZeta);
+    double R = 0.0, Z = 0.0;
+    for (int mn = 0; mn < s.mnmax; ++mn) {
+      const double arg = xm[mn] * theta - xn[mn] * nfp * zeta;
+      R += rmnc[mn] * std::cos(arg);
+      Z += zmns[mn] * std::sin(arg);
+    }
+    const double x = R * std::cos(zeta);
+    const double y = R * std::sin(zeta);
+    double bx = 0.0, by = 0.0, bzc = 0.0;
+    for (int j = 0; j < num_dipoles; ++j) {
+      const double phi_j = 2.0 * M_PI * j / num_dipoles;
+      const double dx = x - R0 * std::cos(phi_j);
+      const double dy = y - R0 * std::sin(phi_j);
+      const double dz = Z;
+      const double r2 = dx * dx + dy * dy + dz * dz;
+      const double r = std::sqrt(r2);
+      const double f = 1.0e-7 / (r2 * r);  // mu0 / (4 pi) / r^3
+      const double mdotn = moment * dz / r;
+      bx += f * 3.0 * mdotn * dx / r;
+      by += f * 3.0 * mdotn * dy / r;
+      bzc += f * (3.0 * mdotn * dz / r - moment);
+    }
+    br[kl] = bx * std::cos(zeta) + by * std::sin(zeta);
+    bp[kl] = -bx * std::sin(zeta) + by * std::cos(zeta);
+    bz[kl] = bzc;
+    const double half_b2 = 0.5 * (bx * bx + by * by + bzc * bzc);
+    ref_sum += half_b2 * half_b2;
+  }
+  const double ref_rms = std::sqrt(ref_sum / s.nZnT);
+
+  MGridProvider mgrid;
+  mgrid.SetFixedMagneticField(br, bp, bz);
+
+  TangentialPartitioning tp(s.nZnT);
+  const int nf = s.ntor;
+  const int mf = s.mpol + 1;
+  const int mnpd = (2 * nf + 1) * (mf + 1);
+  std::vector<double> matrixShare(mnpd * mnpd);
+  std::vector<double> bvecShare(mnpd);
+  std::vector<double> bSqVacShare(s.nZnT);
+  std::vector<double> vacuum_b_r(s.nZnT), vacuum_b_phi(s.nZnT),
+      vacuum_b_z(s.nZnT);
+  std::vector<double> reduce_slots(matrixShare.size());
+  Eigen::PartialPivLU<Eigen::MatrixXd> lu_decomposition;
+  Nestor nestor(&s, &tp, &mgrid, matrixShare, bvecShare, bSqVacShare,
+                &lu_decomposition, vacuum_b_r, vacuum_b_phi, vacuum_b_z,
+                reduce_slots);
+
+  std::vector<double> axis_r(s.nZeta, R0);
+  std::vector<double> axis_z(s.nZeta, 0.0);
+  double bSubUVac = 0.0;
+  double bSubVVac = 0.0;
+  const int kSignOfJacobian = -1;
+  const auto status = nestor.update(rCC, rSS, rSC, rCS, zSC, zCS, zCC, zSS,
+                                    kSignOfJacobian, axis_r, axis_z, &bSubUVac,
+                                    &bSubVVac, /*netToroidalCurrent=*/0.0,
+                                    /*ivacskip=*/0);
+  ASSERT_TRUE(status.ok()) << status.status();
+
+  double err_sum = 0.0;
+  for (int kl = 0; kl < s.nZnT; ++kl) {
+    err_sum += bSqVacShare[kl] * bSqVacShare[kl];
+  }
+  const double err_rms = std::sqrt(err_sum / s.nZnT);
+  // 4.9e-3 with the cross-term sign defect in the analytic add-back
+  EXPECT_LT(err_rms / ref_rms, 1.0e-4)
+      << "|B|^2/2 residual relative to |B_ext|^2/2: " << err_rms / ref_rms;
+}
+
+}  // namespace vmecpp

src/vmecpp/cpp/vmecpp/free_boundary/singular_integrals/singular_integrals.cc

@@ -430,6 +430,17 @@ void SingularIntegrals::performUpdate(const Eigen::VectorXd& bDotN,
 
         } else {
           // analysum2
+          //
+          // T^- and S^- carry sin(mu - |n|v) (mode (m, +|n|)) and T^+ and S^+
+          // carry sin(mu + |n|v) (mode (m, -|n|)): this is the assignment
+          // under which the analytic Fourier coefficients equal those of the
+          // tangent-plane kernels subtracted in RegularizedIntegrals, with the
+          // metric and curvature cross terms as stored. The opposite
+          // assignment (the one PARVMEC's analyt makes by passing slm, tlm,
+          // slp, tlp to analysum2's dummy arguments slp, tlp, slm, tlm) gives
+          // mode (m, n) the coefficient of (m, -n): an error that is first
+          // order in the non-axisymmetric shaping and independent of the
+          // resolution.
 
           for (int kl = tp_.ztMin; kl < tp_.ztMax; ++kl) {
             const int l = kl / s_.nZeta;
@@ -468,10 +479,10 @@ void SingularIntegrals::performUpdate(const Eigen::VectorXd& bDotN,
               const double sinp3 = coeff1 * fb_.cosnv[idx_nk + 3] -
                                    coeff2 * fb_.sinnv[idx_nk + 3];
 
-              buf_m_posn[0] += Tlp[fl][klRel + 0] * c0 * sinp0;
-              buf_m_posn[1] += Tlp[fl][klRel + 1] * c1 * sinp1;
-              buf_m_posn[2] += Tlp[fl][klRel + 2] * c2 * sinp2;
-              buf_m_posn[3] += Tlp[fl][klRel + 3] * c3 * sinp3;
+              buf_m_posn[0] += Tlm[fl][klRel + 0] * c0 * sinp0;
+              buf_m_posn[1] += Tlm[fl][klRel + 1] * c1 * sinp1;
+              buf_m_posn[2] += Tlm[fl][klRel + 2] * c2 * sinp2;
+              buf_m_posn[3] += Tlm[fl][klRel + 3] * c3 * sinp3;
 
               // sin(mu + |n|v) * cmns(l,n,m)
               const double sinm0 = coeff1 * fb_.cosnv[idx_nk + 0] +
@@ -483,29 +494,29 @@ void SingularIntegrals::performUpdate(const Eigen::VectorXd& bDotN,
               const double sinm3 = coeff1 * fb_.cosnv[idx_nk + 3] +
                                    coeff2 * fb_.sinnv[idx_nk + 3];
 
-              buf_m_negn[0] += Tlm[fl][klRel + 0] * c0 * sinm0;
-              buf_m_negn[1] += Tlm[fl][klRel + 1] * c1 * sinm1;
-              buf_m_negn[2] += Tlm[fl][klRel + 2] * c2 * sinm2;
-              buf_m_negn[3] += Tlm[fl][klRel + 3] * c3 * sinm3;
+              buf_m_negn[0] += Tlp[fl][klRel + 0] * c0 * sinm0;
+              buf_m_negn[1] += Tlp[fl][klRel + 1] * c1 * sinm1;
+              buf_m_negn[2] += Tlp[fl][klRel + 2] * c2 * sinm2;
+              buf_m_negn[3] += Tlp[fl][klRel + 3] * c3 * sinm3;
 
               if (fullUpdate) {
                 grpmn_sin[idx_m_posn * numLocal + klRel + 0] +=
-                    Slp[fl][klRel + 0] * sinp0;
+                    Slm[fl][klRel + 0] * sinp0;
                 grpmn_sin[idx_m_posn * numLocal + klRel + 1] +=
-                    Slp[fl][klRel + 1] * sinp1;
+                    Slm[fl][klRel + 1] * sinp1;
                 grpmn_sin[idx_m_posn * numLocal + klRel + 2] +=
-                    Slp[fl][klRel + 2] * sinp2;
+                    Slm[fl][klRel + 2] * sinp2;
                 grpmn_sin[idx_m_posn * numLocal + klRel + 3] +=
-                    Slp[fl][klRel + 3] * sinp3;
+                    Slm[fl][klRel + 3] * sinp3;
 
                 grpmn_sin[idx_m_negn * numLocal + klRel + 0] +=
-                    Slm[fl][klRel + 0] * sinm0;
+                    Slp[fl][klRel + 0] * sinm0;
                 grpmn_sin[idx_m_negn * numLocal + klRel + 1] +=
-                    Slm[fl][klRel + 1] * sinm1;
+                    Slp[fl][klRel + 1] * sinm1;
                 grpmn_sin[idx_m_negn * numLocal + klRel + 2] +=
-                    Slm[fl][klRel + 2] * sinm2;
+                    Slp[fl][klRel + 2] * sinm2;
                 grpmn_sin[idx_m_negn * numLocal + klRel + 3] +=
-                    Slm[fl][klRel + 3] * sinm3;
+                    Slp[fl][klRel + 3] * sinm3;
               }
             }
 
@@ -532,14 +543,14 @@ void SingularIntegrals::performUpdate(const Eigen::VectorXd& bDotN,
                 const double sinp = coeff1 - coeff2;
 
                 const double c = bDotN[klRel] * s_.wInt[l];
-                bvec_sin[idx_m_posn] += Tlp[fl][klRel] * c * sinp;
-                bvec_sin[idx_m_negn] += Tlm[fl][klRel] * c * sinm;
+                bvec_sin[idx_m_posn] += Tlm[fl][klRel] * c * sinp;
+                bvec_sin[idx_m_negn] += Tlp[fl][klRel] * c * sinm;
 
                 if (fullUpdate) {
                   grpmn_sin[idx_m_posn * numLocal + klRel] +=
-                      Slp[fl][klRel] * sinp;
+                      Slm[fl][klRel] * sinp;
                   grpmn_sin[idx_m_negn * numLocal + klRel] +=
-                      Slm[fl][klRel] * sinm;
+                      Slp[fl][klRel] * sinm;
                 }
               }
             }
@@ -573,14 +584,14 @@ void SingularIntegrals::performUpdate(const Eigen::VectorXd& bDotN,
               const double cosp = coeff1 + coeff2;
 
               bvec_cos[idx_m_posn] +=
-                  Tlp[fl][klRel] * bDotN[klRel] * s_.wInt[l] * cosp;
+                  Tlm[fl][klRel] * bDotN[klRel] * s_.wInt[l] * cosp;
               bvec_cos[idx_m_negn] +=
-                  Tlm[fl][klRel] * bDotN[klRel] * s_.wInt[l] * cosm;
+                  Tlp[fl][klRel] * bDotN[klRel] * s_.wInt[l] * cosm;
               if (fullUpdate) {
                 grpmn_cos[idx_m_posn * numLocal + klRel] +=
-                    Slp[fl][klRel] * cosp;
+                    Slm[fl][klRel] * cosp;
                 grpmn_cos[idx_m_negn * numLocal + klRel] +=
-                    Slm[fl][klRel] * cosm;
+                    Slp[fl][klRel] * cosm;
               }
             }
           }  // kl