Skip to content

Commit 6e8836d

Browse files
landremanlandreman
committed
First few SPEC examples now work
1 parent 623c506 commit 6e8836d

9 files changed

+598
-8
lines changed

‎examples/.gitignore

+6
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
wout*.nc
2+
mercier*
3+
jxbout*
4+
spec*.sp
5+
spec*.sp.end
6+
spec*.sp.h5

‎examples/1DOF_Garabedian.sp

+83
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
&physicslist
2+
Igeometry = 3
3+
Istellsym = 1
4+
Lfreebound = 0
5+
phiedge = 2.000000000000000E+00
6+
curtor = 1.038123580200000E-09
7+
curpol = 0.000000000000000E+00
8+
gamma = 0.000000000000000E+00
9+
Nfp = 5
10+
Nvol = 1
11+
Mpol = 3
12+
Ntor = 3
13+
Lrad = 8 4
14+
tflux = 1.000000000000000E+00
15+
pflux = -2.040878894181875E-01
16+
helicity = 1.559429589793997E-03
17+
pscale = 0.000000000000000E+00
18+
Ladiabatic = 0
19+
pressure = 0.000000000000000E+00
20+
adiabatic = 1.000000000000000E+00 0.000000000000000E+00
21+
mu = 0.000000000000000E+00
22+
Lconstraint = 0
23+
pl = 0 0 0
24+
ql = 0 0 0
25+
pr = 0 0 0
26+
qr = 0 0 0
27+
iota = 0.000000000000000E+00 2.809417939338480E-01 3.050000000000000E-01
28+
lp = 0 0 0
29+
lq = 0 0 0
30+
rp = 0 0 0
31+
rq = 0 0 0
32+
oita = 0.000000000000000E+00 2.809417939338480E-01 3.050000000000000E-01
33+
mupftol = 1.000000000000000E-12
34+
mupfits = 128
35+
Rac = 1.0 0.1
36+
Zas = 0.0 0.1
37+
RBC(0,0) = 1.0E+00 ZBS(0,0) = 0.0E+00
38+
RBC(1,0) = 1.0E-01 ZBS(1,0) = 1.0E-01
39+
RBC(0,1) = 1.0E-02 ZBS(0,1) = 1.0E-02
40+
/
41+
&numericlist
42+
Linitialize = 0
43+
Ndiscrete = 2
44+
Nquad = -1
45+
iMpol = -4
46+
iNtor = -4
47+
Lsparse = 0
48+
Lsvdiota = 0
49+
imethod = 3
50+
iorder = 2
51+
iprecon = 1
52+
iotatol = -1.000000000000000E+00
53+
/
54+
&locallist
55+
LBeltrami = 4
56+
Linitgues = 1
57+
/
58+
&globallist
59+
Lfindzero = 2
60+
escale = 0.000000000000000E+00
61+
pcondense = 4.000000000000000E+00
62+
forcetol = 1.000000000000000E-12
63+
c05xtol = 1.000000000000000E-12
64+
c05factor = 1.000000000000000E-04
65+
LreadGF = F
66+
opsilon = 1.000000000000000E+00
67+
epsilon = 1.000000000000000E+00
68+
upsilon = 1.000000000000000E+00
69+
/
70+
&diagnosticslist
71+
odetol = 1.000000000000000E-07
72+
absreq = 1.000000000000000E-08
73+
relreq = 1.000000000000000E-08
74+
absacc = 1.000000000000000E-04
75+
epsr = 1.000000000000000E-08
76+
nPpts = 200
77+
nPtrj = 1
78+
LHevalues = F
79+
LHevectors = F
80+
/
81+
&screenlist
82+
Wpp00aa = T
83+
/

‎examples/input.1DOF_Garabedian

