Making phi reading/writing consistent with other flux functions.

src/_booz_xform/bindings.cpp

@@ -63,6 +63,7 @@ need to call this function directly.
 :param phip: Phip on vmec's full grid. (Defaults to empty Vector)
 :param chi: chi on vmec's full grid. (Defaults to empty Vector)
 :param pres: pressure on vmec's full grid. (Defaults to empty Vector)
+:param phi: phi on vmec's full grid. (Defaults to empty Vector)
 )",
 	 "ns"_a,
 	 "iotas"_a,
@@ -80,7 +81,8 @@ need to call this function directly.
 	 "bsubvmns0"_a,
      "phips"_a=defaultInitPtr,
      "chi"_a=defaultInitPtr,
-     "pres"_a=defaultInitPtr)
+     "pres"_a=defaultInitPtr,
+     "phi"_a=defaultInitPtr)
 
     .def("run", &Booz_xform::run, R"(
 Run the transformation to Boozer coordinates on all the selected

src/_booz_xform/booz_xform.hpp

@@ -444,7 +444,8 @@ namespace booz_xform {
 			Matrix& bsubvmns,
       Vector& phips=defaultInitPtr,
       Vector& chi=defaultInitPtr,
-      Vector& pres=defaultInitPtr);
+      Vector& pres=defaultInitPtr,
+      Vector& phi=defaultInitPtr);
 
     //! Carry out the transformation calculation
     /**

src/_booz_xform/init_from_vmec.cpp

@@ -24,7 +24,8 @@ void Booz_xform::init_from_vmec(int ns,
 				Matrix& bsubvmns0,
         Vector& phip0,
         Vector& chi0,
-        Vector& pres0) {
+        Vector& pres0,
+        Vector& phi0) {
   int j, k;
   ns_in = ns - 1;
 
@@ -38,6 +39,7 @@ void Booz_xform::init_from_vmec(int ns,
       if (phip0.size() != ns) throw std::runtime_error("phip0.size() is not ns");
       if (chi0.size() != ns) throw std::runtime_error("chi0.size() is not ns");
       if (pres0.size() != ns) throw std::runtime_error("pres0.size() is not ns");
+      if (phi0.size() != ns) throw std::runtime_error("phi0.size() is not ns");
   }
   if (xm.size() != mnmax) throw std::runtime_error("Size of xm is not mnmax");
   if (xn.size() != mnmax) throw std::runtime_error("Size of xn is not mnmax");
@@ -100,15 +102,19 @@ void Booz_xform::init_from_vmec(int ns,
       phip.resize(ns_in + 1);
       chi.resize(ns_in + 1);
       pres.resize(ns_in + 1);
+      phi.resize(ns_in + 1);
       for (j = 0; j < ns_in + 1; j++)  {
           phip[j] = -phip0[j]/(2*pi);
           chi[j] = chi0[j];
           pres[j] = pres0[j];
+          phi[j] = phi0[j];
       }
+      toroidal_flux = phi[ns - 1];
   } else {
      phip.resize(0);
      chi.resize(0);
      pres.resize(0);
+     phi.resize(0);
   }
   // By default, prepare to do the Boozer transformation at all
   // half-grid surfaces:

src/_booz_xform/read_wout.cpp

@@ -30,19 +30,18 @@ void Booz_xform::read_wout(std::string filename, bool flux) {
   Vector phip0;
   Vector chi0;
   Vector pres0;
+  Vector phi0;
   if (flux) {
       phip0.resize(ns);
       chi0.resize(ns);
       pres0.resize(ns);
+      phi0.resize(ns);
       nc.get("chi", chi0);
       nc.get("phipf", phip0);
       nc.get("pres",pres0);
+      nc.get("phi", phi0);
   }
 
-  phi.resize(ns);
-  nc.get("phi", phi);
-  toroidal_flux = phi[ns - 1];
-
   Vector iotas;
   iotas.resize(ns);
 
@@ -106,7 +105,7 @@ void Booz_xform::read_wout(std::string filename, bool flux) {
   if (flux) {
       init_from_vmec(ns, iotas, rmnc0, rmns0, zmnc0, zmns0,
            lmnc0, lmns0, bmnc0, bmns0,
-           bsubumnc0, bsubumns0, bsubvmnc0, bsubvmns0, phip0, chi0, pres0);
+           bsubumnc0, bsubumns0, bsubvmnc0, bsubvmns0, phip0, chi0, pres0, phi0);
   } else {
       init_from_vmec(ns, iotas, rmnc0, rmns0, zmnc0, zmns0,
            lmnc0, lmns0, bmnc0, bmns0,

src/_booz_xform/write_boozmn.cpp

@@ -83,21 +83,23 @@ void Booz_xform::write_boozmn(std::string filename) {
   Vector pres_b(ns_in + 1), beta_b(ns_in + 1);
   beta_b.setZero();
   nc.put(radius_dim, "beta_b", pres_b, placeholder_str, "dimensionless");
-  Vector phip_dummy(ns_in + 1), chi_dummy(ns_in + 1);
+  Vector phip_dummy(ns_in + 1), chi_dummy(ns_in + 1), phi_dummy(ns_in+1);
   phip_dummy.setZero();
   chi_dummy.setZero();
+  phi_dummy.setZero();
   if (phip.size()==0) {
       pres_b.setZero();
       nc.put(radius_dim, "phip_b", phip_dummy, placeholder_str, "Tesla * meter^2");
       nc.put(radius_dim, "chi_b", chi_dummy, placeholder_str, "Tesla * meter^2");
       nc.put(radius_dim, "pres_b", pres_b, placeholder_str, "dimensionless");
+      nc.put(radius_dim, "phi_b", phi_dummy, placeholder_str, "Tesla * meter^2");
   } else {
       phip[0] = 0;
       nc.put(radius_dim, "phip_b", phip, "Derivative of toroidal flux (not divided by (2*pi)) with respect to s. This grid generally does not correspond to the radial grid used for other quantities in this file!", "Tesla * meter^2");
       nc.put(radius_dim, "chi_b", chi, "Uniformly spaced grid going from 0 to the boundary poloidal flux (not divided by (2*pi)). This grid generally does not correspond to the radial grid used for other quantities in this file!", "Tesla * meter^2");
       nc.put(radius_dim, "pres_b", pres, "Pressure on full vmec grid. This grid generally does not correspond to the radial grid used for other quantities in this file!", "Pascal");
+      nc.put(radius_dim, "phi_b", phi, "Uniformly spaced grid going from 0 to the boundary toroidal flux (not divided by (2*pi)). This grid generally does not correspond to the radial grid used for other quantities in this file!", "Tesla * meter^2");
   }
-  nc.put(radius_dim, "phi_b", phi, "Uniformly spaced grid going from 0 to the boundary toroidal flux (not divided by (2*pi)). This grid generally does not correspond to the radial grid used for other quantities in this file!", "Tesla * meter^2");
 
   // ND arrays for N > 1:
   std::vector<dim_id_type> bmnc_dim;

tests/test_regression.py

@@ -74,7 +74,7 @@ def test_regression(self):
                 if flux:
                     exclude = ['rmax_b', 'rmin_b', 'betaxis_b', 'version', 'beta_b']
                 else:
-                    exclude = ['rmax_b', 'rmin_b', 'betaxis_b', 'version', 'pres_b', 'beta_b', 'phip_b']
+                    exclude = ['rmax_b', 'rmin_b', 'betaxis_b', 'version', 'pres_b', 'beta_b', 'phip_b', 'phi_b']
                 for var in vars:
                     if var in exclude:
                         continue

tests/test_write_read.py

@@ -60,16 +60,20 @@ def test_write_read(self):
                     np.testing.assert_allclose(b1.iota, b2.iota, rtol=rtol, atol=atol)
                     np.testing.assert_allclose(b1.Boozer_G_all, b2.Boozer_G_all, rtol=rtol, atol=atol)
                     np.testing.assert_allclose(b1.Boozer_I_all, b2.Boozer_I_all, rtol=rtol, atol=atol)
-                    np.testing.assert_allclose(b1.phi,b2.phi, rtol=rtol, atol=atol)
                     if flux:
                         np.testing.assert_allclose(b1.phip,b2.phip, rtol=rtol, atol=atol)
                         np.testing.assert_allclose(b1.chi,b2.chi, rtol=rtol, atol=atol)
                         np.testing.assert_allclose(b1.pres,b2.pres, rtol=rtol, atol=atol)
+                        np.testing.assert_allclose(b1.phi,b2.phi, rtol=rtol, atol=atol)
                     else:
                         np.testing.assert_equal(b2.phip,0)
                         assert len(b1.phip) == 0
                         np.testing.assert_equal(b2.chi,0)
                         assert len(b1.chi) == 0
+                        np.testing.assert_equal(b2.pres,0)
+                        assert len(b1.pres) == 0
+                        np.testing.assert_equal(b2.phi,0)
+                        assert len(b1.phi) == 0
                     os.remove(boozmn_filename)
 
 if __name__ == '__main__':