Merge pull request #19 from bendudson/coils

MultiCoil and ShapedCoil types

01-freeboundary.py

@@ -72,6 +72,7 @@
 # Final plot of equilibrium
 
 axis = eq.plot(show=False)
+eq.tokamak.plot(axis=axis, show=False)
 constrain.plot(axis=axis, show=True)
 
 # Safety factor

docs/creating_equilibria.rst

@@ -72,6 +72,49 @@ to fix the current in some of the coils. This can be done by turning off control
 where the current is in Amps, and is for a coil with a single turn. Setting ``control=False``
 removes the coil from feedback control.
 
+
+Shaped coils
+~~~~~~~~~~~~
+
+The ``ShapedCoil`` class models a coil with uniform current density
+across a specified shape (a polygon). This is done by splitting the
+polygon into triangles (using the `ear clipping method
+<https://en.wikipedia.org/wiki/Polygon_triangulation#Ear_clipping_method>`_)
+and then using Gaussian quadrature to integrate over each triangle (by
+default using 6 points per triangle).  The shape of the coil should be
+specified as a list of ``(R,Z)`` points::
+
+  from freegs.machine import ShapedCoil
+  
+  coils = [("P1L", ShapedCoil([(0.95, -1.15), (0.95, -1.05), (1.05, -1.05), (1.05, -1.15)])),
+           ...]
+
+The above would create a square coil. The number of turns and the
+current in the circuit can be specified::
+
+  coils = [("P1L", ShapedCoil([(0.95, -1.15), (0.95, -1.05), (1.05, -1.05), (1.05, -1.15)],
+                              current = 1e3, turns=20)),
+           ...]
+
+The number of turns doesn't affect the integration, but the total
+current in the coil block is set to ``current * turns``.
+
+Multi strand coils
+~~~~~~~~~~~~~~~~~~
+
+The ``MultiCoil`` class provides a way to model a coil block with turns in specified locations.
+For large numbers of turns this is usually more convenient than creating multiple ``Coil`` objects.
+A ``MultiCoil`` can be used anywhere a ``Coil`` object would be::
+  
+  from freegs.machine import MultiCoil
+  
+  coils = [ ("P2", MultiCoil([0.95, 0.95, 1.05, 1.05], [1.15, 1.05, 1.05, 1.15])) ]
+
+For coils which are wired together as pairs, mirrored in Z, the
+``mirror=True`` keyword can be given.  The ``polarity`` keyword then
+sets the relative sign of the current in the original and mirrored
+coil.
+
 Coil circuits
 ~~~~~~~~~~~~~
 

freegs/coil.py

