Splined bespoke plasma profiles and updated poloidal beta

17-spline-profiles.py

@@ -0,0 +1,141 @@
+#!/usr/bin/env python
+
+import freegs
+import numpy as np
+import matplotlib.pyplot as plt
+
+#########################################
+# Create the machine, which specifies coil locations
+# and equilibrium, specifying the domain to solve over
+
+tokamak = freegs.machine.TestTokamak()
+
+eq = freegs.Equilibrium(tokamak=tokamak,
+                        Rmin=0.1, Rmax=2.0,    # Radial domain
+                        Zmin=-1.0, Zmax=1.0,   # Height range
+                        nx=65, ny=65,          # Number of grid points
+                        boundary=freegs.boundary.freeBoundaryHagenow)  # Boundary condition
+
+
+#########################################
+# Plasma profiles
+
+profiles = freegs.jtor.ConstrainBetapIp(eq,
+                                        3.214806e-02, # Poloidal beta
+                                        2e5, # Plasma current [Amps]
+                                        2.0) # Vacuum f=R*Bt
+
+#########################################
+# Coil current constraints
+#
+# Specify locations of the X-points
+# to use to constrain coil currents
+
+xpoints = [(1.1, -0.6),   # (R,Z) locations of X-points
+           (1.1, 0.8)]
+
+isoflux = [(1.1,-0.6, 1.1,0.6)] # (R1,Z1, R2,Z2) pair of locations
+
+constrain = freegs.control.constrain(xpoints=xpoints, isoflux=isoflux)
+
+#########################################
+# Nonlinear solve
+
+freegs.solve(eq,          # The equilibrium to adjust
+             profiles,    # The toroidal current profile function
+             constrain,
+             show=True)   # Constraint function to set coil currents
+
+# eq now contains the solution
+
+print("Done!")
+
+print("Plasma current: %e Amps" % (eq.plasmaCurrent()))
+print("Plasma pressure on axis: %e Pascals" % (eq.pressure(0.0)))
+print("Poloidal beta: %e" % (eq.poloidalBeta()))
+
+# Now, extract the pprime and ffprime profiles from the solution
+psi_n_data = np.linspace(0.0,1.0,101,endpoint=True)
+pprime_data = eq.pprime(psi_n_data)
+ffprime_data = eq.ffprime(psi_n_data)
+
+# Next, solve again, this time using the above data as a placeholder for bespoke profile data.
+# In user-specific implementations this data may be obtained elsewhere, e.g. from MDSplus or
+# from the output of another code.
+
+print('Using splined profiles for pprime and ffprime')
+
+#########################################
+# Create the machine, which specifies coil locations
+# and equilibrium, specifying the domain to solve over
+
+tokamak2 = freegs.machine.TestTokamak()
+
+eq2 = freegs.Equilibrium(tokamak=tokamak2,
+                        Rmin=0.1, Rmax=2.0,    # Radial domain
+                        Zmin=-1.0, Zmax=1.0,   # Height range
+                        nx=65, ny=65,          # Number of grid points
+                        boundary=freegs.boundary.freeBoundaryHagenow)  # Boundary condition
+
+
+#########################################
+# Plasma profiles
+
+profiles2 = freegs.jtor.BetapIpConstrainedSplineProfiles(eq2,
+                                        3.214806e-02, # Poloidal beta
+                                        2e5, # Plasma current [Amps]
+                                        1.0, # Raxis [m],
+                                        psi_n_data, # Spline data for normalised psi
+                                        pprime_data, # Spline data for pprime
+                                        ffprime_data, # Spline data for ffprime
+                                        2.0) # Vacuum f=R*Bt
+
+#########################################
+# Coil current constraints
+#
+# Specify locations of the X-points
+# to use to constrain coil currents
+
+xpoints = [(1.1, -0.6),   # (R,Z) locations of X-points
+           (1.1, 0.8)]
+
+isoflux = [(1.1,-0.6, 1.1,0.6)] # (R1,Z1, R2,Z2) pair of locations
+
+constrain2 = freegs.control.constrain(xpoints=xpoints, isoflux=isoflux)
+
+#########################################
+# Nonlinear solve
+
+freegs.solve(eq2,          # The equilibrium to adjust
+             profiles2,    # The toroidal current profile function
+             constrain2,
+             show=True)    # Constraint function to set coil currents
+
+# eq2 now contains the solution
+
+print("Done!")
+
+print("Plasma current: %e Amps" % (eq2.plasmaCurrent()))
+print("Plasma pressure on axis: %e Pascals" % (eq2.pressure(0.0)))
+print("Poloidal beta: %e" % (eq2.poloidalBeta()))
+
+# Compare results - parameterised vs splined profiles
+fig, ax = plt.subplots(1,3)
+
+ax[0].contour(eq.R,eq.Z,eq.psi(),levels=[eq.psi_bndry],colors='r')
+ax[0].contour(eq2.R,eq2.Z,eq2.psi(),levels=[eq2.psi_bndry],colors='b')
+ax[0].set_aspect('equal')
+ax[0].set_xlabel('R (m)')
+ax[0].set_ylabel('Z (m)')
+
+ax[1].plot(psi_n_data,eq.pprime(psi_n_data),color='r')
+ax[1].plot(psi_n_data,eq2.pprime(psi_n_data),color='b')
+ax[1].set_xlabel(r'$\psi_{N}$')
+ax[1].set_ylabel(r'pprime($\psi_{N}$)')
+
+ax[2].plot(psi_n_data,eq.ffprime(psi_n_data),color='r')
+ax[2].plot(psi_n_data,eq2.ffprime(psi_n_data),color='b')
+ax[2].set_xlabel(r'$\psi_{N}$')
+ax[2].set_ylabel(r'ffprime($\psi_{N}$)')
+
+plt.show()

