surf ==== .. py:module:: surf Module Contents --------------- .. py:class:: Surface(x0, external_dof_setter, names) ``Surface`` is a base class for various representations of toroidal surfaces. A ``Surface`` is modelled as a function :math:`\Gamma:[0, 1] \times [0, 1] \to R^3` and is evaluated at quadrature points :math:`\{\phi_1, \ldots, \phi_{n_\phi}\}\times\{\theta_1, \ldots, \theta_{n_\theta}\}`. .. py:attribute:: RANGE_FULL_TORUS :value: 'full torus' .. py:attribute:: RANGE_FIELD_PERIOD :value: 'field period' .. py:attribute:: RANGE_HALF_PERIOD :value: 'half period' .. py:method:: from_nphi_ntheta(nphi=61, ntheta=62, range='full torus', nfp=1, **kwargs) :classmethod: Initializes surface classes from the specified number of grid points along toroidal, :math:`\phi`, and poloidal, :math:`\theta`, directions. Additional parameters required for surface initialization could be supplied as keyword arguments. Args: nphi: Number of grid points :math:`\phi_j` in the toroidal angle :math:`\phi`. ntheta: Number of grid points :math:`\theta_i` in the poloidal angle :math:`\theta`. range: Toroidal extent of the :math:`\phi` grid. Set to ``"full torus"`` (or equivalently ``SurfaceRZFourier.RANGE_FULL_TORUS``) to generate quadrature points up to 1 (with no point at 1). Set to ``"field period"`` (or equivalently ``SurfaceRZFourier.RANGE_FIELD_PERIOD``) to generate points up to :math:`1/n_{fp}` (with no point at :math:`1/n_{fp}`). Set to ``"half period"`` (or equivalently ``SurfaceRZFourier.RANGE_HALF_PERIOD``) to generate points up to :math:`1/(2 n_{fp})`, with all grid points shifted by half of the grid spacing in order to provide spectral convergence of integrals. nfp: The number of field periods. kwargs: Additional arguments to initialize the surface classes. Look at the docstrings of the specific class you are interested in. .. py:method:: get_quadpoints(nphi=None, ntheta=None, range=None, nfp=1) :staticmethod: Sets the theta and phi grid points for Surface subclasses. It is typically called in when constructing Surface subclasses. For more information about the arguments ``nphi``, ``ntheta``, ``range``, ``quadpoints_phi``, and ``quadpoints_theta``, see the general documentation on :ref:`surfaces`. Args: nfp: The number of field periods. nphi: Number of grid points :math:`\phi_j` in the toroidal angle :math:`\phi`. ntheta: Number of grid points :math:`\theta_j` in the toroidal angle :math:`\theta`. range: Toroidal extent of the :math:`\phi` grid. Set to ``"full torus"`` (or equivalently ``Surface.RANGE_FULL_TORUS``) to generate points up to 1 (with no point at 1). Set to ``"field period"`` (or equivalently ``Surface.RANGE_FIELD_PERIOD``) to generate points up to :math:`1/n_{fp}` (with no point at :math:`1/n_{fp}`). Set to ``"half period"`` (or equivalently ``Surface.RANGE_HALF_PERIOD``) to generate points up to :math:`1/(2 n_{fp})`, with all grid points shifted by half of the grid spacing in order to provide spectral convergence of integrals. Returns: Tuple containing - **quadpoints_phi**: List of grid points :math:`\phi_j`. - **quadpoints_theta**: List of grid points :math:`\theta_j`. .. py:method:: get_theta_quadpoints(ntheta=None) :staticmethod: Sets the theta grid points for Surface subclasses. Args: ntheta: Number of grid points :math:`\theta_j` in the toroidal angle :math:`\theta`. Returns: List of grid points :math:`\theta_j`. .. py:method:: get_phi_quadpoints(nphi=None, range=None, nfp=1) :staticmethod: Sets the phi grid points for Surface subclasses. Args: nphi: Number of grid points :math:`\phi_j` in the toroidal angle :math:`\phi`. range: Toroidal extent of the :math:`\phi` grid. Set to ``"full torus"`` (or equivalently ``Surface.RANGE_FULL_TORUS``) to generate points up to 1 (with no point at 1). Set to ``"field period"`` (or equivalently ``Surface.RANGE_FIELD_PERIOD``) to generate points up to :math:`1/n_{fp}` (with no point at :math:`1/n_{fp}`). Set to ``"half period"`` (or equivalently ``Surface.RANGE_HALF_PERIOD``) to generate points up to :math:`1/(2 n_{fp})`, with all grid points shifted by half of the grid spacing in order to provide spectral convergence of integrals. nfp: The number of field periods. Returns: List of grid points :math:`\phi_j`. .. py:method:: plot(engine='matplotlib', ax=None, show=True, close=False, axis_equal=True, plot_normal=False, plot_derivative=False, wireframe=True, **kwargs) Plot the surface in 3D using matplotlib/mayavi/plotly. Args: engine: Selects the graphics engine. Currently supported options are ``"matplotlib"`` (default), ``"mayavi"``, and ``"plotly"``. ax: The figure/axis to be plotted on. This argument is useful when plotting multiple objects on the same axes. If equal to the default ``None``, a new axis will be created. show: Whether to call the ``show()`` function of the graphics engine. Should be set to ``False`` if more objects will be plotted on the same axes. close: Whether to close the seams in the surface where the angles jump back to 0. axis_equal: For matplotlib, whether to adjust the scales of the x, y, and z axes so distances in each direction appear equal. plot_normal: Whether to plot the surface normal vectors. Only implemented for mayavi. plot_derivative: Whether to plot the surface derivatives. Only implemented for mayavi. wireframe: Whether to plot the wireframe in Mayavi. kwargs: Any additional arguments to pass to the plotting function, like ``color='r'``. Note: the ``ax`` and ``show`` parameters can be used to plot more than one surface: .. code-block:: python ax = surface1.plot(show=False) ax = surface2.plot(ax=ax, show=False) surface3.plot(ax=ax, show=True) Returns: An axis which could be passed to a further call to the graphics engine so multiple objects are shown together. .. py:method:: to_RZFourier() :abstractmethod: Return a :obj:`SurfaceRZFourier` instance corresponding to the shape of this surface. All subclasses should implement this abstract method. .. py:method:: cross_section(phi, thetas=None, tol=1e-13) Computes an array of points on the cross section with cylindrical angle :math:`\phi`, using the bisection method. The poloidal angles of the points on the cross section are given by ``thetas``. This function assumes that the surface intersection with the :math:`\phi`-plane is a single curve: the surface should only go once around the z-axis, and it should not go back on itself. Args: phi (float): The standard cylindrical angle (toroidal angle) normalized by :math:`2\pi`, i.e. :math:`\phi=0, 1` corresponds to the standard cylindrical angle :math:`0, 2\pi`, respectively. thetas (int, array, optional): An optional argument indicating at which poloidal angle the cross section should be calculated. If ``thetas`` is an ``int``, the cross section will be calculated at the poloidal angles in the array ``np.linspace(0, 1, thetas, endpoint=False)``. ``thetas`` can also be an array of poloidal angles between :math:`0` and :math:`1`. If ``thetas`` is not provided, then the cross section will be calculated at the positions given by ``Surface.quadpoints_theta``. Defaults to ``None``. tol (float): The tolerance for the bisection root-finding. Defaults to ``1e-13``. Returns: ``(ntheta, 3)`` array of the Cartesian coordinates along the cross-section at each of the ``ntheta`` poloidal angles. .. py:method:: cross_section_phi(zetan=0.0, thetas=200) .. py:method:: aspect_ratio() Note: cylindrical coordinates are :math:`(R, \phi, Z)`, where :math:`\phi \in [-\pi,\pi)` and the angles that parametrize the surface are :math:`(\varphi, \theta) \in [0,1)^2` For a given surface, this function computes its aspect ratio using the VMEC definition: .. math:: AR = R_{\text{major}} / R_{\text{minor}} where .. math:: R_{\text{minor}} &= \sqrt{ \overline{A} / \pi } \\ R_{\text{major}} &= \frac{V}{2 \pi^2 R_{\text{minor}}^2} and :math:`V` is the volume enclosed by the surface, and :math:`\overline{A}` is the average cross sectional area. .. py:method:: daspect_ratio_by_dcoeff() Return the derivative of the aspect ratio with respect to the surface coefficients .. py:method:: d2aspect_ratio_by_dcoeff_dcoeff() Return the derivative of the aspect ratio with respect to the surface coefficients .. py:method:: minor_radius() Return the minor radius of the surface using the formula .. math:: R_{\text{minor}} = \sqrt{ \overline{A} / \pi } where :math:`\overline{A}` is the average cross sectional area. .. py:method:: dminor_radius_by_dcoeff() Return the derivative of the minor radius wrt surface coefficients .. py:method:: d2minor_radius_by_dcoeff_dcoeff() Return the second derivative of the minor radius wrt surface coefficients .. py:method:: major_radius() Return the major radius of the surface using the formula .. math:: R_{\text{major}} = \frac{V}{2 \pi^2 R_{\text{minor}}^2} where :math:`\overline{A}` is the average cross sectional area, and :math:`R_{\text{minor}}` is the minor radius of the surface. .. py:method:: dmajor_radius_by_dcoeff() Return the derivative of the major radius wrt surface coefficients .. py:method:: d2major_radius_by_dcoeff_dcoeff() Return the second derivative of the major radius wrt surface coefficients .. py:method:: mean_cross_sectional_area() Note: cylindrical coordinates are :math:`(R, \phi, Z)`, where :math:`\phi \in [-\pi,\pi)` and the angles that parametrize the surface are :math:`(\varphi, \theta) \in [0,1)^2`. The mean cross sectional area is given by the integral .. math:: \overline{A} = \frac{1}{2\pi} \int_{S_{\phi}} ~dS ~d\phi where :math:`S_\phi` is the cross section of the surface at the cylindrical angle :math:`\phi`. Note that :math:`\int_{S_\phi} ~dS` can be rewritten as a line integral .. math:: \int_{S_\phi}~dS &= \int_{S_\phi} ~dR dZ \\ &= \int_{\partial S_\phi} [R,0] \cdot \mathbf n/\|\mathbf n\| ~dl \\ &= \int^1_{0} R \frac{\partial Z}{\partial \theta}~d\theta where :math:`\mathbf n = [n_R, n_Z] = [\partial Z/\partial \theta, -\partial R/\partial \theta]` is the outward pointing normal. Consider the surface in cylindrical coordinates terms of its angles :math:`[R(\varphi,\theta), \phi(\varphi,\theta), Z(\varphi,\theta)]`. The boundary of the cross section :math:`\partial S_\phi` is given by the points :math:`\theta\rightarrow[R(\varphi(\phi,\theta),\theta),\phi, Z(\varphi(\phi,\theta),\theta)]` for fixed :math:`\phi`. The cross sectional area of :math:`S_\phi` becomes .. math:: \int^{1}_{0} R(\varphi(\phi,\theta),\theta) \frac{\partial}{\partial \theta}[Z(\varphi(\phi,\theta),\theta)] ~d\theta Now, substituting this into the formula for the mean cross sectional area, we have .. math:: \overline{A} = \frac{1}{2\pi}\int^{\pi}_{-\pi}\int^{1}_{0} R(\varphi(\phi,\theta),\theta) \frac{\partial}{\partial \theta}[Z(\varphi(\phi,\theta),\theta)] ~d\theta ~d\phi Instead of integrating over cylindrical :math:`\phi`, let's complete the change of variables and integrate over :math:`\varphi` using the mapping: .. math:: [\phi,\theta] \leftarrow [\text{atan2}(y(\varphi,\theta), x(\varphi,\theta)), \theta] After the change of variables, the integral becomes: .. math:: \overline{A} = \frac{1}{2\pi}\int^{1}_{0}\int^{1}_{0} R(\varphi,\theta) \left[\frac{\partial Z}{\partial \varphi} \frac{\partial \varphi}{d \theta} + \frac{\partial Z}{\partial \theta} \right] \text{det} J ~d\theta ~d\varphi where :math:`\text{det}J` is the determinant of the mapping's Jacobian. .. py:method:: dmean_cross_sectional_area_by_dcoeff() Return the derivative of the mean cross sectional area wrt surface coefficients .. py:method:: d2mean_cross_sectional_area_by_dcoeff_dcoeff() Return the second derivative of the mean cross sectional area wrt surface coefficients .. py:method:: arclength_poloidal_angle() Computes a poloidal (angle) coordinate θ on a surface for which the arclength ∂|r|/∂θ is independent of θ in each φ plane. In other words, this function computes the uniform-arclength poloidal coordinate. The returned poloidal coordinate is in the range [0,1), and is used in methods evaluating the adjoint shape gradient. Returns: 2d array of shape ``(numquadpoints_phi, numquadpoints_theta)`` containing the new poloidal angle .. py:method:: interpolate_on_arclength_grid(function, theta_evaluate) Interpolate function onto the theta_evaluate grid in the arclength poloidal angle. This is required for evaluating the adjoint shape gradient for free-boundary calculations. The ``theta_evaluate`` grid may have a different number of poloidal grid points compared to the surface's ``numquadpoints_theta``, but it must have the same number of toroidal grid points as the surface's ``numquadpoints_phi``. This is because we interpolate in theta but not phi. Returns: function_interpolated: 2d array (numquadpoints_phi,numquadpoints_theta) defining interpolated function on arclength angle along curve at constant phi .. py:method:: make_theta_uniform_arclength() Reparameterize the surface in terms of a uniform-arclength poloidal angle. To do the conversion accurately, make sure the surface has both a sufficient number of quadrature points and a sufficiently large number of basis functions. More basis functions may be needed to represent the shape than for the original theta coordinate. .. py:property:: deduced_range The quadpoints of a surface can be anything, but are often set to 'full torus', 'field period' or 'half period'. Since this is not stored in the object, but often useful to know this function deduces the range from the quadpoints