Merge pull request #140 from hiddenSymmetries/graph_surface

Surface and MHD classes converted to graph framework

Author: mbkumar

.github/workflows/extensive_ci.yml

@@ -27,12 +27,12 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: [3.7.10, 3.8.10, 3.9.5] # To sync with coveragerc
         test-type: [unit, integrated]
+        packages: [all, vmec, spec, boost, none]
+        python-version: [3.7.10, 3.8.10, 3.9.5] # To sync with coveragerc
         # spec: [true, false]
         # vmec: [true, false]
         # boost: [true]
-        packages: [none, boost, vmec, spec, all]
         include:
             - python-version: 3.9.5
               test-type: unit

.github/workflows/wheel.yml

@@ -51,8 +51,11 @@ jobs:
         with:
           python-version: ${{ matrix.python }}
 
+      - name: Upgrade pip
+        run: python -m pip install --upgrade pip 
+
       - name: Install wheel builder packages 
-        run: python -m pip install --upgrade pip setuptools wheel delocate
+        run: python -m pip install setuptools wheel delocate
 
       - name: Build and repair wheels
         run: |

.gitignore

@@ -11,6 +11,7 @@ _skbuild/
 build/
 build-cmake/
 dist/
+.eggs/
 pip-wheel-metadata/
 *simsopt.egg-info/
 .DS_Store
@@ -39,9 +40,10 @@ wout*
 mercier*
 jxbout*
 parvmecinfo.txt
-threed*
+*threed*
 vmec_f90wrap.py
 fort.9
+parvmec*.txt
 
 # SPEC specific
 spec*.sp

conda.recipe/conda_build_config.yaml

@@ -1,8 +1,9 @@
 boost:
   - 1.76
 numpy:
+  - 1.19
   - 1.20
 
 pin_run_as_build:
   boost: x.x
-  numpy: x.x
+  numpy: x.x

docs/source/concepts.rst

@@ -203,13 +203,71 @@ A number of quantities are implemented in :obj:`simsopt.geo.curveobjectives` and
 
 The value of the quantity and its derivative with respect to the curve dofs can be obtained by calling e.g., ``CurveLength.J()`` and ``CurveLength.dJ()``.
 
+.. _surfaces:
+
 Surfaces
 ~~~~~~~~
 