freegs/critical.py

@@ -405,7 +405,7 @@ def find_separatrix(
         opoint, xpoint = find_critical(eq.R, eq.Z, psi)
 
     psinorm = (psi - opoint[0][2]) / (eq.psi_bndry - opoint[0][2])
-
+    
     psifunc = interpolate.RectBivariateSpline(eq.R[:, 0], eq.Z[0, :], psinorm)
 
     r0, z0 = opoint[0][0:2]

freegs/equilibrium.py

@@ -10,20 +10,16 @@
 
 from .boundary import fixedBoundary, freeBoundary
 from . import critical
-
 from . import polygons
 
 # Operators which define the G-S equation
 from .gradshafranov import mu0, GSsparse, GSsparse4thOrder
 
 # Multigrid solver
 from . import multigrid
-
 from . import machine
-
 import matplotlib.pyplot as plt
 
-
 class Equilibrium:
     """
     Represents the equilibrium state, including
@@ -126,7 +122,6 @@ def __init__(
 
         self._current = current  # Plasma current
         self.Jtor = None
-
         self._updatePlasmaPsi(psi)  # Needs to be after _pgreen
 
         # Create the solver
@@ -430,6 +425,21 @@ def separatrix(self, npoints=360):
             :, 0:2
         ]
 
+    def psi_surfRZ(self, psiN=0.995, npoints=360):
+        """
+        Returns the R,Z of a flux surface specified by a value of psiN. This flux surface is closed on itself.
+        """
+
+        surf = critical.find_separatrix(self, opoint=None, xpoint=None, ntheta=npoints, psi=None, axis=None, psival=psiN)
+
+        Rsurf = [point[0] for point in surf]
+        Zsurf = [point[1] for point in surf]
+
+        Rsurf.append(Rsurf[0])
+        Zsurf.append(Zsurf[0])
+
+        return np.array(Rsurf), np.array(Zsurf)
+
     def solve(self, profiles, Jtor=None, psi=None, psi_bndry=None):
         """
         Calculate the plasma equilibrium given new profiles
@@ -626,7 +636,7 @@ def _updatePlasmaPsi(self, plasma_psi):
         self.psi_func = interpolate.RectBivariateSpline(
             self.R[:, 0], self.Z[0, :], plasma_psi
         )
-
+        
         # Update the plasma axis and boundary flux as well as mask
         self._updateBoundaryPsi()
 
@@ -1056,32 +1066,46 @@ def internalInductance(self, npoints=360):
         integral = romb(romb(B_polvals_2 * dV))
         return 2 * integral / (mu0 * mu0 * R_geo * Ip * Ip)
 
-    def poloidalBeta(self):
-        """Calculate plasma poloidal beta by integrating the thermal pressure
-        and poloidal magnetic field pressure over the plasma volume."""
+    def flux_surface_averaged_Bpol2(self, psiN=0.995, npoints=360):
+        """
+        Calculates the flux surface averaged value of the square of the poloidal field.
+        """
 
-        R = self.R
-        Z = self.Z
+        # Get R, Z points of the flux surface
+        Rsurf, Zsurf = self.psi_surfRZ(psiN=psiN,npoints=npoints)
 
-        # Produce array of Bpol in (R,Z)
-        B_polvals_2 = self.Br(R, Z) ** 2 + self.Bz(R, Z) ** 2
+        # Get the poloidal field
+        Bpol_surf = self.Bpol(Rsurf,Zsurf)
 
-        dR = R[1, 0] - R[0, 0]
-        dZ = Z[0, 1] - Z[0, 0]
-        dV = 2.0 * np.pi * R * dR * dZ
+        # Get the square of the poloidal field
+        Bpol_surf2 = Bpol_surf**2.0
+
+        # Get dl along the surface
+        dl = np.sqrt(np.diff(Rsurf)**2.0 + np.diff(Zsurf)**2.0)
+        dl = np.insert(dl,0,0.0)
+
+        # Get l along the surface
+        l = np.cumsum(dl)
+
+        # Calculate the flux surface averaged quantity
+        return np.trapz(x=l, y=Bpol_surf2 * Bpol_surf) / np.trapz(x=l, y=np.ones(np.size(l)) * Bpol_surf)
+
+    def poloidalBeta(self):
+        """Return the poloidal beta.
+
+        betaP = 2 * mu0 * <p> / <<Bpol^2>>
+        """
 
         # Normalised psi
         psi_norm = self.psiN()
 
         # Plasma pressure
         pressure = self.pressure(psi_norm)
 
-        if self.mask is not None:  # Only include points in the core
-            dV *= self.mask
+        volume_averaged_pressure = self.calc_volume_averaged(pressure)
+        line_averaged_Bpol2_lcfs = self.flux_surface_averaged_Bpol2(psiN=1.0)
 
-        pressure_integral = romb(romb(pressure * dV))
-        field_integral_pol = romb(romb(B_polvals_2 * dV))
-        return 2 * mu0 * pressure_integral / field_integral_pol
+        return (2.0 * mu0 * volume_averaged_pressure) / line_averaged_Bpol2_lcfs
 
     def poloidalBeta2(self):
         """Return the poloidal beta
@@ -1250,6 +1274,26 @@ def qcyl(self):
 
         return val
 
+    def calc_volume_integrated(self,field):
+        """
+        Calculates the volume integral of the input field.
+        """
+
+        dV = 2.0 * np.pi * self.R * self.dR  *self.dZ
+
+        if self.mask is not None:  # Only include points in the core
+            dV *= self.mask
+
+        return romb(romb(field * dV))
+
+    def calc_volume_averaged(self,field):
+        """
+        Calculates the volume average of the input field.
+        """
+
+        volume_integrated_field = self.calc_volume_integrated(field)
+
+        return volume_integrated_field / self.plasmaVolume()
 
 def refine(eq, nx=None, ny=None):
     """
@@ -1291,7 +1335,6 @@ def refine(eq, nx=None, ny=None):
 
     return result
 
-
 def coarsen(eq):
     """
     Reduce grid resolution, returning a new equilibrium
@@ -1319,7 +1362,6 @@ def coarsen(eq):
 
     return result
 
-
 def newDomain(eq, Rmin=None, Rmax=None, Zmin=None, Zmax=None, nx=None, ny=None):
     """Creates a new Equilibrium, solving in a different domain.
     The domain size (Rmin, Rmax, Zmin, Zmax) and resolution (nx,ny)
@@ -1372,7 +1414,6 @@ def newDomain(eq, Rmin=None, Rmax=None, Zmin=None, Zmax=None, nx=None, ny=None):
 
     return result
 
-
 if __name__ == "__main__":
 
     # Test the different spline interpolation routines