diff --git a/InputFiles/TestCases/G3V01Lm2Fi.info b/InputFiles/TestCases/G3V01Lm2Fi.info new file mode 100644 index 00000000..defc67a1 --- /dev/null +++ b/InputFiles/TestCases/G3V01Lm2Fi.info @@ -0,0 +1 @@ +test case to verify Lconstraint=-2 (vary toroidal flux to satisfy poloidal linking current constraint) and Lbdybnzero=.false. (non-zero value of Bn on boundary of fixed-boundary calculation). diff --git a/InputFiles/TestCases/G3V01Lm2Fi.sp b/InputFiles/TestCases/G3V01Lm2Fi.sp new file mode 100644 index 00000000..9dcc5d6e --- /dev/null +++ b/InputFiles/TestCases/G3V01Lm2Fi.sp @@ -0,0 +1,566 @@ +! auto generated by SPECNamelist.py 2021-10-21 00:01:17 +&physicslist + igeometry = 3 + istellsym = 1 + lfreebound = 0 + phiedge = 1.0 + curtor = 0 + curpol = 0 + gamma = 0.0 + nfp = 2 + nvol = 1 + mpol = 8 + ntor = 8 + lrad = 10 + tflux = 2.326694073505053 + pflux = 0.0 + helicity = 3.477699658548601 + pscale = 0.0 + ladiabatic = 0 + pressure = 1.0 + adiabatic = 1.0 + mu = 0.0 + lconstraint = -2 + pl = 0, 0 + ql = 0, 0 + pr = 0, 0 + qr = 0, 0 + iota = 0.0, 0.280941793933848 + lp = 0, 0 + lq = 0, 0 + rp = 0, 0 + rq = 0, 0 + oita = 0.0, 0.280941793933848 + mupftol = 1e-16 + mupfits = 16 + rpol = 1.0 + rtor = 1.0 + rbc(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.01299829262699682, + -0.0049241737074355385, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.004916498659397338, 0.12888068732779726, + 0.42319286094669983, -0.021257759852879716, 0.003995604635266758, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0008142021375225027, -0.16136229488354392, + 0.13259437283052, -0.0012247422810209936, 0.001859980407135166, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbc(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.044105780299325124, + 0.004897626703269465, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.008429218722148045, 0.11104611530330084, + -0.6189708354320282, 0.028721518892198167, 0.0012547220152077777, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.012746565240090657, -0.02209424370288384, + -0.04303116340918449, 0.011306578430716792, -0.0003626874352627375, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbs(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rbs(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zbc(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rwc(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zws(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rws(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zwc(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + vns(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0032484357678936454, + 0.04759773531023998, 0.07484164088047297, 0.3491068966898682, + 0.07735720204569335, -0.04213990568445082, 0.0, 0.0 + vns(-8:8,1) = 0.0, 0.0, 0.0514342828823717, 0.24622358712752726, 0.001359929358603049, + -0.025834017904898112, -0.06434987281939539, 0.015416556506977486, + -0.1125546708357967, -0.05658841194313129, -0.03786455364343739, + 0.19124227012881556, -0.3645519514104967, -0.006179595178307682, + 0.057362672292951444, 0.0, 0.0 + vns(-8:8,2) = 0.0, 0.0, -0.14481554819901354, -0.42168301273828834, + 0.16233979647432484, 0.31330537429640026, 0.10321400286640497, + 0.05709562946522295, -0.1371688175533019, -0.007672840714461098, + -0.08546794825843257, 0.17723799697567313, 0.15340356518419954, + -0.035080425027089725, -0.02171370713627696, 0.0, 0.0 + vns(-8:8,3) = 0.0, 0.0, 0.15655040702185627, 0.4174834890894755, -0.2642683215937827, + 0.02223199880007334, 0.021430317480327095, 0.07416427163336486, + -0.19931478721188262, -0.06519672896461476, -0.03358927168429972, + -0.26428948720531764, -0.04440912116686778, 0.03456978237671945, + 0.013324200166517247, 0.0, 0.0 + vns(-8:8,4) = 0.0, 0.0, -0.08329033006640713, -0.25694437993929026, + 0.08212352196213508, -0.018307879744856273, -0.01828899126210851, + -0.03882713915994556, -0.0736253412037836, 0.0265586650645998, + 0.1499083616226613, 0.03881004406830072, -0.030764770976077865, + 0.0064636135049393455, 0.006005704806467853, 0.0, 0.0 + vns(-8:8,5) = 0.0, 0.0, 0.008434770759069172, 0.14359133925639994, 0.310005408232767, + -0.5121129424372788, 0.01576343325071893, 0.12870282106583014, + 0.037174669913650815, 0.0007705639677875958, -0.035265435500200705, + 0.01494259944518275, 0.04833979355061212, -0.006121826938287074, + -0.016188729223063716, 0.0, 0.0 + vns(-8:8,6) = 0.0, 0.0, -0.005842569954613564, -0.08871007887173403, + -0.47788806490147845, 0.25892135117550424, 0.045829079344694476, + 0.22170271928799812, 0.11974295868370594, 0.05897964449697521, + -0.015329979047629092, 0.020342062858547474, 0.007038146692507075, + -0.01176393507816326, -0.006930676016976327, 0.0, 0.0 + vns(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + vns(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + vnc(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,1) = 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,2) = 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,3) = 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,4) = 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,5) = 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,6) = 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, + -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0 + vnc(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + vnc(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bns(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,0) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,1) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,2) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,3) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,4) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,5) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,6) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,7) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + bnc(-8:8,8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + rac(0:8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zas(0:8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + ras(0:8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + zac(0:8) = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + lbdybnzero = .false. +/ + +&numericlist + linitialize = 0 + lautoinitbn = 0 + lzerovac = 0 + ndiscrete = 2 + nquad = -1 + impol = -4 + intor = -4 + lsparse = 0 + lsvdiota = 0 + imethod = 3 + iorder = 2 + iprecon = 1 + iotatol = -1.0 + lextrap = 0 + mregular = -1 + lrzaxis = 2 +/ + +&locallist + lbeltrami = 4 + linitgues = 1 +/ + +&globallist + lfindzero = 2 + escale = 0.0 + opsilon = 1.0 + pcondense = 4.0 + epsilon = 1.0 + wpoloidal = 1.0 + upsilon = 1.0 + forcetol = 1e-14 + c05xmax = 1e-06 + c05xtol = 1e-14 + c05factor = 0.0001 + lreadgf = .true. + mfreeits = 0 + gbntol = 1e-06 + gbnbld = 0.666 + vcasingeps = 1e-12 + vcasingtol = 1e-08 + vcasingits = 8 + vcasingper = 1 +/ + +&diagnosticslist + odetol = 1e-07 + nppts = 0 + nptrj = 8 + lhevalues = .false. + lhevectors = .false. + lhmatrix = .false. + lperturbed = 0 + dpp = -1 + dqq = -1 + lcheck = 0 + ltiming = .false. +/ + +&screenlist +/ + 0 0 5.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 -1 5.000000000000000e-01 5.000000000000000e-01 0.000000000000000e+00 0.000000000000000e+00 + 1 0 1.500000000000000e+00 -1.500000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 4 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 6 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 7 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 8 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 9 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 -1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 0 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 10 10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 diff --git a/Makefile b/Makefile index 1b0ba938..12947884 100644 --- a/Makefile +++ b/Makefile @@ -47,7 +47,6 @@ inputlist_d.o: %_d.o: src/inputlist.f90 $(MACROS) @echo '' ############################################################################################################################################################### -# global needs special handling: expansion of CPUVARIABLE, BSCREENLIST and WSCREENLIST using awk (not anymore !!!) global_r.o: %_r.o: inputlist_r.o src/global.f90 $(MACROS) m4 -P $(MACROS) src/global.f90 > global_m.F90 diff --git a/Utilities/pythontools/py_spec/input/spec_namelist.py b/Utilities/pythontools/py_spec/input/spec_namelist.py index 344113cc..b9a9ec9f 100644 --- a/Utilities/pythontools/py_spec/input/spec_namelist.py +++ b/Utilities/pythontools/py_spec/input/spec_namelist.py @@ -216,7 +216,10 @@ def run(self, spec_command="./xspec", filename="spec.sp", force=False, quiet=Fal if run_result.returncode == 0: # the run is successful if not quiet: print("SPEC runs successfully.") - return SPECout(filename + ".h5") + try: + return SPECout(filename + ".h5") + except: + return SPECout(filename[:-3] + ".h5") else: print("SPEC runs unsuccessfully, check terminal output.") return None diff --git a/Utilities/pythontools/py_spec/output/_processing.py b/Utilities/pythontools/py_spec/output/_processing.py index 3247f6d0..277c83b5 100644 --- a/Utilities/pythontools/py_spec/output/_processing.py +++ b/Utilities/pythontools/py_spec/output/_processing.py @@ -1,5 +1,5 @@ import numpy as np - +import gc def get_grid_and_jacobian_and_metric( self, @@ -166,21 +166,235 @@ def metric( return g +def get_B_FFT( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + Nt=1, + Nz=1 +): + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr + ) + + # Lrad = s.input.physics.Lrad[lvol] + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] + + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + Nfp = self.input.physics.Nfp + Ns = len(sarr) + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + + # Zernike polynomial being used + from ._processing import _get_zernike + zernike, dzernike, _ = _get_zernike(sbar, Lrad, Mpol) + + # Obtain Bs + Bs_mn_cos = np.sum( zernike[:, :, im] * (im[nax,:]*Azo.T[nax, :, :] + in_[nax,:]*Ato.T[nax, :, :]) , axis=(1)) + Bs_mn_sin =-np.sum( zernike[:, :, im] * (im[nax,:]*Aze.T[nax, :, :] + in_[nax,:]*Ate.T[nax, :, :]) , axis=(1)) + Bs = invfft_B(Bs_mn_cos, Bs_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain Bt + Bt_mn_cos = -np.sum( dzernike[:, :, im] * Aze.T[nax, :, :], axis=(1)) + Bt_mn_sin = -np.sum( dzernike[:, :, im] * Azo.T[nax, :, :], axis=(1)) + Bt = invfft_B(Bt_mn_cos, Bt_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain Bz + Bz_mn_cos = np.sum( dzernike[:, :, im] * Ate.T[nax, :, :], axis=(1)) + Bz_mn_sin = np.sum( dzernike[:, :, im] * Ato.T[nax, :, :], axis=(1)) + Bz = invfft_B(Bz_mn_cos, Bz_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + return np.array([Bs, Bt, Bz]) / jacobian + +def get_s_der_B_FFT( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + Nt=1, + Nz=1 +): + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr + ) + + # Lrad = s.input.physics.Lrad[lvol] + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] + + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + Nfp = self.input.physics.Nfp + Ns = len(sarr) + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + + # Zernike polynomial being used + from ._processing import _get_zernike + _, dzernike, d2zernike = _get_zernike(sbar, Lrad, Mpol) + + # Obtain dsBs + dsBs_mn_cos = np.sum( dzernike[:, :, im] * (im[nax,:]*Azo.T[nax, :, :] + in_[nax,:]*Ato.T[nax, :, :]) , axis=(1)) + dsBs_mn_sin =-np.sum( dzernike[:, :, im] * (im[nax,:]*Aze.T[nax, :, :] + in_[nax,:]*Ate.T[nax, :, :]) , axis=(1)) + dsBs = invfft_B(dsBs_mn_cos, dsBs_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain dsBt + dsBt_mn_cos = -np.sum( d2zernike[:, :, im] * Aze.T[nax, :, :], axis=(1)) + dsBt_mn_sin = -np.sum( d2zernike[:, :, im] * Azo.T[nax, :, :], axis=(1)) + dsBt = invfft_B(dsBt_mn_cos, dsBt_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain dsBz + dsBz_mn_cos = np.sum( d2zernike[:, :, im] * Ate.T[nax, :, :], axis=(1)) + dsBz_mn_sin = np.sum( d2zernike[:, :, im] * Ato.T[nax, :, :], axis=(1)) + dsBz = invfft_B(dsBz_mn_cos, dsBz_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + return np.array([dsBs, dsBt, dsBz]) / jacobian + +def get_t_der_B_FFT( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + Nt=1, + Nz=1 +): + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr + ) + + # Lrad = s.input.physics.Lrad[lvol] + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] + + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + Nfp = self.input.physics.Nfp + Ns = len(sarr) + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + + # Zernike polynomial being used + from ._processing import _get_zernike + zernike, dzernike, _ = _get_zernike(sbar, Lrad, Mpol) + + # Obtain dtBs + dtBs_mn_sin = -np.sum( zernike[:, :, im] * im[nax,:] * (im[nax,:]*Azo.T[nax, :, :] + in_[nax,:]*Ato.T[nax, :, :]) , axis=(1)) + dtBs_mn_cos = -np.sum( zernike[:, :, im] * im[nax,:] * (im[nax,:]*Aze.T[nax, :, :] + in_[nax,:]*Ate.T[nax, :, :]) , axis=(1)) + dtBs = invfft_B(dtBs_mn_cos, dtBs_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain dtBt + dtBt_mn_sin = np.sum( dzernike[:, :, im] * im[nax,:] * Aze.T[nax, :, :], axis=(1)) + dtBt_mn_cos = -np.sum( dzernike[:, :, im] * im[nax,:] * Azo.T[nax, :, :], axis=(1)) + dtBt = invfft_B(dtBt_mn_cos, dtBt_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain dtBz + dtBz_mn_sin = -np.sum( dzernike[:, :, im] * im[nax,:] * Ate.T[nax, :, :], axis=(1)) + dtBz_mn_cos = np.sum( dzernike[:, :, im] * im[nax,:] * Ato.T[nax, :, :], axis=(1)) + dtBz = invfft_B(dtBz_mn_cos, dtBz_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + return np.array([dtBs, dtBt, dtBz]) / jacobian + + +def get_z_der_B_FFT( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + Nt=1, + Nz=1 +): + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr + ) + + # Lrad = s.input.physics.Lrad[lvol] + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] + + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + Nfp = self.input.physics.Nfp + Ns = len(sarr) + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + + # Zernike polynomial being used + from ._processing import _get_zernike + zernike, dzernike, _ = _get_zernike(sbar, Lrad, Mpol) + + # Obtain dzBs + dzBs_mn_sin = np.sum( zernike[:, :, im] * in_[nax,:] * (im[nax,:]*Azo.T[nax, :, :] + in_[nax,:]*Ato.T[nax, :, :]) , axis=(1)) + dzBs_mn_cos = np.sum( zernike[:, :, im] * in_[nax,:] * (im[nax,:]*Aze.T[nax, :, :] + in_[nax,:]*Ate.T[nax, :, :]) , axis=(1)) + dzBs = invfft_B(dzBs_mn_cos, dzBs_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain dzBt + dzBt_mn_sin = -np.sum( dzernike[:, :, im] * in_[nax,:] * Aze.T[nax, :, :], axis=(1)) + dzBt_mn_cos = np.sum( dzernike[:, :, im] * in_[nax,:] * Azo.T[nax, :, :], axis=(1)) + dzBt = invfft_B(dzBt_mn_cos, dzBt_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + # Obtain dzBz + dzBz_mn_sin = np.sum( dzernike[:, :, im] * in_[nax,:] * Ate.T[nax, :, :], axis=(1)) + dzBz_mn_cos = -np.sum( dzernike[:, :, im] * in_[nax,:] * Ato.T[nax, :, :], axis=(1)) + dzBz = invfft_B(dzBz_mn_cos, dzBz_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz) + + return np.array([dzBs, dzBt, dzBz]) / jacobian + def get_B( self, lvol=0, jacobian=None, - sarr=np.linspace(0, 0, 1), + sarr=np.linspace(1, 1, 1), tarr=np.linspace(0, 0, 1), zarr=np.linspace(0, 0, 1), ): - if jacobian is None: R, Z, jacobian, g = get_grid_and_jacobian_and_metric( self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr ) # Lrad = s.input.physics.Lrad[lvol] + + Ate = self.vector_potential.Ate[lvol] Aze = self.vector_potential.Aze[lvol] Ato = self.vector_potential.Ato[lvol] @@ -206,40 +420,110 @@ def get_B( if LZernike: # Zernike polynomial being used from ._processing import _get_zernike + zernike, dzernike, _ = _get_zernike(sbar, Lrad, Mpol) + else: + # Chebyshev being used + import numpy.polynomial.chebyshev as Cheb + # make basis recombination for cheb + lcoeff = np.arange(0, Lrad + 1) + 1 + + Ate = Ate / lcoeff[None, :] + Aze = Aze / lcoeff[None, :] + Ato = Ato / lcoeff[None, :] + Azo = Azo / lcoeff[None, :] - zernike, dzernike = _get_zernike(sarr, Lrad, Mpol) + Ate[:, 0] = np.sum(Ate * (-1.0) ** lcoeff[None, :], 1) + Aze[:, 0] = np.sum(Aze * (-1.0) ** lcoeff[None, :], 1) + Ato[:, 0] = np.sum(Ato * (-1.0) ** lcoeff[None, :], 1) + Azo[:, 0] = np.sum(Azo * (-1.0) ** lcoeff[None, :], 1) - c = ( - im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] - + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] - ) * cosa[nax, :, :, :] - ( - im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] - + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] - ) * sina[ - nax, :, :, : - ] - Bs = np.sum(zernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) + # Obtain Bs + c = ( + im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] + ) * cosa[nax, :, :, :] - ( + im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] + ) * sina[nax, :, :, :] - c1 = ( - Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] - + Azo.T[:, :, nax, nax] * sina[nax, :, :, :] - ) - Bt = -np.sum(dzernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) + if LZernike: + Bs = np.sum( zernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) + else: + Bs = np.rollaxis(np.sum(Cheb.chebval(sarr, c), axis=0), 2) + + # Obtain Bt + c1 = ( + Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] + + Azo.T[:, :, nax, nax] * sina[nax, :, :, :]) + + if LZernike: + Bt = -np.sum( dzernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) + else: + Bt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1)), axis=0),2) + + # Obtain Bz + c2 = ( + Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] + + Ato.T[:, :, nax, nax] * sina[nax, :, :, :]) + + if LZernike: + Bz = np.sum( dzernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) + else: + Bz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2)), axis=0),2) + + return np.array([Bs, Bt, Bz]) / jacobian + + +def get_s_der_B( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + tarr=np.linspace(0, 0, 1), + zarr=np.linspace(0, 0, 1), +): - c2 = ( - Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] - + Ato.T[:, :, nax, nax] * sina[nax, :, :, :] + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr ) - Bz = np.sum(dzernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) + # Lrad = s.input.physics.Lrad[lvol] + + + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] + + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + # [mn,it,iz] + ang_arg = im[:, nax, nax] * tarr[nax, :, nax] - in_[:, nax, nax] * zarr[nax, nax, :] + cosa = np.cos(ang_arg) + sina = np.sin(ang_arg) + + LZernike = self.input.physics.Igeometry > 1 and lvol == 0 + if LZernike: + # Zernike polynomial being used + from ._processing import _get_zernike + _, dzernike, d2zernike = _get_zernike(sbar, Lrad, Mpol) else: # Chebyshev being used import numpy.polynomial.chebyshev as Cheb - - lcoeff = np.arange(0, Lrad + 1) + 1 # make basis recombination for cheb + lcoeff = np.arange(0, Lrad + 1) + 1 + Ate = Ate / lcoeff[None, :] Aze = Aze / lcoeff[None, :] Ato = Ato / lcoeff[None, :] @@ -250,48 +534,275 @@ def get_B( Ato[:, 0] = np.sum(Ato * (-1.0) ** lcoeff[None, :], 1) Azo[:, 0] = np.sum(Azo * (-1.0) ** lcoeff[None, :], 1) - # Ch ,mn,t ,z - c = ( - im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] - + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] - ) * cosa[nax, :, :, :] - ( - im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] - + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] - ) * sina[ - nax, :, :, : - ] - - Bs = np.rollaxis(np.sum(Cheb.chebval(sarr, c), axis=0), 2) - - c1 = ( - Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] - + Azo.T[:, :, nax, nax] * sina[nax, :, :, :] - ) - Bt = -np.rollaxis( - np.sum( - Cheb.chebval(sarr, Cheb.chebder(c1)), - axis=0, - ), - 2, - ) + # Obtain dsBs + c = ( + im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] + ) * cosa[nax, :, :, :] - ( + im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] + ) * sina[nax, :, :, :] + + if LZernike: + dsBs = np.sum(dzernike[:, :, im, None, None] * c[None, :, :, :, :], axis=(1, 2)) + else: + dsBs = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c)), axis=0), 2) + + # Obtain dsBt + c1 = ( + Aze.T[:, :, nax, nax] * cosa[nax, :, :, :] + + Azo.T[:, :, nax, nax] * sina[nax, :, :, :]) + + if LZernike: + dsBt = -np.sum(d2zernike[:, :, im, None, None] * c1[None, :, :, :, :], axis=(1, 2)) + else: + dsBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1,m=2)), axis=0),2) + + # Obtain dsBz + c2 = ( + Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] + + Ato.T[:, :, nax, nax] * sina[nax, :, :, :]) + + if LZernike: + dsBz = np.sum(d2zernike[:, :, im, None, None] * c2[None, :, :, :, :], axis=(1, 2)) + else: + dsBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2,m=2)), axis=0),2) - c2 = ( - Ate.T[:, :, nax, nax] * cosa[nax, :, :, :] - + Ato.T[:, :, nax, nax] * sina[nax, :, :, :] + return np.array([dsBs, dsBt, dsBz]) / jacobian + + +def get_t_der_B( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + tarr=np.linspace(0, 0, 1), + zarr=np.linspace(0, 0, 1), +): + + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr ) - Bz = np.rollaxis( - np.sum( - Cheb.chebval(sarr, Cheb.chebder(c2)), - axis=0, - ), - 2, + + # Lrad = s.input.physics.Lrad[lvol] + + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] + + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + # [mn,it,iz] + ang_arg = im[:, nax, nax] * tarr[nax, :, nax] - in_[:, nax, nax] * zarr[nax, nax, :] + cosa = np.cos(ang_arg) + sina = np.sin(ang_arg) + + LZernike = self.input.physics.Igeometry > 1 and lvol == 0 + + if LZernike: + # Zernike polynomial being used + from ._processing import _get_zernike + zernike, dzernike, _ = _get_zernike(sbar, Lrad, Mpol) + else: + # Chebyshev being used + import numpy.polynomial.chebyshev as Cheb + # make basis recombination for cheb + lcoeff = np.arange(0, Lrad + 1) + 1 + + Ate = Ate / lcoeff[None, :] + Aze = Aze / lcoeff[None, :] + Ato = Ato / lcoeff[None, :] + Azo = Azo / lcoeff[None, :] + + Ate[:, 0] = np.sum(Ate * (-1.0) ** lcoeff[None, :], 1) + Aze[:, 0] = np.sum(Aze * (-1.0) ** lcoeff[None, :], 1) + Ato[:, 0] = np.sum(Ato * (-1.0) ** lcoeff[None, :], 1) + Azo[:, 0] = np.sum(Azo * (-1.0) ** lcoeff[None, :], 1) + + # Obtain dtBs + ct = im[nax, :, nax, nax] * ( - ( + im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] + ) * sina[nax, :, :, :] - ( + im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] + ) * cosa[nax, :, :, :] ) + + if LZernike: + dtBs = np.sum( zernike[:, :, im, None, None] * ct[None, :, :, :, :], axis=(1, 2)) + else: + dtBs = np.rollaxis(np.sum(Cheb.chebval(sarr, ct), axis=0), 2) + + # Obtain dtBt + c1t = im[nax, :, nax, nax] * ( + - Aze.T[:, :, nax, nax] * sina[nax, :, :, :] + + Azo.T[:, :, nax, nax] * cosa[nax, :, :, :]) + + if LZernike: + dtBt = -np.sum( dzernike[:, :, im, None, None] * c1t[None, :, :, :, :], axis=(1, 2)) + else: + dtBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1t)), axis=0),2) + + # Obtain dtBz + c2t = im[nax, :, nax, nax] * ( + - Ate.T[:, :, nax, nax] * sina[nax, :, :, :] + + Ato.T[:, :, nax, nax] * cosa[nax, :, :, :]) + + if LZernike: + dtBz = np.sum( dzernike[:, :, im, None, None] * c2t[None, :, :, :, :], axis=(1, 2)) + else: + dtBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2t)), axis=0),2) + + return np.array([dtBs, dtBt, dtBz]) / jacobian + + +def get_z_der_B( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + tarr=np.linspace(0, 0, 1), + zarr=np.linspace(0, 0, 1), +): + + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr ) - Bcontrav = np.array([Bs, Bt, Bz]) / jacobian - return Bcontrav + # Lrad = s.input.physics.Lrad[lvol] + + + Ate = self.vector_potential.Ate[lvol] + Aze = self.vector_potential.Aze[lvol] + Ato = self.vector_potential.Ato[lvol] + Azo = self.vector_potential.Azo[lvol] + + mn = Ate.shape[0] + im = self.output.im + in_ = self.output.in_ + Lrad = self.input.physics.Lrad[lvol] + Mpol = self.input.physics.Mpol + + fac = [] + sbar = (sarr + 1) / 2 + + nax = np.newaxis + # [mn,it,iz] + ang_arg = im[:, nax, nax] * tarr[nax, :, nax] - in_[:, nax, nax] * zarr[nax, nax, :] + cosa = np.cos(ang_arg) + sina = np.sin(ang_arg) + + LZernike = self.input.physics.Igeometry > 1 and lvol == 0 + + if LZernike: + # Zernike polynomial being used + from ._processing import _get_zernike + zernike, dzernike, _ = _get_zernike(sbar, Lrad, Mpol) + else: + # Chebyshev being used + import numpy.polynomial.chebyshev as Cheb + # make basis recombination for cheb + lcoeff = np.arange(0, Lrad + 1) + 1 + + Ate = Ate / lcoeff[None, :] + Aze = Aze / lcoeff[None, :] + Ato = Ato / lcoeff[None, :] + Azo = Azo / lcoeff[None, :] + + Ate[:, 0] = np.sum(Ate * (-1.0) ** lcoeff[None, :], 1) + Aze[:, 0] = np.sum(Aze * (-1.0) ** lcoeff[None, :], 1) + Ato[:, 0] = np.sum(Ato * (-1.0) ** lcoeff[None, :], 1) + Azo[:, 0] = np.sum(Azo * (-1.0) ** lcoeff[None, :], 1) + + # Obtain dzBs + cz = in_[nax, :, nax, nax] * ( ( + im[nax, :, nax, nax] * Azo.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ato.T[:, :, nax, nax] + ) * sina[nax, :, :, :] + ( + im[nax, :, nax, nax] * Aze.T[:, :, nax, nax] + + in_[nax, :, nax, nax] * Ate.T[:, :, nax, nax] + ) * cosa[nax, :, :, :] ) + + if LZernike: + dzBs = np.sum( zernike[:, :, im, None, None] * cz[None, :, :, :, :], axis=(1, 2)) + else: + dzBs = np.rollaxis(np.sum(Cheb.chebval(sarr, cz), axis=0), 2) + + # Obtain dzBt + c1z = in_[nax, :, nax, nax] * ( + + Aze.T[:, :, nax, nax] * sina[nax, :, :, :] + - Azo.T[:, :, nax, nax] * cosa[nax, :, :, :]) + + if LZernike: + dzBt = -np.sum( dzernike[:, :, im, None, None] * c1z[None, :, :, :, :], axis=(1, 2)) + else: + dzBt = -np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c1z)), axis=0),2) + + # Obtain dzBz + c2z = in_[nax, :, nax, nax] * ( + + Ate.T[:, :, nax, nax] * sina[nax, :, :, :] + - Ato.T[:, :, nax, nax] * cosa[nax, :, :, :]) + + if LZernike: + dzBz = np.sum( dzernike[:, :, im, None, None] * c2z[None, :, :, :, :], axis=(1, 2)) + else: + dzBz = np.rollaxis(np.sum(Cheb.chebval(sarr, Cheb.chebder(c2z)), axis=0),2) + + return np.array([dzBs, dzBt, dzBz]) / jacobian + + +def get_B_and_der_FFT( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + Nt=1, + Nz=1, +): + + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr) + + Bs, Bt, Bz = get_B_FFT( self, lvol, jacobian, sarr, Nt, Nz) + dsBs, dsBt, dsBz = get_s_der_B_FFT(self, lvol, jacobian, sarr, Nt, Nz) + dtBs, dtBt, dtBz = get_t_der_B_FFT(self, lvol, jacobian, sarr, Nt, Nz) + dzBs, dzBt, dzBz = get_z_der_B_FFT(self, lvol, jacobian, sarr, Nt, Nz) + + Bcontrav_and_der = np.array([Bs, Bt, Bz, dsBs, dsBt, dsBz, dtBs, dtBt, dtBz, dzBs, dzBt, dzBz]) + return Bcontrav_and_der + +def get_B_and_der( + self, + lvol=0, + jacobian=None, + sarr=np.linspace(1, 1, 1), + tarr=np.linspace(0, 0, 1), + zarr=np.linspace(0, 0, 1), +): + + if jacobian is None: + R, Z, jacobian, g = get_grid_and_jacobian_and_metric( + self, lvol=lvol, sarr=sarr, tarr=tarr, zarr=zarr) + Bs, Bt, Bz = get_B( self, lvol, jacobian, sarr, tarr, zarr) + dsBs, dsBt, dsBz = get_s_der_B(self, lvol, jacobian, sarr, tarr, zarr) + dtBs, dtBt, dtBz = get_t_der_B(self, lvol, jacobian, sarr, tarr, zarr) + dzBs, dzBt, dzBz = get_z_der_B(self, lvol, jacobian, sarr, tarr, zarr) -# Bcontrav = get_B(s,lvol=lvol,jacobian=jacobian,sarr=sarr,tarr=tarr,zarr=zarr) + Bcontrav_and_der = np.array([Bs, Bt, Bz, dsBs, dsBt, dsBz, dtBs, dtBt, dtBz, dzBs, dzBt, dzBz]) + return Bcontrav_and_der def get_modB(self, Bcontrav, g): @@ -307,7 +818,7 @@ def get_B_covariant(self, Bcontrav, g): def _get_zernike(sarr, lrad, mpol): """ - Get the value of the zernike polynomials and their derivatives + Get the value of the zernike polynomials, their first and second derivatives Adapted from basefn.f90 """ @@ -316,27 +827,27 @@ def _get_zernike(sarr, lrad, mpol): r = (sarr + 1.0) / 2 rm = np.ones_like(r) # r to the power of m'th rm1 = np.zeros_like(r) # r to the power of m-1'th + rm2 = np.zeros_like(r) # r to the power of m-2'th - zernike = np.zeros([ns, lrad + 1, mpol + 1], dtype=np.float64) + zernike = np.zeros(shape=(ns, lrad + 1, mpol + 1), dtype=np.float64) dzernike = np.zeros_like(zernike) + d2zernike = np.zeros_like(zernike) for m in range(mpol + 1): if lrad >= m: - zernike[:, m, m] = rm - dzernike[:, m, m] = m * rm1 + zernike[:, m, m] = rm + dzernike[:, m, m] = m * rm1 + d2zernike[:, m, m] = m * (m-1) * rm2 if lrad >= m + 2: - zernike[:, m + 2, m] = float(m + 2) * rm * r ** 2 - float(m + 1) * rm - dzernike[:, m + 2, m] = ( - float(m + 2) ** 2 * rm * r - float((m + 1) * m) * rm1 - ) + zernike[:, m + 2, m] = float(m + 2) * rm * r**2 - float(m + 1) * rm + dzernike[:, m + 2, m] = (float(m + 2)**2 * rm * r - float((m + 1) * m) * rm1) + d2zernike[:, m + 2, m] = (float((m + 2)**2 * (m + 1)) * rm - float((m + 1) * m * (m - 1)) * rm2) for n in range(m + 4, lrad + 1, 2): factor1 = float(n) / float(n ** 2 - m ** 2) factor2 = float(4 * (n - 1)) - factor3 = float((n - 2 + m) ** 2) / float(n - 2) + float( - (n - m) ** 2 - ) / float(n) + factor3 = float((n - 2 + m) ** 2) / float(n - 2) + float((n - m) ** 2) / float(n) factor4 = float((n - 2) ** 2 - m ** 2) / float(n - 2) zernike[:, n, m] = factor1 * ( @@ -348,9 +859,16 @@ def _get_zernike(sarr, lrad, mpol): + (factor2 * r ** 2 - factor3) * dzernike[:, n - 2, m] - factor4 * dzernike[:, n - 4, m] ) + d2zernike[:, n, m] = factor1 * ( + 2 * factor2 * zernike[:, n - 2, m] + + 4 * factor2 * r * dzernike[:, n - 2, m] + + (factor2 * r ** 2 - factor3) * d2zernike[:, n - 2, m] + - factor4 * d2zernike[:, n - 4, m] + ) + rm2 = rm1 rm1 = rm - rm = rm * r + rm = rm * r for n in range(2, lrad + 1, 2): zernike[:, n, 0] = zernike[:, n, 0] - (-1) ** (n / 2) @@ -364,11 +882,28 @@ def _get_zernike(sarr, lrad, mpol): (n + 1) / 2 ) + for m in range(mpol + 1): for n in range(m, lrad + 1, 2): - zernike[:, n, m] = zernike[:, n, m] / float(n + 1) - dzernike[:, n, m] = dzernike[:, n, m] / float(n + 1) + zernike[:, n, m] = zernike[:, n, m] / float(n + 1) + dzernike[:, n, m] = dzernike[:, n, m] / float(n + 1) + d2zernike[:, n, m] = d2zernike[:, n, m] / float(n + 1) + + dzernike = dzernike * 0.5 # to account for the factor of half in sbar = (1+s)/2 + d2zernike = d2zernike * 0.25 # to account for the factor of half in sbar = (1+s)/2 + + return zernike, dzernike, d2zernike + + +def invfft_B(B_mn_cos, B_mn_sin, mn, im, in_, Nfp, Ns, Nt, Nz): + B_mn = np.zeros(shape=(Ns,Nt,Nz), dtype=complex) - dzernike = dzernike * 0.5 # to account for the factor of half in sbar = (1+s)/2 + for imn in range(mn): + mm = im[imn] + nn = int(in_[imn]/Nfp) + B_mn[:, mm,-nn] = 0.5*(B_mn_cos[:,imn] - B_mn_sin[:,imn]*1j) + B_mn[:,-mm, nn] = 0.5*(B_mn_cos[:,imn] + B_mn_sin[:,imn]*1j) + B_mn[:,0,0] = 2*B_mn[:,0,0] + B = np.fft.irfft2(B_mn, s=(Nt,Nz))*Nt*Nz - return zernike, dzernike \ No newline at end of file + return B diff --git a/Utilities/pythontools/py_spec/output/spec.py b/Utilities/pythontools/py_spec/output/spec.py index 7eaf7fc2..c5689e63 100644 --- a/Utilities/pythontools/py_spec/output/spec.py +++ b/Utilities/pythontools/py_spec/output/spec.py @@ -34,7 +34,16 @@ class SPECout: grid, jacobian, metric, + get_B_FFT, + get_s_der_B_FFT, + get_t_der_B_FFT, + get_z_der_B_FFT, + get_B_and_der_FFT, get_B, + get_s_der_B, + get_t_der_B, + get_z_der_B, + get_B_and_der, get_modB, get_B_covariant ) diff --git a/src/brcast.f90 b/src/brcast.f90 index d4fc3695..9dfb9de6 100644 --- a/src/brcast.f90 +++ b/src/brcast.f90 @@ -134,7 +134,7 @@ subroutine brcast( lvol ) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! ! if( lvol.gt.Nvol .and. Lconstraint.eq.-1 .and. Wcurent ) then ! 27 Feb 17; - if( lvol.gt.Nvol .and. Wcurent ) then ! 27 Feb 17; + if( (lvol.gt.Nvol .or. Lconstraint .eq.-2) .and. Wcurent ) then ! 27 Feb 17; !write(ounit,'("brcast : " 10x " : myid="i3" ; broadcasting : curtor="es13.5" ; curpol="es13.5" ;")') myid, curtor, curpol RlBCAST( curtor, 1, llmodnp ) RlBCAST( curpol, 1, llmodnp ) diff --git a/src/curent.f90 b/src/curent.f90 index 87cd1255..1f268dee 100644 --- a/src/curent.f90 +++ b/src/curent.f90 @@ -63,7 +63,7 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) use fileunits, only : ounit - use inputlist, only : Wmacros, Wcurent, Lrad + use inputlist, only : Wmacros, Wcurent, Lrad, Lconstraint, Lbdybnzero use cputiming, only : Tcurent @@ -85,8 +85,8 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) INTEGER :: innout, ideriv, ii, ll, Lcurvature, ifail REAL :: lss - REAL :: Bsupt(1:Nt*Nz,-1:2), Bsupz(1:Nt*Nz,-1:2) - + REAL :: Bsupt(1:Nt*Nz,-1:2), Bsupz(1:Nt*Nz,-1:2), Bsups(1:Nt*Nz,-1:2), Bsups_2(1:Nt*Nz,-1:2) + BEGIN(curent) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! @@ -105,6 +105,11 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + if (Lconstraint .eq. -2) then + innout = 1.0 + lss = 1.0 + end if + do ideriv = -1, 2 ! labels derivative of magnetic field wrt enclosed fluxes; 20 Apr 13; if( iflag.eq. 1 .and. ideriv.ne.0 ) cycle ! derivatives of currents are not required; 20 Jun 14; @@ -116,6 +121,24 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) call invfft( mn, im(1:mn), in(1:mn), efmn(1:mn), ofmn(1:mn), cfmn(1:mn), sfmn(1:mn), & Nt, Nz, Bsupz(1:Ntz,ideriv), Bsupt(1:Ntz,ideriv) ) ! map to real space; + if (Lbdybnzero) then + Bsups = zero + else + call build_vector_potential(lvol, innout, ideriv, 0) + + do ii = 1, mn ! loop over Fourier harmonics; 17 May 21; + efmn(ii) = -in(ii)*efmn(ii) ! zeta derivative of Ate + ofmn(ii) = in(ii)*ofmn(ii) ! zeta derivative of Ato + cfmn(ii) = -im(ii)*cfmn(ii) ! theta derivative of Aze + sfmn(ii) = im(ii)*sfmn(ii) ! theta derivative of Azo + enddo ! end of do ii; + + call invfft( mn, im(1:mn), in(1:mn), ofmn(1:mn), cfmn(1:mn), sfmn(1:mn), efmn(1:mn), & + Nt, Nz, Bsups(1:Ntz,ideriv), Bsups_2(1:Ntz,ideriv)) + Bsups(1:Ntz,ideriv) = Bsups(1:Ntz,ideriv) + Bsups_2(1:ntz,ideriv) + + end if + enddo ! end of do ideriv; 31 Jan 13; !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! @@ -136,11 +159,22 @@ subroutine curent( lvol, mn, Nt, Nz, iflag, ldItGp ) if( iflag.eq. 2 .and. ideriv.lt.0 ) cycle ! derivatives of currents wrt geometry is not required; 20 Jun 14; if( iflag.eq.-1 .and. ideriv.gt.0 ) cycle ! derivatives of currents wrt enclosed toroidal and enclosed poloidal flux are not required; 20 Jun 14; - ijreal(1:Ntz) = ( - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,2,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) ) / sg(1:Ntz,0) - ijimag(1:Ntz) = ( - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,3,3,0) ) / sg(1:Ntz,0) + if (Lbdybnzero) then + ijreal(1:Ntz) = ( - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,2,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) ) / sg(1:Ntz,0) + ijimag(1:Ntz) = ( - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,3,3,0) ) / sg(1:Ntz,0) + else + ijreal(1:Ntz) = ( - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,2,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) + Bsups(1:Ntz,ideriv) * guvij(1:Ntz,2,1,0)) / sg(1:Ntz,0) + ijimag(1:Ntz) = ( - Bsupt(1:Ntz,ideriv) * guvij(1:Ntz,2,3,0) + Bsupz(1:Ntz,ideriv) * guvij(1:Ntz,3,3,0) + Bsups(1:Ntz,ideriv) * guvij(1:Ntz,1,3,0)) / sg(1:Ntz,0) + end if + if( ideriv.eq.-1 ) then ! add derivatives of metrics with respect to interface geometry; 15 Sep 16; - ijreal(1:Ntz) = ijreal(1:Ntz) + ( - Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,2,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,2,3,1) ) / sg(1:Ntz,0) - ijimag(1:Ntz) = ijimag(1:Ntz) + ( - Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,3,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,3,3,1) ) / sg(1:Ntz,0) + if (Lbdybnzero) then + ijreal(1:Ntz) = ijreal(1:Ntz) + ( - Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,2,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,2,3,1) ) / sg(1:Ntz,0) + ijimag(1:Ntz) = ijimag(1:Ntz) + ( - Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,3,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,3,3,1) ) / sg(1:Ntz,0) + else + ijreal(1:Ntz) = ijreal(1:Ntz) + ( - Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,2,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,2,3,1) + Bsups(1:Ntz, 0) * guvij(1:Ntz,2,1,1)) / sg(1:Ntz,0) + ijimag(1:Ntz) = ijimag(1:Ntz) + ( - Bsupt(1:Ntz, 0) * guvij(1:Ntz,2,3,1) + Bsupz(1:Ntz, 0) * guvij(1:Ntz,3,3,1) + Bsups(1:Ntz, 0) * guvij(1:Ntz,1,3,1)) / sg(1:Ntz,0) + end if endif ifail = 0 diff --git a/src/global.f90 b/src/global.f90 index 69003fc9..fd19cbcb 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -882,19 +882,27 @@ subroutine build_vector_potential(lvol, iocons, aderiv, tderiv) if( Lcoordinatesingularity ) then mi = im(ii) - do ll = mi, Lrad(lvol),2 ! loop over Zernike polynomials; Lrad is the radial resolution; 01 Jul 19; - ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half - ; ; cfmn(ii) = cfmn(ii) + Aze(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half - if( NOTstellsym ) then ; ofmn(ii) = ofmn(ii) + Ato(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half - ; ; sfmn(ii) = sfmn(ii) + Azo(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half - endif - enddo ! end of do ll; 20 Feb 13; + do ll = mi, Lrad(lvol),2 ! loop over Zernike polynomials; Lrad is the radial resolution; 01 Jul 19; + if( tderiv .eq. 1) then + ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half + ; ; cfmn(ii) = cfmn(ii) + Aze(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half + if( NOTstellsym ) then ; ofmn(ii) = ofmn(ii) + Ato(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half + ; ; sfmn(ii) = sfmn(ii) + Azo(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,1) * half + endif + else + ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) + ; ; cfmn(ii) = cfmn(ii) + Aze(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) + if( NOTstellsym ) then ; ofmn(ii) = ofmn(ii) + Ato(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) + ; ; sfmn(ii) = sfmn(ii) + Azo(lvol,aderiv,ii)%s(ll) * RTT(ll,mi,iocons,0) + endif + endif + enddo ! end of do ll; 20 Feb 13; else do ll = 0, Lrad(lvol) ! loop over Chebyshev polynomials; Lrad is the radial resolution; - ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * TT(ll,iocons,1) ! aderiv labels deriv. wrt mu, pflux; - ; ; cfmn(ii) = cfmn(ii) + Aze(lvol,aderiv,ii)%s(ll) * TT(ll,iocons,1) - if( NOTstellsym ) then ; ofmn(ii) = ofmn(ii) + Ato(lvol,aderiv,ii)%s(ll) * TT(ll,iocons,1) - ; ; sfmn(ii) = sfmn(ii) + Azo(lvol,aderiv,ii)%s(ll) * TT(ll,iocons,1) + ; ; efmn(ii) = efmn(ii) + Ate(lvol,aderiv,ii)%s(ll) * TT(ll,iocons,tderiv) ! aderiv labels deriv. wrt mu, pflux; + ; ; cfmn(ii) = cfmn(ii) + Aze(lvol,aderiv,ii)%s(ll) * TT(ll,iocons,tderiv) + if( NOTstellsym ) then ; ofmn(ii) = ofmn(ii) + Ato(lvol,aderiv,ii)%s(ll) * TT(ll,iocons,tderiv) + ; ; sfmn(ii) = sfmn(ii) + Azo(lvol,aderiv,ii)%s(ll) * TT(ll,iocons,tderiv) endif enddo ! end of do ll; 20 Feb 13; end if ! Lcoordinatesingularity; 01 Jul 19; @@ -1133,14 +1141,14 @@ subroutine check_inputs() cput = GETTIME - write(ounit,1010) cput-cpus, Igeometry, Istellsym, Lreflect + write(ounit,1010) cput-cpus, Igeometry, Istellsym, Lreflect, Lbdybnzero write(ounit,1011) Lfreebound, phiedge, curtor, curpol write(ounit,1012) gamma write(ounit,1013) Nfp, Nvol, Mvol, Mpol, Ntor write(ounit,1014) pscale, Ladiabatic, Lconstraint, mupftol, mupfits write(ounit,1015) Lrad(1:min(Mvol,32)) -1010 format("readin : ",f10.2," : Igeometry=",i3," ; Istellsym=",i3," ; Lreflect="i3" ;") +1010 format("readin : ",f10.2," : Igeometry=",i3," ; Istellsym=",i3," ; Lreflect="i3" ; Lbdybnzero="L2" ;") 1011 format("readin : ", 10x ," : Lfreebound=",i3," ; phiedge="es23.15" ; curtor="es23.15" ; curpol="es23.15" ;") 1012 format("readin : ", 10x ," : gamma="es23.15" ;") 1013 format("readin : ", 10x ," : Nfp=",i3," ; Nvol=",i3," ; Mvol=",i3," ; Mpol=",i3," ; Ntor=",i3," ;") @@ -1167,7 +1175,9 @@ subroutine check_inputs() FATAL( readin, mupftol.le.zero, mupftol is too small ) FATAL( readin, abs(one+gamma).lt.vsmall, 1+gamma appears in denominator in dforce ) !< \todo Please check this; SRH: 27 Feb 18; FATAL( readin, abs(one-gamma).lt.vsmall, 1-gamma appears in denominator in fu00aa ) !< \todo Please check this; SRH: 27 Feb 18; - FATAL( readin, Lconstraint.lt.-1 .or. Lconstraint.gt.3, illegal Lconstraint ) + FATAL( readin, Lconstraint.lt.-2 .or. Lconstraint.gt.3, illegal Lconstraint ) + FATAL( readin, Lconstraint.eq.-2 .and. Nvol.ne.1, Lconstraint=-2 only for single-volume calculation ) + FATAL( readin, .not. Lbdybnzero .and. (Lconstraint.ne.-2 .or. Lfreebound.eq.1) , Lbdybnzero=.false only for fixed-boundary calculation with Lconstraint=-2 ) FATAL( readin, Igeometry.eq.1 .and. rpol.lt.vsmall, poloidal extent of slab too small or negative ) FATAL( readin, Igeometry.eq.1 .and. rtor.lt.vsmall, toroidal extent of slab too small or negative ) @@ -1516,7 +1526,56 @@ subroutine broadcast_inputs if( Wreadin ) then ; cput = GETTIME ; write(ounit,'("readin : ",f10.2," : broadcasting screenlist from ext.sp ;")') cput-cpus endif -! BSCREENLIST ! broadcast screenlist; this is expanded by Makefile; do not remove; + LlBCAST( Wmanual ,1,0) + LlBCAST( Wrzaxis ,1,0) + LlBCAST( Wpackxi ,1,0) + LlBCAST( Wvolume ,1,0) + LlBCAST( Wcoords ,1,0) + LlBCAST( Wbasefn ,1,0) + LlBCAST( Wmemory ,1,0) + LlBCAST( Wmetrix ,1,0) + LlBCAST( Wma00aa ,1,0) + LlBCAST( Wmatrix ,1,0) + LlBCAST( Wspsmat ,1,0) + LlBCAST( Wspsint ,1,0) + LlBCAST( Wmp00ac ,1,0) + LlBCAST( Wma02aa ,1,0) + LlBCAST( Wpackab ,1,0) + LlBCAST( Wtr00ab ,1,0) + LlBCAST( Wcurent ,1,0) + LlBCAST( Wdf00ab ,1,0) + LlBCAST( Wlforce ,1,0) + LlBCAST( Wintghs ,1,0) + LlBCAST( Wmtrxhs ,1,0) + LlBCAST( Wlbpol ,1,0) + LlBCAST( Wbrcast ,1,0) + LlBCAST( Wdfp100 ,1,0) + LlBCAST( Wdfp200 ,1,0) + LlBCAST( Wdforce ,1,0) + LlBCAST( Wnewton ,1,0) + LlBCAST( Wcasing ,1,0) + LlBCAST( Wbnorml ,1,0) + LlBCAST( Wjo00aa ,1,0) + LlBCAST( Wpp00aa ,1,0) + LlBCAST( Wpp00ab ,1,0) + LlBCAST( Wbfield ,1,0) + LlBCAST( Wstzxyz ,1,0) + LlBCAST( Whesian ,1,0) + LlBCAST( Wra00aa ,1,0) + LlBCAST( Wnumrec ,1,0) + LlBCAST( Wdcuhre ,1,0) + LlBCAST( Wminpack,1,0) + LlBCAST( Wiqpack ,1,0) + LlBCAST( Wrksuite,1,0) + LlBCAST( Wi1mach ,1,0) + LlBCAST( Wd1mach ,1,0) + LlBCAST( Wilut ,1,0) + LlBCAST( Witers ,1,0) + LlBCAST( Wsphdf5 ,1,0) + LlBCAST( Wpreset ,1,0) + LlBCAST( Wglobal ,1,0) + LlBCAST( Wxspech ,1,0) + LlBCAST( Wbuild_vector_potential, 1, 0 ) LlBCAST( Wreadin, 1, 0 ) LlBCAST( Wwrtend, 1, 0 ) LlBCAST( Wmacros, 1, 0 ) @@ -1816,7 +1875,57 @@ subroutine wrtend endif write(iunit,'("&screenlist")') -! WSCREENLIST ! write screenlist; this is expanded by Makefile ; do not remove; + if( Wmanual ) write(iunit,'(" Wmanual = ",L1 )') Wmanual + if( Wrzaxis ) write(iunit,'(" Wrzaxis = ",L1 )') Wrzaxis + if( Wpackxi ) write(iunit,'(" Wpackxi = ",L1 )') Wpackxi + if( Wvolume ) write(iunit,'(" Wvolume = ",L1 )') Wvolume + if( Wcoords ) write(iunit,'(" Wcoords = ",L1 )') Wcoords + if( Wbasefn ) write(iunit,'(" Wbasefn = ",L1 )') Wbasefn + if( Wmemory ) write(iunit,'(" Wmemory = ",L1 )') Wmemory + if( Wmetrix ) write(iunit,'(" Wmetrix = ",L1 )') Wmetrix + if( Wma00aa ) write(iunit,'(" Wma00aa = ",L1 )') Wma00aa + if( Wmatrix ) write(iunit,'(" Wmatrix = ",L1 )') Wmatrix + if( Wspsmat ) write(iunit,'(" Wspsmat = ",L1 )') Wspsmat + if( Wspsint ) write(iunit,'(" Wspsint = ",L1 )') Wspsint + if( Wmp00ac ) write(iunit,'(" Wmp00ac = ",L1 )') Wmp00ac + if( Wma02aa ) write(iunit,'(" Wma02aa = ",L1 )') Wma02aa + if( Wpackab ) write(iunit,'(" Wpackab = ",L1 )') Wpackab + if( Wtr00ab ) write(iunit,'(" Wtr00ab = ",L1 )') Wtr00ab + if( Wcurent ) write(iunit,'(" Wcurent = ",L1 )') Wcurent + if( Wdf00ab ) write(iunit,'(" Wdf00ab = ",L1 )') Wdf00ab + if( Wlforce ) write(iunit,'(" Wlforce = ",L1 )') Wlforce + if( Wintghs ) write(iunit,'(" Wintghs = ",L1 )') Wintghs + if( Wmtrxhs ) write(iunit,'(" Wmtrxhs = ",L1 )') Wmtrxhs + if( Wlbpol ) write(iunit,'(" Wlbpol = ",L1 )') Wlbpol + if( Wbrcast ) write(iunit,'(" Wbrcast = ",L1 )') Wbrcast + if( Wdfp100 ) write(iunit,'(" Wdfp100 = ",L1 )') Wdfp100 + if( Wdfp200 ) write(iunit,'(" Wdfp200 = ",L1 )') Wdfp200 + if( Wdforce ) write(iunit,'(" Wdforce = ",L1 )') Wdforce + if( Wnewton ) write(iunit,'(" Wnewton = ",L1 )') Wnewton + if( Wcasing ) write(iunit,'(" Wcasing = ",L1 )') Wcasing + if( Wbnorml ) write(iunit,'(" Wbnorml = ",L1 )') Wbnorml + if( Wjo00aa ) write(iunit,'(" Wjo00aa = ",L1 )') Wjo00aa + if( Wpp00aa ) write(iunit,'(" Wpp00aa = ",L1 )') Wpp00aa + if( Wpp00ab ) write(iunit,'(" Wpp00ab = ",L1 )') Wpp00ab + if( Wbfield ) write(iunit,'(" Wbfield = ",L1 )') Wbfield + if( Wstzxyz ) write(iunit,'(" Wstzxyz = ",L1 )') Wstzxyz + if( Whesian ) write(iunit,'(" Whesian = ",L1 )') Whesian + if( Wra00aa ) write(iunit,'(" Wra00aa = ",L1 )') Wra00aa + if( Wnumrec ) write(iunit,'(" Wnumrec = ",L1 )') Wnumrec + if( Wdcuhre ) write(iunit,'(" Wdcuhre = ",L1 )') Wdcuhre + if( Wminpack ) write(iunit,'(" Wminpack= ",L1 )') Wminpack + if( Wiqpack ) write(iunit,'(" Wiqpack = ",L1 )') Wiqpack + if( Wrksuite ) write(iunit,'(" Wrksuite= ",L1 )') Wrksuite + if( Wi1mach ) write(iunit,'(" Wi1mach = ",L1 )') Wi1mach + if( Wd1mach ) write(iunit,'(" Wd1mach = ",L1 )') Wd1mach + if( Wilut ) write(iunit,'(" Wilut = ",L1 )') Wilut + if( Witers ) write(iunit,'(" Witers = ",L1 )') Witers + if( Wsphdf5 ) write(iunit,'(" Wsphdf5 = ",L1 )') Wsphdf5 + if( Wpreset ) write(iunit,'(" Wpreset = ",L1 )') Wpreset + if( Wglobal ) write(iunit,'(" Wglobal = ",L1 )') Wglobal + if( Wxspech ) write(iunit,'(" Wxspech = ",L1 )') Wxspech + if( Wbuild_vector_potential) write(iunit,'(" Wbuild_vector_potential = ",L1 )') Wbuild_vector_potential + if( Wreadin ) write(iunit,'(" Wreadin = ",L1 )') Wreadin if( Wwrtend ) write(iunit,'(" Wwrtend = ",L1 )') Wwrtend if( Wmacros ) write(iunit,'(" Wmacros = ",L1 )') Wmacros diff --git a/src/inputlist.f90 b/src/inputlist.f90 index b850a501..cc30f1b6 100644 --- a/src/inputlist.f90 +++ b/src/inputlist.f90 @@ -64,6 +64,7 @@ module inputlist !< INTEGER :: Lconstraint = -1 !< selects constraints; primarily used in ma02aa() and mp00ac(). !< INTEGER :: Lreflect = 0 !< =1 reflect the upper and lower bound in slab, =0 do not reflect + LOGICAL :: Lbdybnzero = .true. !< for fixed-boundary, =.true. assume Bsups=0 on boundary, =.false. obtain Bsups from Vns and Vnc input (currently only for Lconstraint=-2) REAL :: Rac( 0:MNtor ) = 0.0 !< stellarator symmetric coordinate axis; REAL :: Zas( 0:MNtor ) = 0.0 !< stellarator symmetric coordinate axis; @@ -640,6 +642,7 @@ module inputlist rpol ,& rtor ,& Lreflect ,& + Lbdybnzero ,& Rac ,& Zas ,& Ras ,& @@ -867,6 +870,7 @@ subroutine initialize_inputs mupfits = 8 Lreflect = 0 + Lbdybnzero = .true. ! numericlist diff --git a/src/ma02aa.f90 b/src/ma02aa.f90 index 4a175dc3..7b6d803b 100644 --- a/src/ma02aa.f90 +++ b/src/ma02aa.f90 @@ -368,13 +368,14 @@ subroutine ma02aa( lvol, NN ) !> !>