Add Equilibrium.plasmaVolume() function

Calculates the volume of the plasma by integrating the volume element
dV = 2*pi*R*dR*dZ with the `romb` method.

03-mast.py

@@ -53,6 +53,8 @@
 
 print("Plasma current: %e Amps" % (eq.plasmaCurrent()))
 print("Pressure on axis: %e Pascals" % (eq.pressure(0.0)))
+print("Plasma poloidal beta: %e" % (eq.poloidalBeta()))
+print("Plasma volume: %e m^3" % (eq.plasmaVolume()))
 
 eq.tokamak.printCurrents()
 

08-mast-upgrade.py

@@ -40,31 +40,12 @@
 isoflux = [(Rx,-Zx, Rmid, 0.0)   # Outboard midplane, lower X-point
            ,(Rx,Zx, Rmid, 0.0)   # Outboard midplane, upper X-point
 
+           # Link inner and outer midplane locations
            ,(Rmid, 0.0, Rin, 0.0)
-
-           # ,(Rx,-Zx, 0.73, -1.6)
-           # ,(Rx, Zx, 0.73,  1.6)
-
-           # ,(Rmid, 0.0, 0.73, -1.6)
-           # ,(Rmid, 0.0, 0.73,  1.6)
-
-           # ,(Rx,-Zx, 0.45, -0.95)
-           # ,(Rx, Zx, 0.45,  0.95)
-
-           # ,(Rin, 0.0, 0.45, -0.95)
-           # ,(Rin, 0.0, 0.45,  0.95)
-           # ,(Rin, 0.0, 0.45, -0.95)
-           # ,(Rin, 0.0, 0.45,  0.95)
-           # ,(Rin, 0.0, 0.45, -0.95)
-           # ,(Rin, 0.0, 0.45,  0.95)
 
-
-            ,(Rx,-Zx, 0.95, -1.77)
-            ,(Rx, Zx, 0.95,  1.77)
-
-           # ,(Rx,-Zx, 1.25, -1.9)
-           # ,(Rx, Zx, 1.25,  1.9)
-
+           # Separatrix in the divertor chamber
+           ,(Rx,-Zx, 0.95, -1.77)
+           ,(Rx, Zx, 0.95,  1.77)
            ]
 
 constrain = freegs.control.constrain(xpoints=xpoints, gamma=1e-5, isoflux=isoflux)
@@ -79,6 +60,10 @@
              constrain,   # Plasma control constraints
              show=False)   # Shows results at each nonlinear iteration
 
+#########################################
+# Now adjust the equilibrium manually
+# 
+
 isoflux = [(Rx,-Zx, Rmid, 0.0)   # Outboard midplane, lower X-point
            ,(Rx,Zx, Rmid, 0.0)   # Outboard midplane, upper X-point
 
@@ -136,6 +121,7 @@
 print("Plasma current: %e Amps" % (eq.plasmaCurrent()))
 print("Pressure on axis: %e Pascals" % (eq.pressure(0.0)))
 print("Poloidal beta: %e" % (eq.poloidalBeta()))
+print("Plasma volume: %e m^3" % (eq.plasmaVolume()))
 
 eq.tokamak.printCurrents()
 

freegs/equilibrium.py

@@ -207,7 +207,42 @@ def poloidalBeta(self):
 
         # Integrate pressure in 2D
         return ((8.*pi)/mu0)*romb(romb(pressure))*dR*dZ / (self.plasmaCurrent()**2)
+
+    def plasmaVolume(self):
+        """Calculate the volume of the plasma in m^3"""
+
+        # Total poloidal flux (plasma + coils)
+        psi = self.psi()
 
+        # Analyse the equilibrium, finding O- and X-points
+        opt, xpt = critical.find_critical(self.R, self.Z, psi)
+        if not opt:
+            raise ValueError("No O-points found!")
+        psi_axis = opt[0][2]
+
+        if xpt:
+            psi_bndry = xpt[0][2]
+            mask = critical.core_mask(self.R, self.Z, psi, opt, xpt)
+        else:
+            # No X-points
+            if psi[0,0] > psi_axis:
+                psi_bndry = np.amax(psi)
+            else:
+                psi_bndry = np.amin(psi)
+            mask = None
+
+        dR = self.R[1,0] - self.R[0,0]
+        dZ = self.Z[0,1] - self.Z[0,0]
+
+        # Volume element
+        dV = 2.*pi*self.R * dR * dZ
+
+        if mask is not None:   # Only include points in the core
+            dV *= mask
+
+        # Integrate volume in 2D
+        return romb(romb(dV))
+
     def plasmaBr(self, R,Z):
         """
         Radial magnetic field due to plasma