Added pres to output.

src/_booz_xform/bindings.cpp

@@ -62,6 +62,7 @@ need to call this function directly.
   be specified.
 :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)
 )",
 	 "ns"_a,
 	 "iotas"_a,
@@ -78,7 +79,8 @@ need to call this function directly.
 	 "bsubvmnc0"_a,
 	 "bsubvmns0"_a,
      "phips"_a=defaultInitPtr,
-     "chi"_a=defaultInitPtr)
+     "chi"_a=defaultInitPtr,
+     "pres"_a=defaultInitPtr)
 
     .def("run", &Booz_xform::run, R"(
 Run the transformation to Boozer coordinates on all the selected
@@ -187,6 +189,9 @@ evaluated on all the magnetic surfaces for which input data was provided.)")
 (1D float array of length ns_in, input) Toroidal flux normalized by 2*pi,
 evaluated on all the magnetic surfaces for which input data was provided.)")
 
+    .def_readwrite("pres", &Booz_xform::pres, R"(
+(1D float array of length ns_in, input) Pressure on full vmec grid.)")
+
     .def_readwrite("rmnc", &Booz_xform::rmnc, R"(
 (2D float array of size mnmax x ns_in, input) cos(m * theta_0 - n *
 zeta_0) Fourier modes of the major radius R.)")

src/_booz_xform/booz_xform.hpp

@@ -401,6 +401,8 @@ namespace booz_xform {
     */
     Vector chi;
 
+    /** Pressure evaluated on the full vmec grid **/
+    Vector pres;
 
     //! Constructor
     /**
@@ -441,7 +443,8 @@ namespace booz_xform {
 			Matrix& bsubvmnc,
 			Matrix& bsubvmns,
       Vector& phips=defaultInitPtr,
-      Vector& chi=defaultInitPtr);
+      Vector& chi=defaultInitPtr,
+      Vector& pres=defaultInitPtr);
 
     //! Carry out the transformation calculation
     /**

src/_booz_xform/init_from_vmec.cpp

@@ -23,7 +23,8 @@ void Booz_xform::init_from_vmec(int ns,
 				Matrix& bsubvmnc0,
 				Matrix& bsubvmns0,
         Vector& phip0,
-        Vector& chi0) {
+        Vector& chi0,
+        Vector& pres0) {
   int j, k;
   ns_in = ns - 1;
 
@@ -36,6 +37,7 @@ void Booz_xform::init_from_vmec(int ns,
   if (!skip_phip0) {
       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 (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");
@@ -97,13 +99,16 @@ void Booz_xform::init_from_vmec(int ns,
   if (!skip_phip0) {
       phip.resize(ns_in + 1);
       chi.resize(ns_in + 1);
+      pres.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];
       }
   } else {
      phip.resize(0);
      chi.resize(0);
+     pres.resize(0);
   }
   // By default, prepare to do the Boozer transformation at all
   // half-grid surfaces:

src/_booz_xform/read_boozmn.cpp

@@ -41,12 +41,14 @@ void Booz_xform::read_boozmn(std::string filename) {
   phi.resize(ns_in+1);
   phip.resize(ns_in+1);
   chi.resize(ns_in+1);
+  pres.resize(ns_in+1);
   nc.get("iota_b", iota_in);
   nc.get("bvco_b", Boozer_G_in);
   nc.get("buco_b", Boozer_I_in);
   nc.get("phi_b", phi);
   nc.get("phip_b", phip);
   nc.get("chi_b", chi);
+  nc.get("pres_b", pres);
   for (j = 0; j < ns_in; j++) {
       iota[j] = iota_in[j+1];
       Boozer_G_all[j] = Boozer_G_in[j+1];

src/_booz_xform/read_wout.cpp

@@ -29,11 +29,14 @@ void Booz_xform::read_wout(std::string filename, bool flux) {
 
   Vector phip0;
   Vector chi0;
+  Vector pres0;
   if (flux) {
       phip0.resize(ns);
       chi0.resize(ns);
+      pres0.resize(ns);
       nc.get("chi", chi0);
       nc.get("phipf", phip0);
+      nc.get("pres",pres0);
   }
 
   phi.resize(ns);
@@ -103,7 +106,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);
+           bsubumnc0, bsubumns0, bsubvmnc0, bsubvmns0, phip0, chi0, pres0);
   } else {
       init_from_vmec(ns, iotas, rmnc0, rmns0, zmnc0, zmns0,
            lmnc0, lmns0, bmnc0, bmns0,

src/_booz_xform/write_boozmn.cpp

@@ -81,21 +81,21 @@ void Booz_xform::write_boozmn(std::string filename) {
   // available here. Just write vectors of zeros for
   // backwards-compatibility of the boozmn.nc files.
   Vector pres_b(ns_in + 1), beta_b(ns_in + 1);
-  pres_b.setZero();
-  nc.put(radius_dim, "pres_b", pres_b, placeholder_str, "dimensionless");
   beta_b.setZero();
   nc.put(radius_dim, "beta_b", pres_b, placeholder_str, "dimensionless");
-  Vector phip_dummy(ns_in + 1);
-  Vector chi_dummy(ns_in + 1);
+  Vector phip_dummy(ns_in + 1), chi_dummy(ns_in + 1);
   phip_dummy.setZero();
   chi_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");
   } 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");
 

tests/test_regression.py

@@ -72,7 +72,7 @@ def test_regression(self):
                 vars = f.variables.keys()
                 # These variables will not match:
                 if flux:
-                    exclude = ['rmax_b', 'rmin_b', 'betaxis_b', 'version', 'pres_b', 'beta_b']
+                    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']
                 for var in vars:

tests/test_write_read.py

@@ -64,6 +64,7 @@ def test_write_read(self):
                     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)
                     else:
                         np.testing.assert_equal(b2.phip,0)
                         assert len(b1.phip) == 0