-The second large class of geometric objects are surfaces.
-A :obj:`simsopt.geo.surface.Surface` is modelled as a function :math:`\Gamma:[0, 1] \times [0, 1] \to \mathbb{R}^3` and is evaluated at quadrature points :math:`\{\phi_1, \ldots, \phi_{n_\phi}\}\times\{\theta_1, \ldots, \theta_{n_\theta}\}`.
-
-The usage is similar to that of the :obj:`~simsopt.geo.curve.Curve` class:
+The second large class of geometric objects are surfaces.  A
+:obj:`simsopt.geo.surface.Surface` is modelled as a function
+:math:`\Gamma:[0, 1] \times [0, 1] \to \mathbb{R}^3` and is evaluated
+at quadrature points :math:`\{\phi_1, \ldots,
+\phi_{n_\phi}\}\times\{\theta_1, \ldots, \theta_{n_\theta}\}`.  Here,
+:math:`\phi` is the toroidal angle and :math:`\theta` is the poloidal
+angle. Note that :math:`\phi` and :math:`\theta` go up to 1, not up to
+:math:`2 \pi`!
+
+In practice, you almost never use the base
+:obj:`~simsopt.geo.surface.Surface` class.  Rather, you typically use
+one of the subclasses corresponding to a specific parameterization.
+Presently, the available subclasses are
+:obj:`~simsopt.geo.surfacerzfourier.SurfaceRZFourier`,
+:obj:`~simsopt.geo.surfacegarabedian.SurfaceGarabedian`,
+:obj:`~simsopt.geo.surfacehenneberg.SurfaceHenneberg`,
+:obj:`~simsopt.geo.surfacexyzfourier.SurfaceXYZFourier`,
+and
+:obj:`~simsopt.geo.surfacexyztensorfourier.SurfaceXYZTensorFourier`.
+In many cases you can convert a surface from one type to another by going through
+:obj:`~simsopt.geo.surfacerzfourier.SurfaceRZFourier`, as most surface types have
+``to_RZFourier()`` and ``from_RZFourier()`` methods.
+Note that :obj:`~simsopt.geo.surfacerzfourier.SurfaceRZFourier`
+corresponds to the surface parameterization used internally in the VMEC and SPEC codes.
+However when using these codes in simsopt, any of the available surface subclasses
+can be used to represent the surfaces, and simsopt will automatically handle the conversion
+to :obj:`~simsopt.geo.surfacerzfourier.SurfaceRZFourier` when running the code.
+
+The points :math:`\phi_j` and :math:`\theta_j` are used for evaluating
+the position vector and its derivatives, for computing integrals, and
+for plotting, and there are several available methods to specify these
+points.  For :math:`\theta_j`, you typically specify a keyword
+argument ``ntheta`` to the constructor when instantiating a surface
+class. This results in a grid of ``ntheta`` uniformly spaced points
+between 0 and 1, with no endpoint at 1. Alternatively, you can specify
+a list or array of points to the ``quadpoints_theta`` keyword argument
+when instantiating a surface class, specifying the :math:`\theta_j`
+directly.  If both ``ntheta`` and ``quadpoints_theta`` are specified,
+an exception will be raised.  For the :math:`\phi` coordinate, you
+sometimes want points up to 1 (the full torus), sometimes up to
+:math:`1/n_{fp}` (one field period), and sometimes up to :math:`1/(2
+n_{fp})` (half a field period). These three cases can be selected by
+setting the ``range`` keyword argument of the surface subclasses to
+``"full torus"``, ``"field period"``, or ``"half period"``.
+Equivalently, you can set ``range`` to the constants
+``S.RANGE_FULL_TORUS``, ``S.RANGE_FIELD_PERIOD``, or
+``S.RANGE_HALF_PERIOD``, where ``S`` can be
+:obj:`simsopt.geo.surface.Surface` or any of its subclasses.  For all
+of these cases, the ``nphi`` keyword argument can be set to the
+desired number of :math:`\phi` grid points. Alternatively, you can
+pass a list or array to the ``quadpoints_phi`` keyword argument of the
+constructor for any Surface subclass to specify the :math:`\phi_j`
+points directly.  An exception will be raised if both ``nphi`` and
+``quadpoints_phi`` are specified.  For more information about these
+arguments, see the
+:obj:`~simsopt.geo.surfacerzfourier.SurfaceRZFourier` API
+documentation.
+
+The methods available to each surface class are similar to those of
+the :obj:`~simsopt.geo.curve.Curve` class:
 
 - ``Surface.gamma()``: returns a ``(n_phi, n_theta, 3)`` array containing :math:`\Gamma(\phi_i, \theta_j)` for :math:`i\in\{1, \ldots, n_\phi\}, j\in\{1, \ldots, n_\theta\}`, i.e. returns a list of XYZ coordinates on the surface.
 - ``Surface.gammadash1()``: returns a ``(n_phi, n_theta, 3)`` array containing :math:`\partial_\phi \Gamma(\phi_i, \theta_j)` for :math:`i\in\{1, \ldots, n_\phi\}, j\in\{1, \ldots, n_\theta\}`.

docs/source/getting_started.rst

@@ -99,12 +99,16 @@ can be installed for your user only::
 If you want to build SIMSOPT locally with the optional dependencies,
 you can run
 
+.. code-block::
+
     pip install --user -e .[MPI,SPEC]
 
 However, if you're using a zsh terminal (example: latest Macbook versions),
 you'll need to run instead
 
-    pip install --user -e .\[MPI,SPEC]\
+.. code-block::
+
+    pip install --user -e ".[MPI,SPEC]"
 
 
 From docker container

examples/1_Simple/graph_surf_vol_area.py

@@ -1,13 +1,15 @@
 #!/usr/bin/env python3
 
-from simsopt.geo.graph_surface import SurfaceRZFourier
+from simsopt.geo.surfacerzfourier import SurfaceRZFourier
 from simsopt.objectives.graph_least_squares import LeastSquaresProblem
 from simsopt.solve.graph_serial import least_squares_serial_solve
 """
 Optimize the minor radius and elongation of an axisymmetric torus to
 obtain a desired volume and area using the graph method.
 """
 
+print("Running 1_Simple/graph_surf_vol_area.py")
+print("=======================================")
 desired_volume = 0.6
 desired_area = 8.0
 
@@ -20,43 +22,8 @@
 # from the space of independent variables by setting their 'fixed'
 # property to True.
 surf.fix('rc(0,0)')
-surf.get_return_fn_names()
 
 # Approach 1