+4-4
Original file line numberDiff line numberDiff line change
@@ -42,11 +42,11 @@
4242
PIOTA_TYPE = 'power_series'
4343
PCURR_TYPE = 'power_series'
4444
!----- Axis Parameters -----
45-
RAXIS_CC = 1 1.0e-1
46-
ZAXIS_CS = 0 1.0e-1
45+
RAXIS_CC = 1.0 1.0E-1
46+
ZAXIS_CS = 0.0 1.0E-1
4747
!----- Boundary Parameters -----
4848
! n comes before m!
49-
RBC(0,0) = 1.0E+00 ZBS(0,0) = 0.0000E+00
50-
RBC(1,0) = 1.0E-01 ZBS(1,0) = 1.0e-01
49+
RBC(0,0) = 1.0E+00 ZBS(0,0) = 0.0E+00
50+
RBC(1,0) = 1.0E-01 ZBS(1,0) = 1.0E-01
5151
RBC(0,1) = 1.0E-02 ZBS(0,1) = 1.0E-02
5252
/

‎examples/stellopt_scenarios_1DOF_circularCrossSection_varyAxis_targetIota_spec

+69
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
#!/usr/bin/env python3
2+
3+
import sys
4+
sys.path.append('../src')
5+
from mpi4py import MPI
6+
import numpy as np
7+
import logging
8+
from simsopt import Spec, SurfaceGarabedian, LeastSquaresProblem, \
9+
least_squares_serial_solve
10+
11+
"""
12+
This script implements the "1DOF_circularCrossSection_varyAxis_targetIota"
13+
example from
14+
https://github.com/landreman/stellopt_scenarios
15+
16+
This example demonstrates optimizing a surface shape using the
17+
Garabedian representation instead of VMEC's RBC/ZBS representation.
18+
This optimization problem has one independent variable, the Garabedian
19+
Delta_{m=1, n=-1} coefficient, representing the helical excursion of
20+
the magnetic axis. The objective function is (iota - iota_target)^2,
21+
where iota is measured on the magnetic axis.
22+
23+
Details of the optimum and a plot of the objective function landscape
24+
can be found here:
25+
https://github.com/landreman/stellopt_scenarios/tree/master/1DOF_circularCrossSection_varyAxis_targetIota
26+
"""
27+
28+
logging.basicConfig(level=logging.DEBUG)
29+
30+
# Create a Spec object:
31+
equil = Spec('1DOF_Garabedian.sp')
32+
# If the xspec executable is not in PATH, the path to the executable
33+
# should be specified as in the following line:
34+
# equil = Spec('1DOF_Garabedian.sp', exe='/Users/mattland/SPEC/xspec')
35+
36+
# Start with a default surface, which is axisymmetric with major
37+
# radius 1 and minor radius 0.1.
38+
39+
# We will optimize in the space of Garabedian coefficients rather than
40+
# RBC/ZBS coefficients. To do this, we convert the boundary to the
41+
# Garabedian representation:
42+
surf = equil.boundary.to_Garabedian()
43+
equil.boundary = surf
44+
45+
# VMEC parameters are all fixed by default, while surface parameters
46+
# are all non-fixed by default. You can choose which parameters are
47+
# optimized by setting their 'fixed' attributes.
48+
surf.all_fixed()
49+
surf.set_fixed('Delta(1,-1)', False)
50+
51+
# Each function we want in the objective function is then equipped
52+
# with a shift and weight, to become a term in a least-squares
53+
# objective function. A list of terms are combined to form a
54+
# nonlinear-least-squares problem.
55+
desired_iota = 0.41 # Sign was + for VMEC
56+
prob = LeastSquaresProblem([(equil.iota, desired_iota, 1)])
57+
58+
# Solve the minimization problem. We can choose whether to use a
59+
# derivative-free or derivative-based algorithm.
60+
least_squares_serial_solve(prob, grad=False)
61+
62+
print("At the optimum,")
63+
print(" Delta(m=1,n=-1) = ", surf.get_Delta(1, -1))
64+
print(" iota on axis = ", equil.iota())
65+
print(" objective function = ", prob.objective())
66+
67+
assert np.abs(surf.get_Delta(1, -1) - 0.08575) < 1.0e-4
68+
assert np.abs(equil.iota() - desired_iota) < 1.0e-5
69+
assert prob.objective() < 1.0e-15

‎examples/stellopt_scenarios_1DOF_circularCrossSection_varyR0_targetVolume_spec