@@ -0,0 +1,255 @@
+"""
+Poloidal field coil
+
+Used in machine to define coils. Can also be a base class for other coil types.
+
+License
+-------
+
+Copyright 2016-2019 Ben Dudson, University of York. Email: benjamin.dudson@york.ac.uk
+
+This file is part of FreeGS.
+
+FreeGS is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+FreeGS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with FreeGS.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from .gradshafranov import Greens, GreensBr, GreensBz, mu0
+import numpy as np
+import numbers
+
+class AreaCurrentLimit:
+    """
+    Calculate the coil area based on a fixed current density limit
+    """
+    def __init__(self, current_density = 3.5e9):
+        """
+        current_density   - Maximum current density in A/m^2
+        
+        Limits in general depend on the magnetic field
+        Typical values Nb3Sn ~ 3.5e9 A/m^2
+        https://doi.org/10.1016/0167-899X(86)90010-8
+        
+        """
+        self._current_density = current_density
+
+    def __call__(self, coil):
+        """
+        Return the area in m^2, given a Coil object
+        """
+        return abs(coil.current * coil.turns) / self._current_density
+
+class Coil:
+    """
+    Represents a poloidal field coil
+
+    public members
+    --------------
+    
+    R, Z - Location of the coil
+    current - current in the coil in Amps
+    turns   - Number of turns
+    control - enable or disable control system
+    area    - Cross-section area in m^2
+
+    The total toroidal current carried by the coil is current * turns
+    """
+
+    # A dtype for converting to Numpy array and storing in HDF5 files
+    dtype = np.dtype([
+        (str("R"), np.float64),
+        (str("Z"), np.float64),
+        (str("current"), np.float64),
+        (str("turns"), np.int),
+        (str("control"), np.bool),
+    ])
+
+    def __init__(self, R, Z, current=0.0, turns=1, control=True, area=AreaCurrentLimit()):
+        """
+        R, Z - Location of the coil
+        
+        current - current in each turn of the coil in Amps
+        turns   - Number of turns. Total coil current is current * turns
+        control - enable or disable control system
+        area    - Cross-section area in m^2
+
+        Area can be a fixed value (e.g. 0.025 for 5x5cm coil), or can be specified
+        using a function which takes a coil as an input argument.
+        To specify a current density limit, use:
+      
+        area = AreaCurrentLimit(current_density)
+        
+        where current_density is in A/m^2. The area of the coil will be recalculated
+        as the coil current is changed.
+
+        The most important effect of the area is on the coil self-force:
+        The smaller the area the larger the hoop force for a given current.
+        """
+        self.R = R
+        self.Z = Z
+
+        self.current = current
+        self.turns = turns
+        self.control = control
+        self.area = area
+
+    def psi(self, R, Z):
+        """
+        Calculate poloidal flux at (R,Z)
+        """
+        return self.controlPsi(R,Z) * self.current
+
+    def createPsiGreens(self, R, Z):
+        """
+        Calculate the Greens function at every point, and return
+        array. This will be passed back to evaluate Psi in
+        calcPsiFromGreens()
+        """
+        return self.controlPsi(R,Z)
+
+    def calcPsiFromGreens(self, pgreen):
+        """
+        Calculate plasma psi from Greens functions and current
+        """
+        return self.current * pgreen
+
+    def Br(self, R, Z):
+        """
+        Calculate radial magnetic field Br at (R,Z)
+        """
+        return self.controlBr(R,Z) * self.current
+
+    def Bz(self, R, Z):
+        """
+        Calculate vertical magnetic field Bz at (R,Z)
+        """
+        return self.controlBz(R,Z) * self.current
+
+    def controlPsi(self, R, Z):
+        """
+        Calculate poloidal flux at (R,Z) due to a unit current
+        """
+        return Greens(self.R, self.Z, R, Z) * self.turns
+
+    def controlBr(self, R, Z):
+        """
+        Calculate radial magnetic field Br at (R,Z) due to a unit current
+        """
+        return GreensBr(self.R,self.Z, R, Z) * self.turns
+
+    def controlBz(self, R, Z):
+        """
+        Calculate vertical magnetic field Bz at (R,Z) due to a unit current
+        """
+        return GreensBz(self.R,self.Z, R, Z) * self.turns
+
+    def getForces(self, equilibrium):
+        """
+        Calculate forces on the coils in Newtons
+        
+        Returns an array of two elements: [ Fr, Fz ]
+
+        
+        Force on coil due to its own current:
+            Lorentz self‐forces on curved current loops
+            Physics of Plasmas 1, 3425 (1998); https://doi.org/10.1063/1.870491
+            David A. Garren and James Chen
+        """
+        current = self.current # current per turn
+        total_current = current * self.turns # Total toroidal current
+
+        # Calculate field at this coil due to all other coils
+        # and plasma. Need to zero this coil's current
+        self.current = 0.0
+        Br = equilibrium.Br(self.R, self.Z)
+        Bz = equilibrium.Bz(self.R, self.Z)
+        self.current = current
+
+        # Assume circular cross-section for hoop (self) force
+        minor_radius = np.sqrt(self.area / np.pi)
+
+        # Self inductance factor, depending on internal current
+        # distribution. 0.5 for uniform current, 0 for surface current
+        self_inductance = 0.5
+
+        # Force per unit length.
+        # In cgs units f = I^2/(c^2 * R) * (ln(8*R/a) - 1 + xi/2)
+        # In SI units f = mu0 * I^2 / (4*pi*R) * (ln(8*R/a) - 1 + xi/2)
+        self_fr = (mu0 * total_current**2 / (4.*np.pi*self.R)) * (np.log(8.*self.R/minor_radius) - 1 + self_inductance/2.)
+
+        Ltor = 2*np.pi*self.R  # Length of coil
+        return np.array([ (total_current * Bz  + self_fr) * Ltor, # Jphi x Bz = Fr, self force always outwards
+                          -total_current * Br * Ltor]) # Jphi x Br = - Fz
+
+    def __repr__(self):
+        return ("Coil(R={0}, Z={1}, current={2:.1f}, turns={3}, control={4})"
+                .format(self.R, self.Z, self.current, self.turns, self.control))
+
+    def __eq__(self, other):
+        return (self.R == other.R
+                and self.Z == other.Z
+                and self.current == other.current
+                and self.turns == other.turns
+                and self.control == other.control)
+
+    def __ne__(self, other):
+        return not self == other
+
+    def to_numpy_array(self):
+        """
+        Helper method for writing output
+        """
+        return np.array((self.R, self.Z, self.current, self.turns, self.control),
+                        dtype=self.dtype)
+
+    @classmethod
+    def from_numpy_array(cls, value):
+        if value.dtype != cls.dtype:
+            raise ValueError("Can't create {this} from dtype: {got} (expected: {dtype})"
+                             .format(this=type(cls), got=value.dtype, dtype=cls.dtype))
+        return Coil(*value[()])
+
+    @property
+    def area(self):
+        """
+        The cross-section area of the coil in m^2
+        """
+        if isinstance(self._area, numbers.Number):
+            assert self._area > 0
+            return self._area
+        # Calculate using functor
+        area = self._area(self)
+        assert area > 0
+        return area
+
+    @area.setter
+    def area(self, area):
+        self._area = area
+
+    def plot(self, axis=None, show=False):
+        """
+        Plot the coil location, using axis if given
+    
+        The area of the coil is used to set the radius
+        """
+        minor_radius = np.sqrt(self.area / np.pi)
+
+        import matplotlib.pyplot as plt
+
+        if axis is None:
+            fig = plt.figure()
+            axis = fig.add_subplot(111)
+
+        circle = plt.Circle((self.R, self.Z), minor_radius, color='b')
+        axis.add_artist(circle)
+        return axis