-
-prob1 = LeastSquaresProblem(depends_on=surf,
-                            goals=[desired_area, desired_volume],
-                            weights=[1, 1])
-least_squares_serial_solve(prob1)
-
-print("At the optimum using approach 1,")
-print(" rc(m=1,n=0) = ", surf.get_rc(1, 0))
-print(" zs(m=1,n=0) = ", surf.get_zs(1, 0))
-print(" volume = ", surf.volume())
-print(" area = ", surf.area())
-print(" objective function = ", prob1.objective())
-print(" -------------------------\n\n")
-
-
-# Approach 2
-
-surf2 = SurfaceRZFourier()
-surf2.fix('rc(0,0)')
-prob2 = LeastSquaresProblem(depends_on=surf2,
-                            opt_return_fns=['area', 'volume'],
-                            goals=[desired_area, desired_volume],
-                            weights=[1, 1])
-least_squares_serial_solve(prob2)
-print("At the optimum using approach 2,")
-print(" rc(m=1,n=0) = ", surf2.get_rc(1, 0))
-print(" zs(m=1,n=0) = ", surf2.get_zs(1, 0))
-print(" volume = ", surf2.volume())
-print(" area = ", surf2.area())
-print(" objective function = ", prob2.objective())
-print(" -------------------------\n\n")
-
-
-# Approach 3
 surf3 = SurfaceRZFourier()
 surf3.fix('rc(0,0)')
 prob3 = LeastSquaresProblem(funcs_in=[surf3.area, surf3.volume],
@@ -71,7 +38,7 @@
 print(" objective function = ", prob3.objective())
 print(" -------------------------\n\n")
 
-# Approach 4
+# Approach 2
 surf4 = SurfaceRZFourier()
 surf4.fix('rc(0,0)')
 prob4 = LeastSquaresProblem.from_tuples([(surf4.area, desired_area, 1),
@@ -86,21 +53,6 @@
 print(" objective function = ", prob4.objective())
 print(" -------------------------\n\n")
 
-# Each target function is then equipped with a shift and weight, to
-# become a term in a least-squares objective function
-#term1 = (surf.volume, desired_volume, 1)
-#term2 = (surf.area,   desired_area,   1)
-
-# A list of terms are combined to form a nonlinear-least-squares
-# problem.
-#prob = LeastSquaresProblem([term1, term2])
-
-# Solve the minimization problem:
-#least_squares_serial_solve(prob)
 
-#print("At the optimum,")
-#print(" rc(m=1,n=0) = ", surf.get_rc(1, 0))
-#print(" zs(m=1,n=0) = ", surf.get_zs(1, 0))
-#print(" volume = ", surf.volume())
-#print(" area = ", surf.area())
-#print(" objective function = ", prob.objective())
+print("End of 1_Simple/graph_surf_vol_area.py")
+print("=======================================")

examples/1_Simple/just_a_quadratic.py

@@ -13,6 +13,8 @@
 # To print out diagnostic information along the way, uncomment this
 # next line:
 #logging.basicConfig(level=logging.INFO)
+print("Running 1_Simple/just_a_quadratic.py")
+print("====================================")
 
 # Define some Target objects that depend on Parameter objects. In the
 # future these functions would involve codes like VMEC, but for now we
@@ -44,3 +46,5 @@
 print("An optimum was found at x=", iden1.x, ", y=", iden2.x, \
       ", z=", iden3.x)
 print("The minimum value of the objective function is ", prob.objective())
+print("End of 1_Simple/just_a_quadratic.py")
+print("====================================")

examples/1_Simple/logger_example.py

@@ -9,6 +9,8 @@
 
 # Serial logging
 initialize_logging(filename='serial.log')
+print("Running 1_Simple/logger_example.py")
+print("==================================")
 for i in range(2):
     logging.info("Hello (times %i) from serial job" % (i+1))
 
@@ -23,3 +25,5 @@
     initialize_logging(mpi=True, filename='mpi.log')
     for i in range(2):
         logging.warning("Hello (times %i) from mpi job" % (i+1))
+print("End of 1_Simple/logger_example.py")
+print("==================================")

examples/1_Simple/minimize_curve_length.py

@@ -10,7 +10,8 @@
 Minimize the length of a curve, holding the 0-frequency Fourier mode fixed.
 The result should be a circle.
 """
-
+print("Running 1_Simple/minimize_curve_length.py")
+print("==============================================")
 # Create a curve:
 nquadrature = 100
 nfourier = 4
@@ -58,3 +59,5 @@
 print(' Final curve length:    ', obj.J())
 print(' Expected final length: ', 2 * np.pi * x0[0])
 print(' objective function: ', prob.objective())
+print("End of 1_Simple/minimize_curve_length.py")
+print("==============================================")