Fix to normalization of vacuum well.

src/simsopt/geo/surfaceobjectives.py

@@ -957,7 +957,7 @@ def compute(self):
         B = B.reshape((nphi, ntheta, 3))
         modB = np.sqrt(B[:, :, 0]**2 + B[:, :, 1]**2 + B[:, :, 2]**2)
         G = self.boozer_surface.res['G']
-        jac = G/modB**2
+        jac = -2*np.pi*G/modB**2
 
         self._J = np.mean(jac) 
 
@@ -968,7 +968,7 @@ def compute(self):
 
         dJ_by_dB = self.dJ_by_dB().reshape((-1, 3))
         dJ_by_dcoils = self.biotsavart.B_vjp(dJ_by_dB)
-        dJ_by_dG = self._J/G
+        dJ_by_dG = -2*np.pi*self._J/G
 
         # tack on dJ_diota = 0, dJ_dG to the end of dJ_ds
         dJ_ds = np.concatenate((self.dJ_by_dsurfacecoefficients(), [0., dJ_by_dG]))
@@ -992,7 +992,7 @@ def dJ_by_dB(self):
 
         dmodB_dB = B / modB[..., None]
         dnum_by_dB = - 2 * G * dmodB_dB /(modB[...,None]**3 * nphi * ntheta) 
-        return dnum_by_dB 
+        return -2*np.pi*dnum_by_dB 
 
     def dJ_by_dsurfacecoefficients(self):
         """
@@ -1015,7 +1015,7 @@ def dJ_by_dsurfacecoefficients(self):
         dmodB_dc = (B[:, :, 0, None] * dB_dc[:, :, 0, :] + B[:, :, 1, None] * dB_dc[:, :, 1, :] + B[:, :, 2, None] * dB_dc[:, :, 2, :])/modB[:, :, None]
 
         dnum_dc = np.mean(- 2 * G * dmodB_dc/(modB[...,None]**3), axis=(0,1 )) 
-        return dnum_dc 
+        return -2*np.pi*dnum_dc 
 
 class Iotas(Optimizable):
     """

tests/geo/test_surface_objectives.py

@@ -341,34 +341,35 @@ def df(dofs):
 
 class DifferentialVolumeTests(unittest.TestCase):
     def test_diff_vol(self):
-        nphi = 10
-        ntheta = 10
-        mpol = 6
-        ntor = 6
-        maxiter = 100
-        bs, boozer_surface_0 = get_boozer_surface(label="ToroidalFlux",nphi=nphi,ntheta=ntheta,mpol=mpol,ntor=ntor,maxiter=maxiter,lab_scale=0.01)
-        target_0 = boozer_surface_0.targetlabel
-        iota = boozer_surface_0.res['iota']
-        G0 = boozer_surface_0.res['G']
-        # dtheta/du = 2pi 
-        diffvol_0 = DifferentialVolume(boozer_surface_0, bs)
-        print(diffvol_0.J()*(2*np.pi))
-
-        for eps in [1e-4,1e-3,1e-2,1e-1,0.5]:
-            for label in ["Volume", "ToroidalFlux"]:
-                bs_1, boozer_surface_1 = get_boozer_surface(label=label,nphi=nphi,ntheta=ntheta,mpol=mpol,ntor=ntor,maxiter=maxiter,lab_scale=0.01*(1+eps))
-                bs_2, boozer_surface_2 = get_boozer_surface(label=label,nphi=nphi,ntheta=ntheta,mpol=mpol,ntor=ntor,maxiter=maxiter,lab_scale=0.01*(1-eps))
-
-                surf_1 = boozer_surface_1.surface
-                surf_2 = boozer_surface_2.surface 
-                vol_1 = Volume(surf_1)
-                vol_2 = Volume(surf_2)
-                tf_1 = ToroidalFlux(surf_1,bs_1)
-                tf_2 = ToroidalFlux(surf_2,bs_2)
-
-                diff_flux_norm = (tf_1.J()-tf_2.J())/(2*np.pi)
-                diffvol = (vol_1.J()-vol_2.J())/diff_flux_norm
-                print(diffvol)
+        nphi = 30
+        ntheta = 30
+        mpol = 7
+        ntor = 7
+        maxiter = 1000
+        diff_vols = []
+        diff_vols_eps = []
+        lab_scale = 1.0
+        eps = 1e-4 
+        label = "Volume"
+        bs, boozer_surface_0 = get_boozer_surface(label=label,nphi=nphi,ntheta=ntheta,mpol=mpol,ntor=ntor,maxiter=maxiter,lab_scale=lab_scale)
+
+        surf_0 = boozer_surface_0.surface
+        tf_0 = ToroidalFlux(surf_0,bs).J()
+        vol_0 = Volume(surf_0).J()
+        diffvol_0 = DifferentialVolume(boozer_surface_0, bs).J()
+
+        epsilons=np.power(2., -np.asarray(range(3, 7)))
+        for eps in epsilons:
+            bs_1, boozer_surface_1 = get_boozer_surface(label=label,nphi=nphi,ntheta=ntheta,mpol=mpol,ntor=ntor,maxiter=maxiter,lab_scale=lab_scale*(1-eps))
+
+            surf_1 = boozer_surface_1.surface
+            vol_1 = Volume(surf_1).J()
+            tf_1 = ToroidalFlux(surf_1,bs_1).J()
+
+            diff_flux_norm = (tf_1-tf_0)/(2*np.pi)
+            diffvol = (vol_1-vol_0)/diff_flux_norm
+            err = np.abs(diffvol - diffvol_0)/np.abs(diffvol_0)
+            assert err < 1e-3
 
     def test_diffvol_derivative(self):
         """