Compare a coarse Nestor run against a fine only-coils run

Author: jonathanschilling

examples/compare_nestor_only_coils.py

@@ -12,16 +12,20 @@
 folder = Path("examples/data")
 
 vmecpp_nestor = h5py.File(folder / "w7x_free_bdy_vac_nestor.out.h5", "r")
-nfp = vmecpp_nestor["/wout/nfp"][()]
-ns = vmecpp_nestor["/wout/ns"][()]
-mnmax = vmecpp_nestor["/wout/mnmax"][()]
-xm = vmecpp_nestor["/wout/xm"][()]
-xn = vmecpp_nestor["/wout/xn"][()]
+ref_nfp = vmecpp_nestor["/wout/nfp"][()]
+ref_ns = vmecpp_nestor["/wout/ns"][()]
+ref_mnmax = vmecpp_nestor["/wout/mnmax"][()]
+ref_xm = vmecpp_nestor["/wout/xm"][()]
+ref_xn = vmecpp_nestor["/wout/xn"][()]
 ref_rmnc = vmecpp_nestor["/wout/rmnc"][()]
 ref_zmns = vmecpp_nestor["/wout/zmns"][()]
 
 vmecpp_only_coils = h5py.File(folder / "w7x_free_bdy_vac_only_coils.out.h5", "r")
-# assume the rest is consistent
+tst_nfp = vmecpp_only_coils["/wout/nfp"][()]
+tst_ns = vmecpp_only_coils["/wout/ns"][()]
+tst_mnmax = vmecpp_only_coils["/wout/mnmax"][()]
+tst_xm = vmecpp_only_coils["/wout/xm"][()]
+tst_xn = vmecpp_only_coils["/wout/xn"][()]
 tst_rmnc = vmecpp_only_coils["/wout/rmnc"][()]
 tst_zmns = vmecpp_only_coils["/wout/zmns"][()]
 
@@ -30,21 +34,26 @@
 theta = np.linspace(0.0, 2.0 * np.pi, ntheta)
 for phi_degrees in [0, 18, 36]:
     phi = np.deg2rad(phi_degrees)
-    kernel = np.outer(xm, theta) - np.outer(xn, phi)
-    cos_kernel = np.cos(kernel)
-    sin_kernel = np.sin(kernel)
+    
+    ref_kernel = np.outer(ref_xm, theta) - np.outer(ref_xn, phi)
+    ref_cos_kernel = np.cos(ref_kernel)
+    ref_sin_kernel = np.sin(ref_kernel)
+
+    tst_kernel = np.outer(tst_xm, theta) - np.outer(tst_xn, phi)
+    tst_cos_kernel = np.cos(tst_kernel)
+    tst_sin_kernel = np.sin(tst_kernel)
 
     plt.figure()
-    for j in [0, 2**2, 4**2, 6**2, 50]:
-        ref_r = np.dot(ref_rmnc[j, :], cos_kernel)
-        ref_z = np.dot(ref_zmns[j, :], sin_kernel)
+    for j in [0, 2**2, 4**2, 6**2, -1]:
+        ref_r = np.dot(ref_rmnc[j, :], ref_cos_kernel)
+        ref_z = np.dot(ref_zmns[j, :], ref_sin_kernel)
         if j == 0:
             plt.plot(ref_r, ref_z, "ro", label="Nestor")
         else:
             plt.plot(ref_r, ref_z, "r-", lw=2)
 
-        tst_r = np.dot(tst_rmnc[j, :], cos_kernel)
-        tst_z = np.dot(tst_zmns[j, :], sin_kernel)
+        tst_r = np.dot(tst_rmnc[j, :], tst_cos_kernel)
+        tst_z = np.dot(tst_zmns[j, :], tst_sin_kernel)
         if j == 0:
             plt.plot(tst_r, tst_z, "bx", label="only coils")
         else:

examples/data/w7x_free_bdy_vac_only_coils.json

@@ -1,12 +1,12 @@
 {
   "lasym" : false,
   "nfp" : 5,
-  "mpol" : 7,
-  "ntor" : 6,
+  "mpol" : 20,
+  "ntor" : 12,
   "ntheta" : 0,
   "nzeta" : 36,
-  "ns_array" : [ 51 ],
-  "ftol_array" : [ 1.0E-12 ],
+  "ns_array" : [ 99 ],
+  "ftol_array" : [ 1.0E-15 ],
   "niter_array" : [ 50000 ],
   "delt" : 1.0,
   "tcon0" : 1.0,
@@ -23,8 +23,8 @@
   "ac" : [ 0.0, 1.0 ],
   "curtor" : 0.0,
   "bloat" : 1.0,
-  "raxis_c" : [ 5.6343, 0.35209, 0.0, 0.0, 0.0, 0.0, 0.0 ],
-  "zaxis_s" : [ 0.0, -0.29578, 0.0, 0.0, 0.0, 0.0, 0.0 ],
+  "raxis_c" : [ 5.6343, 0.35209, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ],
+  "zaxis_s" : [ 0.0, -0.29578, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ],
   "rbc" : [ {
     "n" : 0,
     "m" : 0,