+75
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,75 @@
1+
#!/usr/bin/env python3
2+
3+
import logging
4+
import numpy as np
5+
from simsopt import Spec, LeastSquaresProblem, least_squares_serial_solve, SurfaceRZFourier
6+
7+
"""
8+
This script implements the "1DOF_circularCrossSection_varyR0_targetVolume"
9+
example from
10+
https://github.com/landreman/stellopt_scenarios
11+
12+
This optimization problem has one independent variable, representing
13+
the mean major radius. The problem also has one objective: the plasma
14+
volume. There is not actually any need to run an equilibrium code like
15+
VMEC since the objective function can be computed directly from the
16+
boundary shape. But this problem is a fast way to test the
17+
optimization infrastructure with VMEC.
18+
19+
Details of the optimum and a plot of the objective function landscape
20+
can be found here:
21+
https://github.com/landreman/stellopt_scenarios/tree/master/1DOF_circularCrossSection_varyR0_targetVolume
22+
"""
23+
24+
logging.basicConfig(level=logging.DEBUG)
25+
26+
# Create a Spec object:
27+
equil = Spec()
28+
# If the xspec executable is not in PATH, the path to the executable
29+
# should be specified as in the following line:
30+
# equil = Spec(exe='/Users/mattland/SPEC/xspec')
31+
32+
# Start with a default surface, which is axisymmetric with major
33+
# radius 1 and minor radius 0.1.
34+
equil.boundary = SurfaceRZFourier(nfp=5, mpol=5, ntor=4)
35+
surf = equil.boundary
36+
37+
# Set the initial boundary shape. Here is one syntax:
38+
surf.set('rc(0,0)', 1.0)
39+
# Here is another syntax:
40+
surf.set_rc(0, 1, 0.1)
41+
surf.set_zs(0, 1, 0.1)
42+
43+
surf.set_rc(1, 0, 0.1)
44+
surf.set_zs(1, 0, 0.1)
45+
46+
print('rc:', surf.rc)
47+
print('zs:', surf.zs)
48+
#exit(0)
49+
50+
# VMEC parameters are all fixed by default, while surface parameters
51+
# are all non-fixed by default. You can choose which parameters are
52+
# optimized by setting their 'fixed' attributes.
53+
surf.all_fixed()
54+
surf.set_fixed('rc(0,0)', False)
55+
56+
# Each Target is then equipped with a shift and weight, to become a
57+
# term in a least-squares objective function. A list of terms are
58+
# combined to form a nonlinear-least-squares problem.
59+
desired_volume = 0.15
60+
prob = LeastSquaresProblem([(equil.volume, desired_volume, 1)])
61+
62+
# Solve the minimization problem. We can choose whether to use a
63+
# derivative-free or derivative-based algorithm.
64+
least_squares_serial_solve(prob, grad=True)
65+
66+
print("At the optimum,")
67+
print(" rc(m=0,n=0) = ", surf.get_rc(0, 0))
68+
print(" volume, according to SPEC = ", equil.volume())
69+
print(" volume, according to Surface = ", surf.volume())
70+
print(" objective function = ", prob.objective())
71+
72+
assert np.abs(surf.get_rc(0, 0) - 0.7599088773175) < 1.0e-5
73+
assert np.abs(equil.volume() - 0.15) < 1.0e-6
74+
assert np.abs(surf.volume() - 0.15) < 1.0e-6
75+
assert prob.objective() < 1.0e-15

‎src/simsopt/core/surface.py

+2-4
Original file line numberDiff line numberDiff line change
@@ -150,11 +150,9 @@ class SurfaceRZFourier(Surface):
150150
is any poloidal angle.
151151
"""
152152
def __init__(self, nfp=1, stelsym=True, mpol=1, ntor=0):
153+
mpol = int(mpol)
154+
ntor = int(ntor)
153155
# Perform some validation.
154-
if not isinstance(mpol, int):
155-
raise TypeError("mpol must have type int")
156-
if not isinstance(ntor, int):
157-
raise TypeError("ntor must have type int")
158156
if mpol < 1:
159157
raise ValueError("mpol must be at least 1")
160158
if ntor < 0:

‎src/simsopt/mhd/__init__.py

+1
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
from .vmec import vmec_found, Vmec
2+
from .spec import Spec
23

34
#try:
45
# import vmec

0 commit comments

Comments
 (0)