diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 2dd748a1..76738b54 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -69,6 +69,12 @@ jobs:
export OMP_NUM_THREADS=1
mpiexec -n 2 --allow-run-as-root ${SPEC_PATH}/xspec G3V08L3Fr.001.sp
python3 -m py_spec.ci.test compare.h5 G3V08L3Fr.001.sp.h5 --tol 1e-10
+ - name: current_constraint_free_boundary_grid
+ run: |
+ cd ${SPEC_PATH}/ci/G3V08L3FrGrid
+ export OMP_NUM_THREADS=1
+ mpiexec -n 2 --allow-run-as-root $SPEC_PATH/xspec G3V08L3FrGrid.001.sp
+ python3 -m py_spec.ci.test compare.h5 G3V08L3FrGrid.001.sp.h5 --tol 1e-10
diff --git a/.github/workflows/build_cmake.yml b/.github/workflows/build_cmake.yml
index 60496cd5..d710f539 100644
--- a/.github/workflows/build_cmake.yml
+++ b/.github/workflows/build_cmake.yml
@@ -75,4 +75,10 @@ jobs:
export OMP_NUM_THREADS=1
mpiexec -n 2 --allow-run-as-root $SPEC_PATH/install/bin/xspec G3V08L3Fr.001.sp
python3 -m py_spec.ci.test compare.h5 G3V08L3Fr.001.sp.h5 --tol 1e-10
+ - name: current_constraint_free_boundary_grid
+ run: |
+ cd ${SPEC_PATH}/ci/G3V08L3FrGrid
+ export OMP_NUM_THREADS=1
+ mpiexec -n 2 --allow-run-as-root $SPEC_PATH/install/bin/xspec G3V08L3FrGrid.001.sp
+ python3 -m py_spec.ci.test compare.h5 G3V08L3FrGrid.001.sp.h5 --tol 1e-10
diff --git a/ci/G3V08L3FrGrid/G3V08L3FrGrid.001.sp b/ci/G3V08L3FrGrid/G3V08L3FrGrid.001.sp
new file mode 100644
index 00000000..a434bdb0
--- /dev/null
+++ b/ci/G3V08L3FrGrid/G3V08L3FrGrid.001.sp
@@ -0,0 +1,234 @@
+&physicslist
+ Igeometry = 3
+ Istellsym = 1
+ Lfreebound = 1
+ phiedge = 1.000000000000000E+00
+ curtor = -1.187347930492713E-01
+ curpol = 4.881440045660601E+00
+ gamma = 0.000000000000000E+00
+ Nfp = 1
+ Nvol = 8
+ Mpol = 4
+ Ntor = 0
+ Lrad = 8 4 4 4 4 4 4 4 4
+ tflux = 1.562500000000000E-02 6.250000000000000E-02 1.406250000000000E-01 2.500000000000001E-01 3.906250000000002E-01 5.625000000000001E-01 7.656250000000000E-01 1.000000000000000E+00 2.091910274411301E+00
+ pflux = 0.000000000000000E+00 -2.229051459811318E-02 -5.797221981274356E-02 -1.048465195637780E-01 -1.598342781540251E-01 -2.189770364009819E-01 -2.774434484065295E-01 -3.295699704562904E-01 -5.063579643744914E-01
+ helicity = -1.185233159533041E-04 4.352416347064213E-06 2.247969531459685E-05 6.349418292781349E-05 1.365469458693930E-04 2.508072061865464E-04 4.154977940805720E-04 6.398895758928918E-04 2.108373382858917E-02
+ pscale = 1.256637076121100E-06
+ Ladiabatic = 0
+ pressure = 1.240234375000000E-01 1.201171875000000E-01 1.123046875000000E-01 1.005859375000000E-01 8.496093749999964E-02 6.542968750000001E-02 4.199218749999999E-02 1.464843750000000E-02 0.000000000000000E+00
+ adiabatic = 1.240234375000000E-01 1.201171875000000E-01 1.123046875000000E-01 1.005859375000000E-01 8.496093749999964E-02 6.542968750000001E-02 4.199218749999999E-02 1.464843750000000E-02 0.000000000000000E+00
+ mu = -2.505972582937235E-01 -2.423553768104948E-01 -2.258046346086362E-01 -2.007986034655939E-01 -1.671013994805084E-01 -1.244808992090318E-01 -7.325294758141898E-02 -1.695730962214170E-02 0.000000000000000E+00
+ Ivolume = -3.917350779014575E-03 -1.528292257746691E-02 -3.293210346280451E-02 -5.490505155132925E-02 -7.841584093870176E-02 -9.982370193182791E-02 -1.147156777262325E-01 -1.187033527604032E-01 -1.187033527604032E-01
+ Isurf = -1.493130820478333E-09 -1.032749836620444E-08 -3.437174496228381E-08 -8.418403698765600E-08 -1.749271148531898E-07 -3.359465605741845E-07 -6.368339541893372E-07 -6.214707460650634E-07 0.000000000000000E+00
+ Lconstraint = 3
+ pl = 0 0 0 0 0 0 0 0 0 0
+ ql = 0 0 0 0 0 0 0 0 0 0
+ pr = 0 0 0 0 0 0 0 0 0 0
+ qr = 0 0 0 0 0 0 0 0 0 0
+ iota = -2.603125000000000E-01 -4.853125000000000E-01 -4.712500000000000E-01 -4.478125000000000E-01 -4.150000000000000E-01 -3.728125000000000E-01 -3.212500000000000E-01 -2.603125000000000E-01 -1.900000000000000E-01 0.000000000000000E+00
+ lp = 0 0 0 0 0 0 0 0 0 0
+ lq = 0 0 0 0 0 0 0 0 0 0
+ rp = 0 0 0 0 0 0 0 0 0 0
+ rq = 0 0 0 0 0 0 0 0 0 0
+ oita = -2.603125000000000E-01 -4.853125000000000E-01 -4.712500000000000E-01 -4.478125000000000E-01 -4.150000000000000E-01 -3.728125000000000E-01 -3.212500000000000E-01 -2.603125000000000E-01 -1.900000000000000E-01 0.000000000000000E+00
+ mupftol = 1.000000000000000E-16
+ mupfits = 64
+ Rac = 4.095909789607840E+00
+ Zas = 0.000000000000000E+00
+ Ras = 0.000000000000000E+00
+ Zac = 0.000000000000000E+00
+Rbc(0,0) = 4.041135400356037E+00 Zbs(0,0) = 0.000000000000000E+00 Rbs(0,0) = 0.000000000000000E+00 Zbc(0,0) = 0.000000000000000E+00
+Rbc(0,1) = 1.038178094199237E+00 Zbs(0,1) = -1.561487214734595E+00 Rbs(0,1) = 0.000000000000000E+00 Zbc(0,1) = 0.000000000000000E+00
+Rbc(0,2) = 2.098569438650291E-02 Zbs(0,2) = 7.517833895284363E-02 Rbs(0,2) = 0.000000000000000E+00 Zbc(0,2) = 0.000000000000000E+00
+Rbc(0,3) = -1.556652677426258E-02 Zbs(0,3) = 1.367069428209015E-02 Rbs(0,3) = 0.000000000000000E+00 Zbc(0,3) = 0.000000000000000E+00
+Rbc(0,4) = 2.661700854262414E-03 Zbs(0,4) = -5.115432971414732E-03 Rbs(0,4) = 0.000000000000000E+00 Zbc(0,4) = 0.000000000000000E+00
+Rbc(0,5) = 4.785899100446688E-04 Zbs(0,5) = -2.367922532281356E-05 Rbs(0,5) = 0.000000000000000E+00 Zbc(0,5) = 0.000000000000000E+00
+Rbc(0,6) = -4.277170158982305E-04 Zbs(0,6) = 4.797132309430473E-04 Rbs(0,6) = 0.000000000000000E+00 Zbc(0,6) = 0.000000000000000E+00
+Rbc(0,7) = 5.543778295070250E-05 Zbs(0,7) = -8.072762565267903E-05 Rbs(0,7) = 0.000000000000000E+00 Zbc(0,7) = 0.000000000000000E+00
+Rbc(0,8) = 4.979036938393830E-05 Zbs(0,8) = -3.718554509294367E-05 Rbs(0,8) = 0.000000000000000E+00 Zbc(0,8) = 0.000000000000000E+00
+Rbc(0,9) = -1.796765283799183E-05 Zbs(0,9) = 1.536933967615198E-05 Rbs(0,9) = 0.000000000000000E+00 Zbc(0,9) = 0.000000000000000E+00
+Rbc(0,10) = -2.377661849947362E-06 Zbs(0,10) = -2.338589168549145E-06 Rbs(0,10) = 0.000000000000000E+00 Zbc(0,10) = 0.000000000000000E+00
+Rbc(0,11) = 3.551786680572178E-06 Zbs(0,11) = -4.338940150422311E-08 Rbs(0,11) = 0.000000000000000E+00 Zbc(0,11) = 0.000000000000000E+00
+Rbc(0,12) = -1.701788334984984E-06 Zbs(0,12) = 1.400088453818361E-06 Rbs(0,12) = 0.000000000000000E+00 Zbc(0,12) = 0.000000000000000E+00
+Rbc(0,13) = 4.333339215270045E-07 Zbs(0,13) = -1.119726848153810E-06 Rbs(0,13) = 0.000000000000000E+00 Zbc(0,13) = 0.000000000000000E+00
+Rbc(0,14) = 3.049252529910454E-07 Zbs(0,14) = -1.318268530889247E-09 Rbs(0,14) = 0.000000000000000E+00 Zbc(0,14) = 0.000000000000000E+00
+Rbc(0,15) = -3.604407231982447E-07 Zbs(0,15) = 3.913690673654217E-07 Rbs(0,15) = 0.000000000000000E+00 Zbc(0,15) = 0.000000000000000E+00
+Rbc(0,16) = 7.584655355530582E-08 Zbs(0,16) = -1.402111970447130E-07 Rbs(0,16) = 0.000000000000000E+00 Zbc(0,16) = 0.000000000000000E+00
+Rbc(0,17) = 8.002500606685238E-08 Zbs(0,17) = -6.474815340032061E-08 Rbs(0,17) = 0.000000000000000E+00 Zbc(0,17) = 0.000000000000000E+00
+Rbc(0,18) = -5.045745254025505E-08 Zbs(0,18) = 5.615707341427055E-08 Rbs(0,18) = 0.000000000000000E+00 Zbc(0,18) = 0.000000000000000E+00
+Rbc(0,19) = -4.742044039673486E-09 Zbs(0,19) = 3.180905965801948E-10 Rbs(0,19) = 0.000000000000000E+00 Zbc(0,19) = 0.000000000000000E+00
+Rbc(0,20) = 1.365790513800586E-08 Zbs(0,20) = -1.327752900952558E-08 Rbs(0,20) = 0.000000000000000E+00 Zbc(0,20) = 0.000000000000000E+00
+Rbc(0,21) = -2.491780260235909E-09 Zbs(0,21) = 3.226308626500275E-09 Rbs(0,21) = 0.000000000000000E+00 Zbc(0,21) = 0.000000000000000E+00
+Rbc(0,22) = -2.258175910797723E-09 Zbs(0,22) = 2.085280130938560E-09 Rbs(0,22) = 0.000000000000000E+00 Zbc(0,22) = 0.000000000000000E+00
+Rbc(0,23) = 9.035538509789723E-10 Zbs(0,23) = -1.103246648737506E-09 Rbs(0,23) = 0.000000000000000E+00 Zbc(0,23) = 0.000000000000000E+00
+Rbc(0,24) = 2.521910786093891E-10 Zbs(0,24) = -1.931113582046620E-10 Rbs(0,24) = 0.000000000000000E+00 Zbc(0,24) = 0.000000000000000E+00
+Rbc(0,25) = -1.544586175302791E-10 Zbs(0,25) = 2.366410258976099E-10 Rbs(0,25) = 0.000000000000000E+00 Zbc(0,25) = 0.000000000000000E+00
+Rbc(0,26) = -3.364106577785311E-11 Zbs(0,26) = -3.467444317281823E-12 Rbs(0,26) = 0.000000000000000E+00 Zbc(0,26) = 0.000000000000000E+00
+Rbc(0,27) = 5.002585191818008E-12 Zbs(0,27) = -3.355899105022525E-11 Rbs(0,27) = 0.000000000000000E+00 Zbc(0,27) = 0.000000000000000E+00
+Rbc(0,28) = 1.651166828997631E-11 Zbs(0,28) = 5.085494430080715E-12 Rbs(0,28) = 0.000000000000000E+00 Zbc(0,28) = 0.000000000000000E+00
+Rbc(0,29) = 1.360704678907421E-11 Zbs(0,29) = 2.322075285404156E-12 Rbs(0,29) = 0.000000000000000E+00 Zbc(0,29) = 0.000000000000000E+00
+Rbc(0,30) = -2.041475056207957E-11 Zbs(0,30) = -8.039344124784536E-12 Rbs(0,30) = 0.000000000000000E+00 Zbc(0,30) = 0.000000000000000E+00
+Rbc(0,31) = -7.244505903823517E-12 Zbs(0,31) = 4.432435918368827E-13 Rbs(0,31) = 0.000000000000000E+00 Zbc(0,31) = 0.000000000000000E+00
+Rbc(0,32) = 3.134916792674994E-11 Zbs(0,32) = 2.843559466998492E-11 Rbs(0,32) = 0.000000000000000E+00 Zbc(0,32) = 0.000000000000000E+00
+Rwc(0,0) = 3.999000000000000E+00 Zws(0,0) = 0.000000000000000E+00 Rws(0,0) = 0.000000000000000E+00 Zwc(0,0) = 0.000000000000000E+00
+Rwc(0,1) = 1.800000000000000E+00 Zws(0,1) = -1.800000000000000E+00 Rws(0,1) = 0.000000000000000E+00 Zwc(0,1) = 0.000000000000000E+00
+Rwc(0,2) = 0.000000000000000E+00 Zws(0,2) = 0.000000000000000E+00 Rws(0,2) = 0.000000000000000E+00 Zwc(0,2) = 0.000000000000000E+00
+Rwc(0,3) = 0.000000000000000E+00 Zws(0,3) = 0.000000000000000E+00 Rws(0,3) = 0.000000000000000E+00 Zwc(0,3) = 0.000000000000000E+00
+Rwc(0,4) = 0.000000000000000E+00 Zws(0,4) = 0.000000000000000E+00 Rws(0,4) = 0.000000000000000E+00 Zwc(0,4) = 0.000000000000000E+00
+Rwc(0,5) = 0.000000000000000E+00 Zws(0,5) = 0.000000000000000E+00 Rws(0,5) = 0.000000000000000E+00 Zwc(0,5) = 0.000000000000000E+00
+Rwc(0,6) = 0.000000000000000E+00 Zws(0,6) = 0.000000000000000E+00 Rws(0,6) = 0.000000000000000E+00 Zwc(0,6) = 0.000000000000000E+00
+Rwc(0,7) = 0.000000000000000E+00 Zws(0,7) = 0.000000000000000E+00 Rws(0,7) = 0.000000000000000E+00 Zwc(0,7) = 0.000000000000000E+00
+Rwc(0,8) = 0.000000000000000E+00 Zws(0,8) = 0.000000000000000E+00 Rws(0,8) = 0.000000000000000E+00 Zwc(0,8) = 0.000000000000000E+00
+Rwc(0,9) = 0.000000000000000E+00 Zws(0,9) = 0.000000000000000E+00 Rws(0,9) = 0.000000000000000E+00 Zwc(0,9) = 0.000000000000000E+00
+Rwc(0,10) = 0.000000000000000E+00 Zws(0,10) = 0.000000000000000E+00 Rws(0,10) = 0.000000000000000E+00 Zwc(0,10) = 0.000000000000000E+00
+Rwc(0,11) = 0.000000000000000E+00 Zws(0,11) = 0.000000000000000E+00 Rws(0,11) = 0.000000000000000E+00 Zwc(0,11) = 0.000000000000000E+00
+Rwc(0,12) = 0.000000000000000E+00 Zws(0,12) = 0.000000000000000E+00 Rws(0,12) = 0.000000000000000E+00 Zwc(0,12) = 0.000000000000000E+00
+Rwc(0,13) = 0.000000000000000E+00 Zws(0,13) = 0.000000000000000E+00 Rws(0,13) = 0.000000000000000E+00 Zwc(0,13) = 0.000000000000000E+00
+Rwc(0,14) = 0.000000000000000E+00 Zws(0,14) = 0.000000000000000E+00 Rws(0,14) = 0.000000000000000E+00 Zwc(0,14) = 0.000000000000000E+00
+Rwc(0,15) = 0.000000000000000E+00 Zws(0,15) = 0.000000000000000E+00 Rws(0,15) = 0.000000000000000E+00 Zwc(0,15) = 0.000000000000000E+00
+Rwc(0,16) = 0.000000000000000E+00 Zws(0,16) = 0.000000000000000E+00 Rws(0,16) = 0.000000000000000E+00 Zwc(0,16) = 0.000000000000000E+00
+Rwc(0,17) = 0.000000000000000E+00 Zws(0,17) = 0.000000000000000E+00 Rws(0,17) = 0.000000000000000E+00 Zwc(0,17) = 0.000000000000000E+00
+Rwc(0,18) = 0.000000000000000E+00 Zws(0,18) = 0.000000000000000E+00 Rws(0,18) = 0.000000000000000E+00 Zwc(0,18) = 0.000000000000000E+00
+Rwc(0,19) = 0.000000000000000E+00 Zws(0,19) = 0.000000000000000E+00 Rws(0,19) = 0.000000000000000E+00 Zwc(0,19) = 0.000000000000000E+00
+Rwc(0,20) = 0.000000000000000E+00 Zws(0,20) = 0.000000000000000E+00 Rws(0,20) = 0.000000000000000E+00 Zwc(0,20) = 0.000000000000000E+00
+Rwc(0,21) = 0.000000000000000E+00 Zws(0,21) = 0.000000000000000E+00 Rws(0,21) = 0.000000000000000E+00 Zwc(0,21) = 0.000000000000000E+00
+Rwc(0,22) = 0.000000000000000E+00 Zws(0,22) = 0.000000000000000E+00 Rws(0,22) = 0.000000000000000E+00 Zwc(0,22) = 0.000000000000000E+00
+Rwc(0,23) = 0.000000000000000E+00 Zws(0,23) = 0.000000000000000E+00 Rws(0,23) = 0.000000000000000E+00 Zwc(0,23) = 0.000000000000000E+00
+Rwc(0,24) = 0.000000000000000E+00 Zws(0,24) = 0.000000000000000E+00 Rws(0,24) = 0.000000000000000E+00 Zwc(0,24) = 0.000000000000000E+00
+Rwc(0,25) = 0.000000000000000E+00 Zws(0,25) = 0.000000000000000E+00 Rws(0,25) = 0.000000000000000E+00 Zwc(0,25) = 0.000000000000000E+00
+Rwc(0,26) = 0.000000000000000E+00 Zws(0,26) = 0.000000000000000E+00 Rws(0,26) = 0.000000000000000E+00 Zwc(0,26) = 0.000000000000000E+00
+Rwc(0,27) = 0.000000000000000E+00 Zws(0,27) = 0.000000000000000E+00 Rws(0,27) = 0.000000000000000E+00 Zwc(0,27) = 0.000000000000000E+00
+Rwc(0,28) = 0.000000000000000E+00 Zws(0,28) = 0.000000000000000E+00 Rws(0,28) = 0.000000000000000E+00 Zwc(0,28) = 0.000000000000000E+00
+Rwc(0,29) = 0.000000000000000E+00 Zws(0,29) = 0.000000000000000E+00 Rws(0,29) = 0.000000000000000E+00 Zwc(0,29) = 0.000000000000000E+00
+Rwc(0,30) = 0.000000000000000E+00 Zws(0,30) = 0.000000000000000E+00 Rws(0,30) = 0.000000000000000E+00 Zwc(0,30) = 0.000000000000000E+00
+Rwc(0,31) = 0.000000000000000E+00 Zws(0,31) = 0.000000000000000E+00 Rws(0,31) = 0.000000000000000E+00 Zwc(0,31) = 0.000000000000000E+00
+Rwc(0,32) = 0.000000000000000E+00 Zws(0,32) = 0.000000000000000E+00 Rws(0,32) = 0.000000000000000E+00 Zwc(0,32) = 0.000000000000000E+00
+Vns(0,0) = 0.000000000000000E+00 Bns(0,0) = 0.000000000000000E+00 Vnc(0,0) = 0.000000000000000E+00 Bnc(0,0) = 0.000000000000000E+00
+Vns(0,1) = -3.794877659624626E-02 Bns(0,1) = 3.216851566291370E-02 Vnc(0,1) = 0.000000000000000E+00 Bnc(0,1) = 0.000000000000000E+00
+Vns(0,2) = -3.922424672833572E-02 Bns(0,2) = -5.715031890662380E-03 Vnc(0,2) = 0.000000000000000E+00 Bnc(0,2) = 0.000000000000000E+00
+Vns(0,3) = 9.308639007382408E-03 Bns(0,3) = -5.739913236142605E-04 Vnc(0,3) = 0.000000000000000E+00 Bnc(0,3) = 0.000000000000000E+00
+Vns(0,4) = -4.028451409448938E-03 Bns(0,4) = 8.204157785357483E-04 Vnc(0,4) = 0.000000000000000E+00 Bnc(0,4) = 0.000000000000000E+00
+Vns(0,5) = 3.313589443914825E-03 Bns(0,5) = -2.953431815673370E-05 Vnc(0,5) = 0.000000000000000E+00 Bnc(0,5) = 0.000000000000000E+00
+Vns(0,6) = -5.532607453983990E-04 Bns(0,6) = -1.668386253319423E-04 Vnc(0,6) = 0.000000000000000E+00 Bnc(0,6) = 0.000000000000000E+00
+Vns(0,7) = 3.428218994380060E-05 Bns(0,7) = 3.234579790990208E-05 Vnc(0,7) = 0.000000000000000E+00 Bnc(0,7) = 0.000000000000000E+00
+Vns(0,8) = 8.337672863031870E-05 Bns(0,8) = 3.995575596057115E-05 Vnc(0,8) = 0.000000000000000E+00 Bnc(0,8) = 0.000000000000000E+00
+Vns(0,9) = -1.474119717292630E-05 Bns(0,9) = -1.639792815785964E-05 Vnc(0,9) = 0.000000000000000E+00 Bnc(0,9) = 0.000000000000000E+00
+Vns(0,10) = 6.218008236436993E-05 Bns(0,10) = -9.455366957759755E-06 Vnc(0,10) = 0.000000000000000E+00 Bnc(0,10) = 0.000000000000000E+00
+Vns(0,11) = -5.844833566026721E-05 Bns(0,11) = 7.056193435980942E-06 Vnc(0,11) = 0.000000000000000E+00 Bnc(0,11) = 0.000000000000000E+00
+Vns(0,12) = 2.524161994086079E-05 Bns(0,12) = 1.943931173292392E-06 Vnc(0,12) = 0.000000000000000E+00 Bnc(0,12) = 0.000000000000000E+00
+Vns(0,13) = -8.506003434154704E-06 Bns(0,13) = -2.865574135919580E-06 Vnc(0,13) = 0.000000000000000E+00 Bnc(0,13) = 0.000000000000000E+00
+Vns(0,14) = 2.346227552948437E-06 Bns(0,14) = -1.801492747585585E-07 Vnc(0,14) = 0.000000000000000E+00 Bnc(0,14) = 0.000000000000000E+00
+Vns(0,15) = 7.951661703768458E-07 Bns(0,15) = 1.113786392891701E-06 Vnc(0,15) = 0.000000000000000E+00 Bnc(0,15) = 0.000000000000000E+00
+Vns(0,16) = -1.745853525635911E-06 Bns(0,16) = -1.383304516545098E-07 Vnc(0,16) = 0.000000000000000E+00 Bnc(0,16) = 0.000000000000000E+00
+Vns(0,17) = 1.078835737470322E-06 Bns(0,17) = -4.125555850075048E-07 Vnc(0,17) = 0.000000000000000E+00 Bnc(0,17) = 0.000000000000000E+00
+Vns(0,18) = -3.939953204899884E-07 Bns(0,18) = 1.321286689483563E-07 Vnc(0,18) = 0.000000000000000E+00 Bnc(0,18) = 0.000000000000000E+00
+Vns(0,19) = 1.217144462573904E-07 Bns(0,19) = 1.430195017347739E-07 Vnc(0,19) = 0.000000000000000E+00 Bnc(0,19) = 0.000000000000000E+00
+Vns(0,20) = -1.902155301364819E-08 Bns(0,20) = -8.078837184863531E-08 Vnc(0,20) = 0.000000000000000E+00 Bnc(0,20) = 0.000000000000000E+00
+Vns(0,21) = -2.648040478570356E-08 Bns(0,21) = -4.453358659624485E-08 Vnc(0,21) = 0.000000000000000E+00 Bnc(0,21) = 0.000000000000000E+00
+Vns(0,22) = 2.779878514149091E-08 Bns(0,22) = 4.207837878333624E-08 Vnc(0,22) = 0.000000000000000E+00 Bnc(0,22) = 0.000000000000000E+00
+Vns(0,23) = -1.375607517619595E-08 Bns(0,23) = 1.111051743069157E-08 Vnc(0,23) = 0.000000000000000E+00 Bnc(0,23) = 0.000000000000000E+00
+Vns(0,24) = 5.170475964997598E-09 Bns(0,24) = -1.995212214569733E-08 Vnc(0,24) = 0.000000000000000E+00 Bnc(0,24) = 0.000000000000000E+00
+Vns(0,25) = -1.832238887978118E-09 Bns(0,25) = -1.115699430794461E-09 Vnc(0,25) = 0.000000000000000E+00 Bnc(0,25) = 0.000000000000000E+00
+Vns(0,26) = 1.567041550057399E-10 Bns(0,26) = 8.799897560316338E-09 Vnc(0,26) = 0.000000000000000E+00 Bnc(0,26) = 0.000000000000000E+00
+Vns(0,27) = 4.746636964325032E-10 Bns(0,27) = -1.119834921835541E-09 Vnc(0,27) = 0.000000000000000E+00 Bnc(0,27) = 0.000000000000000E+00
+Vns(0,28) = -3.901756370229985E-10 Bns(0,28) = -3.623166913106286E-09 Vnc(0,28) = 0.000000000000000E+00 Bnc(0,28) = 0.000000000000000E+00
+Vns(0,29) = 1.874047401481348E-10 Bns(0,29) = 1.147748890119685E-09 Vnc(0,29) = 0.000000000000000E+00 Bnc(0,29) = 0.000000000000000E+00
+Vns(0,30) = -8.053571015120110E-11 Bns(0,30) = 1.379073565906986E-09 Vnc(0,30) = 0.000000000000000E+00 Bnc(0,30) = 0.000000000000000E+00
+Vns(0,31) = 2.769273393977595E-11 Bns(0,31) = -7.551711113767656E-10 Vnc(0,31) = 0.000000000000000E+00 Bnc(0,31) = 0.000000000000000E+00
+Vns(0,32) = 9.623711250383996E-13 Bns(0,32) = -4.709530705132861E-10 Vnc(0,32) = 0.000000000000000E+00 Bnc(0,32) = 0.000000000000000E+00
+/
+&numericlist
+ Linitialize = 0
+ Lzerovac = 0
+ Ndiscrete = 12
+ Nquad = -12
+ iMpol = -12
+ iNtor = -12
+ Lsparse = 0
+ Lsvdiota = 1
+ imethod = 3
+ iorder = 2
+ iprecon = 1
+ iotatol = -1.000000000000000E+00
+ Lextrap = 0
+ Mregular = -1
+ LautoinitBn = 0
+ Lvcgrid = 1
+/
+&locallist
+ LBeltrami = 4
+ Linitgues = 1
+ epsGMRES = 1e-15
+/
+&globallist
+ Lfindzero = 2
+ escale = 0.000000000000000E+00
+ opsilon = 1.000000000000000E+00
+ pcondense = 4.000000000000000E+00
+ epsilon = 1.000000000000000E-01
+ wpoloidal = 1.000000000000000E+00
+ upsilon = 1.000000000000000E-01
+ forcetol = 1.000000000000000E-13
+ c05xmax = 1.000000000000000E-06
+ c05xtol = 1.000000000000000E-13
+ c05factor = 1.000000000000000E-04
+ LreadGF = T
+ mfreeits = 100
+ gBntol = 1.000000000000000E-06
+ gBnbld = 0.000000000000000E+00
+ vcasingeps = 0.000000000000000E+00
+ vcasingtol = 1.000000000000000E-12
+ vcasingits = 1000
+ vcasingper = 1
+ vcNt = 360
+ vcNz = 412
+/
+&diagnosticslist
+ odetol = 1.000000000000000E-07
+ nPpts = 200
+ nPtrj = 3 3 3 3 3 3 3 3 3
+ LHevalues = F
+ LHevectors = F
+ LHmatrix = F
+ Lperturbed = 0
+ dpp = -1
+ dqq = -1
+ Lcheck = 0
+ Ltiming = F
+/
+&screenlist
+/
+ 0 0 4.093334585563770E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 4.092446004207707E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 4.090650609897955E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 4.087528763035849E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 4.082401351391629E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 4.074235151822619E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 4.061433938496663E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 4.041135400356037E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00
+ 1 0 1.365592450388462E-01 -1.900398693544442E-01 0.000000000000000E+00 0.000000000000000E+00 2.730078337399156E-01 -3.798783139565216E-01 0.000000000000000E+00 0.000000000000000E+00 4.085003472248540E-01 -5.703176833819588E-01 0.000000000000000E+00 0.000000000000000E+00 5.423583664687349E-01 -7.618722797471389E-01 0.000000000000000E+00 0.000000000000000E+00 6.736364994916200E-01 -9.553690039318099E-01 0.000000000000000E+00 0.000000000000000E+00 8.008988051688418E-01 -1.152289760829545E+00 0.000000000000000E+00 0.000000000000000E+00 9.221469177407754E-01 -1.354973794988916E+00 0.000000000000000E+00 0.000000000000000E+00 1.038178094199237E+00 -1.561487214734595E+00 0.000000000000000E+00 0.000000000000000E+00
+ 2 0 -1.554813308941516E-04 1.204541228834543E-03 0.000000000000000E+00 0.000000000000000E+00 -6.510695836783213E-04 4.891490533915513E-03 0.000000000000000E+00 0.000000000000000E+00 -1.199614828150767E-03 1.084529981955521E-02 0.000000000000000E+00 0.000000000000000E+00 -1.437734841859040E-03 1.891420691173961E-02 0.000000000000000E+00 0.000000000000000E+00 -7.653046650349492E-04 2.892168625277512E-02 0.000000000000000E+00 0.000000000000000E+00 1.765219565789863E-03 4.085080723993717E-02 0.000000000000000E+00 0.000000000000000E+00 7.790321294680378E-03 5.543276433193616E-02 0.000000000000000E+00 0.000000000000000E+00 2.098569438650291E-02 7.517833895284363E-02 0.000000000000000E+00 0.000000000000000E+00
+ 3 0 -1.583494662191774E-04 1.696168858215460E-04 0.000000000000000E+00 0.000000000000000E+00 -7.574125089835899E-04 8.099487780575853E-04 0.000000000000000E+00 0.000000000000000E+00 -1.438566893981201E-03 1.336535599498400E-03 0.000000000000000E+00 0.000000000000000E+00 -2.240119053221266E-03 1.685155509498641E-03 0.000000000000000E+00 0.000000000000000E+00 -3.203573890044584E-03 1.811218149433840E-03 0.000000000000000E+00 0.000000000000000E+00 -4.402042812300572E-03 1.750289471346588E-03 0.000000000000000E+00 0.000000000000000E+00 -6.422105575519354E-03 2.448565971247674E-03 0.000000000000000E+00 0.000000000000000E+00 -1.556652677426258E-02 1.367069428209015E-02 0.000000000000000E+00 0.000000000000000E+00
+ 4 0 3.097134291174607E-05 -3.895173513682294E-05 0.000000000000000E+00 0.000000000000000E+00 1.425657025694941E-04 -1.824249772845022E-04 0.000000000000000E+00 0.000000000000000E+00 3.437261639917861E-04 -4.218158729888696E-04 0.000000000000000E+00 0.000000000000000E+00 6.545779028408568E-04 -7.468611596646806E-04 0.000000000000000E+00 0.000000000000000E+00 1.103216584871905E-03 -1.151434081030379E-03 0.000000000000000E+00 0.000000000000000E+00 1.715077801299083E-03 -1.651396285646003E-03 0.000000000000000E+00 0.000000000000000E+00 2.423153876698911E-03 -2.381434902410343E-03 0.000000000000000E+00 0.000000000000000E+00 2.661700854262414E-03 -5.115432971414732E-03 0.000000000000000E+00 0.000000000000000E+00
+ 5 0 8.193539404358615E-06 -9.929829168638184E-06 0.000000000000000E+00 0.000000000000000E+00 2.395461189927865E-05 -2.492453141690717E-05 0.000000000000000E+00 0.000000000000000E+00 2.682863868497100E-05 -2.191354231703629E-05 0.000000000000000E+00 0.000000000000000E+00 1.038493305970101E-05 7.174737833127007E-06 0.000000000000000E+00 0.000000000000000E+00 -3.141333305293069E-05 6.640619977581583E-05 0.000000000000000E+00 0.000000000000000E+00 -9.921997934822642E-05 1.569225710439943E-04 0.000000000000000E+00 0.000000000000000E+00 -1.589690080599841E-04 2.812119634769483E-04 0.000000000000000E+00 0.000000000000000E+00 4.785899100446688E-04 -2.367922532281356E-05 0.000000000000000E+00 0.000000000000000E+00
+ 6 0 -2.090036284807985E-06 2.592257442260072E-06 0.000000000000000E+00 0.000000000000000E+00 -9.773967554198338E-06 1.171691906367510E-05 0.000000000000000E+00 0.000000000000000E+00 -2.170681073173839E-05 2.511305448924686E-05 0.000000000000000E+00 0.000000000000000E+00 -3.599515498910216E-05 3.908448776317121E-05 0.000000000000000E+00 0.000000000000000E+00 -5.176809526529286E-05 4.966714118129494E-05 0.000000000000000E+00 0.000000000000000E+00 -7.350631719203419E-05 5.626020406190869E-05 0.000000000000000E+00 0.000000000000000E+00 -1.326531108634506E-04 8.322231194216734E-05 0.000000000000000E+00 0.000000000000000E+00 -4.277170158982305E-04 4.797132309430473E-04 0.000000000000000E+00 0.000000000000000E+00
+ 7 0 -3.952973731518395E-07 4.295656472737399E-07 0.000000000000000E+00 0.000000000000000E+00 -6.834664774369463E-07 5.867441389282066E-07 0.000000000000000E+00 0.000000000000000E+00 6.341349817908607E-07 -1.389568593410887E-06 0.000000000000000E+00 0.000000000000000E+00 4.686885235532178E-06 -6.118587259671001E-06 0.000000000000000E+00 0.000000000000000E+00 1.244833366619190E-05 -1.354656476025734E-05 0.000000000000000E+00 0.000000000000000E+00 2.630828859535331E-05 -2.346908263614810E-05 0.000000000000000E+00 0.000000000000000E+00 5.219450183729755E-05 -3.853704988041856E-05 0.000000000000000E+00 0.000000000000000E+00 5.543778295070250E-05 -8.072762565267903E-05 0.000000000000000E+00 0.000000000000000E+00
+ 8 0 1.771841544593403E-07 -2.047089339184605E-07 0.000000000000000E+00 0.000000000000000E+00 8.193257094581209E-07 -9.430795857116912E-07 0.000000000000000E+00 0.000000000000000E+00 1.703172984254098E-06 -1.809440443474689E-06 0.000000000000000E+00 0.000000000000000E+00 2.273201259817857E-06 -2.211464920681497E-06 0.000000000000000E+00 0.000000000000000E+00 1.782030585508849E-06 -1.319863665740324E-06 0.000000000000000E+00 0.000000000000000E+00 -5.689455875398658E-07 1.810986132489908E-06 0.000000000000000E+00 0.000000000000000E+00 -2.953719358817113E-06 6.725227979328033E-06 0.000000000000000E+00 0.000000000000000E+00 4.979036938393830E-05 -3.718554509294367E-05 0.000000000000000E+00 0.000000000000000E+00
+ 9 0 3.363315956425363E-08 -3.784831883171022E-08 0.000000000000000E+00 0.000000000000000E+00 -1.518093382108205E-09 1.282270664620430E-09 0.000000000000000E+00 0.000000000000000E+00 -2.602387005345098E-07 3.208566916096165E-07 0.000000000000000E+00 0.000000000000000E+00 -7.329209207871484E-07 8.377705998811014E-07 0.000000000000000E+00 0.000000000000000E+00 -1.277452217173141E-06 1.302484355210795E-06 0.000000000000000E+00 0.000000000000000E+00 -1.636088767249279E-06 1.181372244363894E-06 0.000000000000000E+00 0.000000000000000E+00 -2.800764854716669E-06 5.222890829103212E-07 0.000000000000000E+00 0.000000000000000E+00 -1.796765283799183E-05 1.536933967615198E-05 0.000000000000000E+00 0.000000000000000E+00
+ 10 0 -1.938707786273855E-08 2.225673443022073E-08 0.000000000000000E+00 0.000000000000000E+00 -8.421705052763909E-08 9.733618387100236E-08 0.000000000000000E+00 0.000000000000000E+00 -1.338825621763157E-07 1.443018205535782E-07 0.000000000000000E+00 0.000000000000000E+00 -1.147143741873534E-07 9.246724129509841E-08 0.000000000000000E+00 0.000000000000000E+00 6.556195463735405E-08 -1.358031213960642E-07 0.000000000000000E+00 0.000000000000000E+00 5.112313247368755E-07 -6.518968440854984E-07 0.000000000000000E+00 0.000000000000000E+00 1.489009221256456E-06 -1.961467704030756E-06 0.000000000000000E+00 0.000000000000000E+00 -2.377661849947362E-06 -2.338589168549145E-06 0.000000000000000E+00 0.000000000000000E+00
+ 11 0 -2.480375200784727E-09 2.597929417570308E-09 0.000000000000000E+00 0.000000000000000E+00 4.032360369560990E-09 -6.719227558479267E-09 0.000000000000000E+00 0.000000000000000E+00 4.716764951286963E-08 -5.224160262077968E-08 0.000000000000000E+00 0.000000000000000E+00 1.019970496111408E-07 -1.089419447724858E-07 0.000000000000000E+00 0.000000000000000E+00 1.296288318476605E-07 -1.161296865827370E-07 0.000000000000000E+00 0.000000000000000E+00 4.571641195471396E-08 6.583554980744962E-08 0.000000000000000E+00 0.000000000000000E+00 -1.793853325907110E-07 8.041261331489185E-07 0.000000000000000E+00 0.000000000000000E+00 3.551786680572178E-06 -4.338940150422311E-08 0.000000000000000E+00 0.000000000000000E+00
+ 12 0 2.175315531662830E-09 -2.413332973426469E-09 0.000000000000000E+00 0.000000000000000E+00 9.758575146184670E-09 -1.037822873787438E-08 0.000000000000000E+00 0.000000000000000E+00 9.299435405526318E-09 -1.097981104528593E-08 0.000000000000000E+00 0.000000000000000E+00 -6.924687923417384E-09 8.629215556796198E-09 0.000000000000000E+00 0.000000000000000E+00 -4.330016788581699E-08 4.828436869747123E-08 0.000000000000000E+00 0.000000000000000E+00 -1.123742177173461E-07 1.029950061867391E-07 0.000000000000000E+00 0.000000000000000E+00 -3.408958011815419E-07 2.143441767289206E-07 0.000000000000000E+00 0.000000000000000E+00 -1.701788334984984E-06 1.400088453818361E-06 0.000000000000000E+00 0.000000000000000E+00
+ 13 0 2.215841696661710E-10 -2.439242166251912E-10 0.000000000000000E+00 0.000000000000000E+00 -1.442596273513140E-09 1.671023695346469E-09 0.000000000000000E+00 0.000000000000000E+00 -5.948599125839664E-09 6.880562502747422E-09 0.000000000000000E+00 0.000000000000000E+00 -1.056858356539000E-08 1.124629852633377E-08 0.000000000000000E+00 0.000000000000000E+00 -4.036909260369123E-09 1.765368100770664E-09 0.000000000000000E+00 0.000000000000000E+00 3.973908455939783E-08 -5.154316066994999E-08 0.000000000000000E+00 0.000000000000000E+00 2.225045538991286E-07 -2.735080457099680E-07 0.000000000000000E+00 0.000000000000000E+00 4.333339215270045E-07 -1.119726848153810E-06 0.000000000000000E+00 0.000000000000000E+00
+ 14 0 -2.732736934847214E-10 3.044989451281387E-10 0.000000000000000E+00 0.000000000000000E+00 -1.022949542877850E-09 1.109945162631398E-09 0.000000000000000E+00 0.000000000000000E+00 -1.003646533827945E-09 9.350573346625307E-10 0.000000000000000E+00 0.000000000000000E+00 2.372028123072966E-09 -2.266570558757585E-09 0.000000000000000E+00 0.000000000000000E+00 6.846094083120364E-09 -6.373773569979845E-09 0.000000000000000E+00 0.000000000000000E+00 6.454417577513094E-09 -1.745418522629735E-09 0.000000000000000E+00 0.000000000000000E+00 -9.625870922988825E-09 4.176386378900927E-08 0.000000000000000E+00 0.000000000000000E+00 3.049252529910454E-07 -1.318268530889247E-09 0.000000000000000E+00 0.000000000000000E+00
+ 15 0 -1.534022402974906E-11 1.599233043351777E-11 0.000000000000000E+00 0.000000000000000E+00 2.781128992821719E-10 -3.037767356028012E-10 0.000000000000000E+00 0.000000000000000E+00 9.302826502045416E-10 -1.053755210769576E-09 0.000000000000000E+00 0.000000000000000E+00 1.024407175377766E-09 -1.247091348714154E-09 0.000000000000000E+00 0.000000000000000E+00 -1.159779446699762E-09 9.173686375935570E-10 0.000000000000000E+00 0.000000000000000E+00 -9.551827807144180E-09 8.945305774231066E-09 0.000000000000000E+00 0.000000000000000E+00 -4.780515429380643E-08 4.342688629398271E-08 0.000000000000000E+00 0.000000000000000E+00 -3.604407231982447E-07 3.913690673654217E-07 0.000000000000000E+00 0.000000000000000E+00
+ 16 0 3.474764301130512E-11 -3.810788024039506E-11 0.000000000000000E+00 0.000000000000000E+00 1.081569525299195E-10 -1.178623678559804E-10 0.000000000000000E+00 0.000000000000000E+00 1.507815670091551E-11 2.814088982465353E-11 0.000000000000000E+00 0.000000000000000E+00 -3.373380636877571E-10 4.177280245818081E-10 0.000000000000000E+00 0.000000000000000E+00 -5.818597193311137E-10 6.264161364052450E-10 0.000000000000000E+00 0.000000000000000E+00 1.846449434927652E-09 -1.994667620876201E-09 0.000000000000000E+00 0.000000000000000E+00 1.915570495705607E-08 -2.095831764960398E-08 0.000000000000000E+00 0.000000000000000E+00 7.584655355530582E-08 -1.402111970447130E-07 0.000000000000000E+00 0.000000000000000E+00
+ 17 0 1.085911184381699E-12 -1.194197673505632E-12 0.000000000000000E+00 0.000000000000000E+00 -4.071545074096268E-11 4.629849603642142E-11 0.000000000000000E+00 0.000000000000000E+00 -1.435817120584105E-10 1.419470173136273E-10 0.000000000000000E+00 0.000000000000000E+00 -1.308998216717278E-10 1.309323635935785E-10 0.000000000000000E+00 0.000000000000000E+00 2.153802160607438E-10 -1.849239206654281E-10 0.000000000000000E+00 0.000000000000000E+00 9.985011283015466E-10 -8.484645725242429E-10 0.000000000000000E+00 0.000000000000000E+00 2.889674418127198E-09 -1.698521099248319E-09 0.000000000000000E+00 0.000000000000000E+00 8.002500606685238E-08 -6.474815340032061E-08 0.000000000000000E+00 0.000000000000000E+00
+ 18 0 -4.787020056125859E-12 5.235296022469060E-12 0.000000000000000E+00 0.000000000000000E+00 -1.296108191380588E-11 1.306060343442856E-11 0.000000000000000E+00 0.000000000000000E+00 2.561095225761027E-11 -2.439457103666858E-11 0.000000000000000E+00 0.000000000000000E+00 6.686623770592315E-11 -8.347857886082847E-11 0.000000000000000E+00 0.000000000000000E+00 2.068860208818236E-11 -4.294448142234546E-11 0.000000000000000E+00 0.000000000000000E+00 -5.723363694873955E-10 5.256664195074986E-10 0.000000000000000E+00 0.000000000000000E+00 -4.399065911524824E-09 4.040825702730734E-09 0.000000000000000E+00 0.000000000000000E+00 -5.045745254025505E-08 5.615707341427055E-08 0.000000000000000E+00 0.000000000000000E+00
+ 19 0 2.476873300596411E-14 -6.223707320797854E-14 0.000000000000000E+00 0.000000000000000E+00 7.256937541235760E-12 -7.831135456103413E-12 0.000000000000000E+00 0.000000000000000E+00 1.387183155980786E-11 -1.564020180383458E-11 0.000000000000000E+00 0.000000000000000E+00 1.280596473101765E-11 -6.408190552182444E-12 0.000000000000000E+00 0.000000000000000E+00 -8.841120052404723E-12 1.867415588949782E-11 0.000000000000000E+00 0.000000000000000E+00 4.266686110756091E-11 -4.463973345344577E-11 0.000000000000000E+00 0.000000000000000E+00 8.990013736183562E-10 -9.002312259807897E-10 0.000000000000000E+00 0.000000000000000E+00 -4.742044039673486E-09 3.180905965801948E-10 0.000000000000000E+00 0.000000000000000E+00
+ 20 0 6.585982586468189E-13 -6.952749427988881E-13 0.000000000000000E+00 0.000000000000000E+00 9.773666173224589E-13 -1.019705638111360E-12 0.000000000000000E+00 0.000000000000000E+00 -4.100262881941560E-12 4.594781442457957E-12 0.000000000000000E+00 0.000000000000000E+00 -1.703896688112470E-11 1.609590101326194E-11 0.000000000000000E+00 0.000000000000000E+00 -9.689499420886095E-12 9.660811560169936E-12 0.000000000000000E+00 0.000000000000000E+00 5.421046421575961E-11 -4.729985940927915E-11 0.000000000000000E+00 0.000000000000000E+00 4.153224248180884E-10 -3.597092546498331E-10 0.000000000000000E+00 0.000000000000000E+00 1.365790513800586E-08 -1.327752900952558E-08 0.000000000000000E+00 0.000000000000000E+00
+ 21 0 -2.952884055119933E-14 3.617645261095944E-14 0.000000000000000E+00 0.000000000000000E+00 -1.088295533240058E-12 1.152503172466951E-12 0.000000000000000E+00 0.000000000000000E+00 -1.723369243598130E-12 1.852826615044261E-12 0.000000000000000E+00 0.000000000000000E+00 2.945564171797096E-12 -2.872803287257627E-12 0.000000000000000E+00 0.000000000000000E+00 5.231214672257249E-12 -7.848294636687318E-12 0.000000000000000E+00 0.000000000000000E+00 -1.916319688025303E-11 1.247640975825311E-11 0.000000000000000E+00 0.000000000000000E+00 -2.645440681019831E-10 2.374658095235364E-10 0.000000000000000E+00 0.000000000000000E+00 -2.491780260235909E-09 3.226308626500275E-09 0.000000000000000E+00 0.000000000000000E+00
+ 22 0 -8.929242006219457E-14 9.278028236399091E-14 0.000000000000000E+00 0.000000000000000E+00 -6.014484856522291E-14 7.317480931896816E-14 0.000000000000000E+00 0.000000000000000E+00 7.617151702519985E-13 -8.782889354051124E-13 0.000000000000000E+00 0.000000000000000E+00 1.543489778992982E-12 -1.494243261708172E-12 0.000000000000000E+00 0.000000000000000E+00 9.505692029839608E-13 4.982624452159275E-13 0.000000000000000E+00 0.000000000000000E+00 -1.113617606456332E-12 5.709544796467534E-12 0.000000000000000E+00 0.000000000000000E+00 2.204424974949836E-11 -1.092846679253227E-11 0.000000000000000E+00 0.000000000000000E+00 -2.258175910797723E-09 2.085280130938560E-09 0.000000000000000E+00 0.000000000000000E+00
+ 23 0 9.065251905033686E-15 -1.057018252700562E-14 0.000000000000000E+00 0.000000000000000E+00 1.538207946190563E-13 -1.736438587689405E-13 0.000000000000000E+00 0.000000000000000E+00 2.050907381771085E-13 -1.972069459151243E-13 0.000000000000000E+00 0.000000000000000E+00 -6.252819283697561E-13 5.161366095121224E-13 0.000000000000000E+00 0.000000000000000E+00 -1.899643748393911E-12 1.400202703373992E-12 0.000000000000000E+00 0.000000000000000E+00 1.194360232669751E-12 -2.713320128957647E-12 0.000000000000000E+00 0.000000000000000E+00 2.886471198826756E-11 -3.247311795762251E-11 0.000000000000000E+00 0.000000000000000E+00 9.035538509789723E-10 -1.103246648737506E-09 0.000000000000000E+00 0.000000000000000E+00
+ 24 0 1.174612960756333E-14 -1.126844825722886E-14 0.000000000000000E+00 0.000000000000000E+00 6.428471189592457E-15 4.121195061811011E-15 0.000000000000000E+00 0.000000000000000E+00 -1.550527038455116E-13 1.780269818272596E-13 0.000000000000000E+00 0.000000000000000E+00 -6.245661410852694E-14 1.608349533184838E-13 0.000000000000000E+00 0.000000000000000E+00 6.170481470025763E-13 -4.975726505732894E-13 0.000000000000000E+00 0.000000000000000E+00 8.542320932582662E-13 -8.356487722371056E-13 0.000000000000000E+00 0.000000000000000E+00 -8.981862246110923E-12 7.897137988346914E-12 0.000000000000000E+00 0.000000000000000E+00 2.521910786093891E-10 -1.931113582046620E-10 0.000000000000000E+00 0.000000000000000E+00
+ 25 0 -2.166387071713295E-15 1.972254371917784E-15 0.000000000000000E+00 0.000000000000000E+00 -2.363246027753549E-14 2.910482012503409E-14 0.000000000000000E+00 0.000000000000000E+00 9.306965958091072E-16 1.079709189070911E-14 0.000000000000000E+00 0.000000000000000E+00 5.945201662935406E-14 -6.932799626674669E-14 0.000000000000000E+00 0.000000000000000E+00 9.600159976192443E-14 -3.624079104633833E-14 0.000000000000000E+00 0.000000000000000E+00 -3.816807755619390E-13 6.541144083591703E-13 0.000000000000000E+00 0.000000000000000E+00 -2.142468931378792E-12 4.193529029515934E-12 0.000000000000000E+00 0.000000000000000E+00 -1.544586175302791E-10 2.366410258976099E-10 0.000000000000000E+00 0.000000000000000E+00
+ 26 0 3.241830658928220E-16 2.281988194284258E-15 0.000000000000000E+00 0.000000000000000E+00 1.426111066501638E-15 -8.492089249742721E-15 0.000000000000000E+00 0.000000000000000E+00 2.809085453347072E-14 -3.472509833006675E-14 0.000000000000000E+00 0.000000000000000E+00 5.261241409916813E-14 -2.751340745814023E-14 0.000000000000000E+00 0.000000000000000E+00 -4.477741587415061E-14 1.005572928496178E-13 0.000000000000000E+00 0.000000000000000E+00 -1.297097421854425E-13 1.937595907801317E-13 0.000000000000000E+00 0.000000000000000E+00 1.124654201489716E-12 -1.963423999666651E-12 0.000000000000000E+00 0.000000000000000E+00 -3.364106577785311E-11 -3.467444317281823E-12 0.000000000000000E+00 0.000000000000000E+00
+ 27 0 -5.874688808145129E-16 -7.860823976910845E-16 0.000000000000000E+00 0.000000000000000E+00 -1.843673154696316E-14 -1.970052499723189E-14 0.000000000000000E+00 0.000000000000000E+00 -6.780510092318230E-14 -4.510776111366325E-14 0.000000000000000E+00 0.000000000000000E+00 -1.595061702913071E-13 -7.701016144066863E-14 0.000000000000000E+00 0.000000000000000E+00 -2.206245467069110E-13 -1.977715804744331E-13 0.000000000000000E+00 0.000000000000000E+00 -2.168467170920278E-13 -5.496836829004818E-13 0.000000000000000E+00 0.000000000000000E+00 -5.951762309309079E-14 -1.301954290404940E-12 0.000000000000000E+00 0.000000000000000E+00 5.002585191818008E-12 -3.355899105022525E-11 0.000000000000000E+00 0.000000000000000E+00
+ 28 0 8.753291649675167E-15 4.417457468014839E-15 0.000000000000000E+00 0.000000000000000E+00 1.167249812200004E-13 9.027937584224611E-14 0.000000000000000E+00 0.000000000000000E+00 2.339826587015167E-13 1.824106774113203E-13 0.000000000000000E+00 0.000000000000000E+00 3.842626749737954E-13 2.779093957170420E-13 0.000000000000000E+00 0.000000000000000E+00 5.859698103369001E-13 4.047272234312309E-13 0.000000000000000E+00 0.000000000000000E+00 8.237217458372439E-13 6.698263902450226E-13 0.000000000000000E+00 0.000000000000000E+00 7.782154803264459E-14 1.252818713871486E-12 0.000000000000000E+00 0.000000000000000E+00 1.651166828997631E-11 5.085494430080715E-12 0.000000000000000E+00 0.000000000000000E+00
+ 29 0 -5.178014471779466E-16 6.239003975371151E-16 0.000000000000000E+00 0.000000000000000E+00 2.362443967581400E-14 3.530765776663440E-14 0.000000000000000E+00 0.000000000000000E+00 8.353414011476994E-14 1.104244001141589E-13 0.000000000000000E+00 0.000000000000000E+00 1.940312039993623E-13 2.434251246135747E-13 0.000000000000000E+00 0.000000000000000E+00 3.836392842190700E-13 4.783134084706228E-13 0.000000000000000E+00 0.000000000000000E+00 6.983505760143940E-13 9.039163236984156E-13 0.000000000000000E+00 0.000000000000000E+00 1.151173606433346E-12 1.503735146668129E-12 0.000000000000000E+00 0.000000000000000E+00 1.360704678907421E-11 2.322075285404156E-12 0.000000000000000E+00 0.000000000000000E+00
+ 30 0 1.873544924725119E-14 5.214628508722986E-15 0.000000000000000E+00 0.000000000000000E+00 -1.343907098461939E-13 -2.321594740300153E-13 0.000000000000000E+00 0.000000000000000E+00 -2.865409504240518E-13 -4.859507519109229E-13 0.000000000000000E+00 0.000000000000000E+00 -4.324437598797033E-13 -7.688210320926267E-13 0.000000000000000E+00 0.000000000000000E+00 -5.505184994220200E-13 -1.107023122409338E-12 0.000000000000000E+00 0.000000000000000E+00 -5.975418909147068E-13 -1.552738912828717E-12 0.000000000000000E+00 0.000000000000000E+00 -1.402929132145298E-12 -2.309899860505471E-12 0.000000000000000E+00 0.000000000000000E+00 -2.041475056207957E-11 -8.039344124784536E-12 0.000000000000000E+00 0.000000000000000E+00
+ 31 0 -2.984437747709454E-16 2.942273500679816E-16 0.000000000000000E+00 0.000000000000000E+00 -1.831272809928283E-14 7.595538604007984E-15 0.000000000000000E+00 0.000000000000000E+00 -5.056293551476471E-14 1.938222013421069E-14 0.000000000000000E+00 0.000000000000000E+00 -9.345973819600038E-14 3.270441943304734E-14 0.000000000000000E+00 0.000000000000000E+00 -1.389717069072252E-13 4.034984081233282E-14 0.000000000000000E+00 0.000000000000000E+00 -1.823359078228813E-13 3.306895226949764E-14 0.000000000000000E+00 0.000000000000000E+00 -4.192326743942745E-13 1.459389997661594E-13 0.000000000000000E+00 0.000000000000000E+00 -7.244505903823517E-12 4.432435918368827E-13 0.000000000000000E+00 0.000000000000000E+00
+ 32 0 1.020254460643745E-13 5.891011588199194E-14 0.000000000000000E+00 0.000000000000000E+00 3.017405694420945E-12 2.423937654357218E-12 0.000000000000000E+00 0.000000000000000E+00 6.039690796171587E-12 4.865016592549510E-12 0.000000000000000E+00 0.000000000000000E+00 9.347184400495976E-12 7.497972577921347E-12 0.000000000000000E+00 0.000000000000000E+00 1.322047865753964E-11 1.050309156002523E-11 0.000000000000000E+00 0.000000000000000E+00 1.813047176492616E-11 1.417697917545544E-11 0.000000000000000E+00 0.000000000000000E+00 2.467404888814358E-11 1.939861363640770E-11 0.000000000000000E+00 0.000000000000000E+00 3.134916792674994E-11 2.843559466998492E-11 0.000000000000000E+00 0.000000000000000E+00
diff --git a/ci/G3V08L3FrGrid/compare.h5 b/ci/G3V08L3FrGrid/compare.h5
new file mode 100644
index 00000000..b6664412
Binary files /dev/null and b/ci/G3V08L3FrGrid/compare.h5 differ
diff --git a/spec_refs.bib b/spec_refs.bib
index 3b2d6d06..54a12e45 100644
--- a/spec_refs.bib
+++ b/spec_refs.bib
@@ -331,6 +331,21 @@ @article{y2017_loizu
language = "en"
}
+@article{y2020_hudson,
+ doi = {10.1088/1361-6587/ab9a61},
+ url = {https://dx.doi.org/10.1088/1361-6587/ab9a61},
+ year = {2020},
+ month = {jul},
+ publisher = {IOP Publishing},
+ volume = {62},
+ number = {8},
+ pages = {084002},
+ author = {S R Hudson and J Loizu and C Zhu and Z S Qu and C Nührenberg and S Lazerson and C B Smiet and M J Hole},
+ title = {Free-boundary MRxMHD equilibrium calculations using the stepped-pressure equilibrium code},
+ journal = {Plasma Physics and Controlled Fusion},
+ abstract = {The stepped-pressure equilibrium code (SPEC) (Hudson et al 2012 Phys. Plasmas 19, 112 502) is extended to enable free-boundary multi-region relaxed magnetohydrodynamic (MRxMHD) equilibrium calculations. The vacuum field surrounding the plasma inside an arbitrary 'computational boundary', is computed, and the virtual-casing principle is used iteratively to compute the normal field on so that the equilibrium is consistent with an externally produced magnetic field. Recent modifications to SPEC are described, such as the use of Chebyshev polynomials to describe the radial dependence of the magnetic vector potential, and a variety of free-boundary verification calculations are presented.}
+}
+
@article{y2023_baillod,
title = "{Equilibrium\textit{$\beta$}-limits} dependence on bootstrap
current in classical stellarators",
diff --git a/src/bnorml.f90 b/src/bnorml.f90
index 4a47e8fb..4aaffbc1 100644
--- a/src/bnorml.f90
+++ b/src/bnorml.f90
@@ -31,7 +31,8 @@
!l tex (Note that the computational boundary does not change, so this needs only to be determined once.)
!l tex \item At each point on the computational boundary (i.e., on the discrete grid),
!l tex \link{casing} is used to compute the plasma field using the virtual casing principle.
-!l tex I think that \link{casing} returns the field in Cartesian coordinates, i.e., ${\bf B} = B_x {\bf i} + B_y {\bf j} + B_z {\bf k}$.
+!l tex \link{casing} returns the normal field ${\bf B} \cdot {\bf n}$ on the computational boundary, with ${\bf n}$ being the unit normal vector to the
+!l tex computational boundary. Internally it computes the field in Cartesian coordinates, i.e., ${\bf B} = B_x {\bf i} + B_y {\bf j} + B_z {\bf k}$
!l tex \item In toroidal geometry, the vector transformation from Cartesian to cylindrical is given by
!l tex \be \begin{array}{cccccccccccccccccccccc}
!l tex B^R & = & & + B_x \cos \z & + & B_y \sin \z & & \\
@@ -51,22 +52,19 @@
!> \ingroup grp_free-boundary
!>
!> **free-boundary constraint**
-!>
+!>
!> - The normal field at the computational boundary, \f$\partial {\cal D}\f$, should be equal to
!> \f$\left({\bf B}_P + {\bf B}_C\right)\cdot {\bf e}_\theta \times {\bf e}_\zeta\f$,
!> where \f${\bf B}_P\f$ is the "plasma" field (produced by internal plasma currents) and is computed using virtual casing,
!> and \f${\bf B}_C\f$ is the "vacuum" field (produced by the external coils) and is given on input.
!> - The plasma field, \f${\bf B}_P\f$, can only be computed after the equilibrium is determined,
-!> but this information is required to compute the equilibrium to begin with; and so there is an iteration involved.
+!> but this information is required to compute the equilibrium to begin with; and so there is an iteration involved \cite y2020_hudson .
!> - Suggested values of the vacuum field can be self generated; see xspech() for more documentation on this.
!>
!>
!> **compute the normal field on a regular grid on the computational boundary**
!>
!> - For each point on the compuational boundary, casing() is called to compute the normal field produced by the plasma currents.
-!> - \todo There is a very clumsy attempt to parallelize this which could be greatly improved.
-!>
-!>
!> - An FFT gives the required Fourier harmonics.
!>
!> \see casing.f90
@@ -81,11 +79,11 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
use constants, only : zero, half, one, two, pi, pi2, ten
- use numerical, only : small
+ use numerical, only : small, vsmall
use fileunits, only : ounit, lunit
- use inputlist, only : Wmacros, Wbnorml, Igeometry, Lcheck, vcasingtol, vcasingper, Lrad
+ use inputlist, only : Wmacros, Wbnorml, Igeometry, Lcheck, vcasingtol, vcasingits, vcasingper, Lrad, Lvcgrid, vcNz, vcNt, Nfp
use cputiming, only : Tbnorml
@@ -96,7 +94,8 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
im, in, Ate, Aze, Ato, Azo, &
Nt, Nz, cfmn, sfmn, &
ijreal, ijimag, jireal, jiimag, &
- globaljk, tetazeta, virtualcasingfactor, gteta, gzeta, Dxyz, Nxyz
+ globaljk, virtualcasingfactor, gteta, gzeta, Dxyz, Nxyz, &
+ Jxyz, Pbxyz, YESstellsym, prevcstride
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
@@ -105,12 +104,12 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
INTEGER, intent(in) :: mn, Ntz
REAL , intent(out) :: efmn(1:mn), ofmn(1:mn)
- INTEGER :: lvol, Lcurvature, Lparallel, ii, jj, kk, jk, ll, kkmodnp, jkmodnp, ifail, id01daf, nvccalls, icasing, ideriv
- REAL :: lss, zeta, teta, cszeta(0:1), tetalow, tetaupp, absacc, gBn
- REAL :: Jxyz(1:Ntz,1:3), Bxyz(1:Ntz,1:3), dAt(1:Ntz), dAz(1:Ntz), distance(1:Ntz)
+ INTEGER :: lvol, Lparallel, ii, jj, kk, jk, ll, kkmodnp, jkmodnp, ifail, id01daf, nvccalls, icasing
+ REAL :: zeta, teta, gBn
+ REAL :: Bxyz(1:Ntz,1:3), distance(1:Ntz)
+ REAL :: absvcerr, relvcerr, resulth, resulth2, resulth4, deltah4h2, deltah2h
+ INTEGER :: vcstride, Nzwithsym, vcNznfp
- !REAL :: vcintegrand, zetalow, zetaupp
-! external :: vcintegrand, zetalow, zetaupp
BEGIN(bnorml)
@@ -121,11 +120,9 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
ijreal(1:Ntz) = zero ! normal plasma field; 15 Oct 12;
- !ijimag(1:Ntz) = zero
-
- !jireal(1:Ntz) = zero
- !jiimag(1:Ntz) = zero
+ ijimag(1:Ntz) = zero
+ vcNznfp = vcNz * nfp
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
#ifdef DEBUG
@@ -137,13 +134,128 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
#endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+ ! In stellerator symmetric geometry, the virtual casing field also exhibits symmetriy which we can exploit
+ if (YESstellsym) then
+ ! we have to compute approximately half points due to symmetry, but including the middle and the edge points! (Closed interval [0, pi])
+ Nzwithsym = min(Nz, (Nz + 3) / 2 )
+ else
+ Nzwithsym = Nz
+ endif
+
+#ifdef COMPARECASING
+! When comparing results, both methods should always run
+if (.true.) then
+#else
+if ( Lvcgrid.eq.1 ) then
+#endif
+ ! Precompute Jxyz(:,1:3) and the corresponding positions on the high resolution plasma boundary
+ Pbxyz(:,1:3) = zero
+ Jxyz(:,1:3) = zero
+ !$OMP PARALLEL DO SHARED(Pbxyz, Jxyz) PRIVATE(jk, jj, kk, teta, zeta) COLLAPSE(2)
+ do kk = 0, vcNznfp-1 ;
+ do jj = 0, vcNt-1 ;
+ zeta = kk * pi2 / vcNznfp
+ teta = jj * pi2 / vcNt ;
+ jk = 1 + jj + kk*vcNt
+
+ ! Each MPI rank only computes every a 1/ncpu surfacecurrent() calls
+ select case( Lparallel )
+ case( 0 ) ! Lparallel = 0
+ if( myid.ne.modulo(kk,ncpu) ) cycle
+ case( 1 ) ! Lparallel = 1
+ if( myid.ne.modulo(jk-1,ncpu) ) cycle
+ case default ! Lparallel
+ FATAL( bnorml, .true., invalid Lparallel in parallelization loop )
+ end select
+
+ ! pbxyz and jxyz are both [out] parameters
+ call surfacecurrent( teta, zeta, Pbxyz(jk,1:3), Jxyz(jk,1:3) )
+ enddo
+ enddo
+
+ ! MPI reductions for positions and currents to accumulate them on all ranks (valid because initialized to zero)
+ ! and Broadcast the total currents and evaluation points back to all ranks
+ call MPI_Allreduce(MPI_IN_PLACE, Pbxyz(:,1:3), 3*vcNt*vcNznfp, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_SPEC, ierr )
+ if (ierr.ne.MPI_SUCCESS) then
+ FATAL( bnorml, .true., error in MPI_Allreduce for Pbxyz )
+ endif
+ call MPI_Allreduce(MPI_IN_PLACE, Jxyz(:,1:3), 3*vcNt*vcNznfp, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_SPEC, ierr )
+ if (ierr.ne.MPI_SUCCESS) then
+ FATAL( bnorml, .true., error in MPI_Allreduce for Jxyz )
+ endif
+
+ ! iterate over resolutions of the virtual casing grid to get an estimate of the accuracy. Write the result into ijimag
+ deltah2h = 1e12
+ !> The surface currents get precomputed at all points, but maybe we don't have to integrate over the whole grid.
+ !> Start with a large stride over the plasmaboundary (= how many points to skip at each step), then get progressively finer.
+ !> Do at least vcasingits iterations (sqrt(vcasingits) resolution per dimension) and halve the stride until the accuracy is reached.
+ prevcstride = min(prevcstride, INT(0.5*log(real(vcNt*vcNznfp/vcasingits))/log(2.0)))
+ prevcstride = max(prevcstride, 1) !> At least one step in the resolution cascade is required to determine an estimate of the current accuracy.
+ do vcstride = prevcstride, 0, -1
+ !$OMP PARALLEL DO SHARED(Dxyz, Nxyz, Pbxyz, Jxyz, ijreal, ijimag) FIRSTPRIVATE(jk, gBn) COLLAPSE(2)
+ do kk = 0, Nzwithsym-1 ;
+ do jj = 0, Nt-1 ;
+ jk = 1 + jj + kk*Nt
+
+ ! kk = 0 terms are also symmetryic (relevant e.g. for tokamaks)
+ if ((Nz.eq.1) .and. (jj.gt.((Nt+1)/2))) cycle
+
+ ! Each MPI rank only computes every a 1/ncpu surfacecurrent() calls
+ ! Identical MPI parallelization scheme as for Lvcgrid=0
+ select case( Lparallel )
+ case( 0 ) ! Lparallel = 0
+ if( myid.ne.modulo(kk,ncpu) ) cycle
+ case( 1 ) ! Lparallel = 1
+ if( myid.ne.modulo(jk-1,ncpu) ) cycle
+ case default ! Lparallel;
+ FATAL( bnorml, .true., invalid Lparallel in parallelization loop )
+ end select ! end of select case( Lparallel )
+
+ gBn = zero
+ globaljk = jk ! only to check against the precomputed values against dvcfield
+ call casinggrid( Dxyz(1:3,jk), Nxyz(1:3,jk), Pbxyz, Jxyz, 2**vcstride, gBn)
+
+ ijimag(jk) = ijreal(jk) ! previous solution (lower resolution)
+ ijreal(jk) = gBn ! current solution (higher resolution)
+ enddo ! end of do jj
+ enddo ! end of do kk
+
+ deltah4h2 = deltah2h
+ deltah2h = maxval(abs(ijimag - ijreal)) ! mean delta between the h and h/2 solutions
+
+ ! Order of the integration method: log(deltah4h2/deltah2h)/log(2.0) = 1
+ absvcerr = deltah2h
+ relvcerr = 2 * deltah2h / maxval(abs(ijreal)) ! overestimate the relative error by a factor of two
+ if (myid.eq.0) then
+ write(ounit, '("bnorml : ", 10x ," : relvcerr = ",es13.5," ; absvcerr = ",es13.5," ; resolution reduction = ",i8)') relvcerr, absvcerr, 2**vcstride
+ ! print *, "Convergence order: ", log(deltah4h2/deltah2h)/log(2.0)
+ endif
+ ! Tolerance was already reached, exit the loop. (Different threads may exit at different times)
+ if (absvcerr.lt.vcasingtol .and. relvcerr.lt.vcasingtol) then
+ prevcstride = vcstride
+ exit
+ endif
+ enddo ! end of do vcstride
+#ifdef COMPARECASING
+ ! To compare with the casing() implementation, copy the results into ijimag
+ if(Lvcgrid.eq.1) then
+ ijimag(1:Ntz) = ijreal(1:Ntz)
+ endif
+endif
+! When comparing results, both methods should always run
+if (.true.) then
+#else
+else ! if not Lvcgrid
+#endif
- do kk = 0, Nz-1 ; zeta = kk * pi2nfp / Nz
+ do kk = 0, Nzwithsym-1 ;
+ zeta = kk * pi2 / Nz
- if( Igeometry.eq.3 ) then ; cszeta(0:1) = (/ cos(zeta), sin(zeta) /)
- endif
+ do jj = 0, Nt-1 ;
+ teta = jj * pi2 / Nt ; jk = 1 + jj + kk*Nt
- do jj = 0, Nt-1 ; teta = jj * pi2 / Nt ; jk = 1 + jj + kk*Nt
+ ! kk = 0 terms are also symmetryic (relevant e.g. for tokamaks)
+ if ((Nz.eq.1) .and. (jj.gt.((Nt+1)/2))) cycle
globaljk = jk ! this is global; passed through to vcintegrand & casing;
@@ -156,35 +268,18 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
FATAL( bnorml, .true., invalid Lparallel in parallelization loop )
end select ! end of select case( Lparallel ) ; 09 Mar 17;
- tetazeta(1:2) = (/ teta, zeta /) ! this is global; passed through to zetalow & zetaupp; 14 Apr 17;
-
-!#ifdef COMPARECASING
-!
-! tetalow = tetazeta(1) - vcasingper * pi ; tetaupp = tetazeta(1) + vcasingper * pi ; absacc = vcasingtol
-!
-! id01daf = 1 ; call D01DAF( tetalow, tetaupp, zetalow, zetaupp, vcintegrand, absacc, gBn, nvccalls, id01daf ) ! 04 May 17;
-!
-! ijimag(jk) = gBn
-!
-!#endif
-
- WCALL( bnorml, casing, ( teta, zeta, gBn, icasing ) ) ! tetazeta is global; 26 Apr 17;
-
+ WCALL( bnorml, casing, ( teta, zeta, gBn, icasing ) )
ijreal(jk) = gBn
-!#ifdef COMPARECASING
-! write(ounit,1000) myid, zeta, teta, ijreal(jk), ijimag(jk), ijreal(jk)-ijimag(jk)
-!#endif
-!
-!#ifdef CASING
-! write(ounit,1000) myid, zeta, teta, ijreal(jk), ijimag(jk), ijreal(jk)-ijimag(jk)
-!#endif
-
-1000 format("bnorml : ", 10x ," : myid=",i3," : \z =",f6.3," ; \t =",f6.3," ; B . x_t x x_z =",2f22.15," ; ":"err =",es13.5," ;")
+#ifdef COMPARECASING
+ ! write(ounit,1000) myid, zeta, teta, ijreal(jk), ijimag(jk), ijreal(jk)-ijimag(jk)
+ write(ounit,1000) myid, zeta, teta, ijreal(jk), ijimag(jk), ijreal(jk)-ijimag(jk)
+ 1000 format("bnorml : ", 10x ," : myid=",i3," : \z =",f6.3," ; \t =",f6.3," ; B . x_t x x_z =",2f22.15," ; ":"err =",es13.5," ;")
+#endif
enddo ! end of do jj;
enddo ! end of do kk;
-
+endif ! end of if (Lvcgrid )
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
1001 format("bnorml : ", 10x ," : "a1" : (t,z) = ("f8.4","f8.4" ) ; gBn=",f23.15," ; ":" error =",f23.15" ;")
@@ -206,11 +301,6 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
case( 0 ) ! Lparallel = 0 ; 09 Mar 17;
RlBCAST(ijreal(1+kk*Nt:Nt+kk*Nt),Nt,kkmodnp) ! plasma; 03 Apr 13;
- !RlBCAST(ijimag(1+kk*Nt:Nt+kk*Nt),Nt,kkmodnp)
-
- !RlBCAST(jireal(1+kk*Nt:Nt+kk*Nt),Nt,kkmodnp)
- !RlBCAST(jiimag(1+kk*Nt:Nt+kk*Nt),Nt,kkmodnp)
-
case( 1 ) ! Lparallel = 1 ; 09 Mar 17;
do jj = 0, Nt-1
@@ -220,10 +310,6 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
jkmodnp = modulo(jk-1,ncpu)
RlBCAST(ijreal(jk),1,jkmodnp) ! plasma; 03 Apr 13;
- !RlBCAST(ijimag(jk),1,jkmodnp)
-
- !RlBCAST(jireal(jk),1,jkmodnp)
- !RlBCAST(jiimag(jk),1,jkmodnp)
enddo
@@ -235,6 +321,25 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
enddo ! 11 Oct 12;
+!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+ ! In stellerator symmetric geometry, the virtual casing field also exhibits symmetriy which we can exploit
+ if (YESstellsym) then
+ ! We are dealing with half open intervals in theta and phi [0, 2pi[ so we need to skip the first row and column when mirroring
+ if (Nz.eq.1) then
+ do jj = 0, Nt/2;
+ jk = 1 + jj
+ ijreal(Ntz + 1 - jk) = -ijreal(jk+1)
+ enddo
+ endif
+
+ do kk = 0, Nzwithsym-2;
+ do jj = 0, Nt-1;
+ jk = 1 + jj + kk*Nt
+
+ ijreal(Ntz + 1 - jk) = -ijreal(jk+Nt+1)
+ enddo
+ enddo
+ endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
ijreal(1:Ntz) = ijreal(1:Ntz) * virtualcasingfactor
diff --git a/src/casing.f90 b/src/casing.f90
index 8bfd29b0..fa9afb6f 100644
--- a/src/casing.f90
+++ b/src/casing.f90
@@ -65,37 +65,6 @@
!> The minimum number of function evaluations is provided by the \c inputvar \c vcasingits.
!>
!>
-!> **Calculation of integrand**
-!>
-!>
-!> - An adaptive integration is used to compute the integrals.
-!> Consequently, the magnetic field tangential to the plasma boundary is required at an arbitrary point.
-!> This is computed, as always, from \f${\bf B} = \nabla \times {\bf A}\f$, and this provides \f${\bf B} = B^\theta {\bf e}_\theta + B^\zeta {\bf e}_\zeta\f$.
-!> Recall that \f$B^s=0\f$ by construction on the plasma boundary.
-!> \todo It would be MUCH faster to only require the tangential field on a regular grid!!!
-!>
-!>
-!> - Then, the metric elements \f$g_{\theta\theta}\f$, \f$g_{\theta\zeta}\f$ and \f$g_{\zeta\zeta}\f$ are computed.
-!> These are used to "lower" the components of the magnetic field, \f${\bf B} = B_\theta \nabla \theta + B_\zeta \nabla \zeta\f$.
-!> \todo Please check why \f$B_s\f$ is not computed. Is it because \f$B_s \nabla s \times {\bf n} = 0\f$ ?
-!>
-!>
-!> - The distance between the "evaluate" point, \f$(X,Y,Z)\f$, and the given point on the surface, \f$(x,y,z)\f$ is computed.
-!> - If the computational boundary becomes too close to the plasma boundary, the distance is small and this causes problems for the numerics.
-!> I have tried to regularize this problem by introducing \f$\epsilon \equiv \f$\c inputvar \c vcasingeps.
-!> Let the "distance" be
-!> \f{eqnarray}{ D \equiv \sqrt{(X-x)^2 + (Y-y)^2 + (Z-Z)^2} + \epsilon^2.
-!> \f}
-!> - On taking the limit that \f$\epsilon \rightarrow 0\f$, the virtual casing integrand is
-!> \f{eqnarray}{ \verb+vcintegrand+ \equiv ( B_x n_x + B_y n_y + B_z n_z ) ( 1 + 3 \epsilon^2 / D^2 ) / D^3, \label{eq:integrand_casing}
-!> \f}
-!> where the normal vector is \f${\bf n} \equiv n_x {\bf i} + n_y {\bf j} + n_z {\bf k}\f$.
-!> The normal vector, \c Nxyz , to the computational boundary (which does not change) is computed in preset().
-!> \todo This needs to be revised.
-!>
-!>
-!>
-!>
!> @param[in] teta \f$\theta\f$
!> @param[in] zeta \f$\zeta\f$
!> @param[out] gBn \f$ \sqrt g {\bf B} \cdot {\bf n}\f$
@@ -106,8 +75,6 @@ subroutine casing( teta, zeta, gBn, icasing )
use constants, only : zero, pi, pi2
- use numerical, only :
-
use fileunits, only : ounit, vunit
use inputlist, only : Wmacros, Wcasing, vcasingtol, vcasingits, vcasingper
@@ -120,13 +87,13 @@ subroutine casing( teta, zeta, gBn, icasing )
LOCALS
- REAL, intent(in) :: teta, zeta ! arbitrary location; Cartesian;
- REAL, intent(out) :: gBn ! magnetic field; Cartesian;
+ REAL, intent(in) :: teta, zeta ! theta and zeta evaluation point on the computational boundary, corresponding to the cartesian evaluation point Dxyz(1:3,globaljk)
+ REAL, intent(out) :: gBn ! B.n magnetic field
INTEGER, intent(out) :: icasing
INTEGER, parameter :: Ndim = 2, Nfun = 1
- INTEGER :: ldim, lfun, minpts, maxpts, Lrwk, idcuhre, jk, irestart, funcls, key, num, maxsub
+ INTEGER :: ldim, lfun, minpts, maxpts, Lrwk, idcuhre, irestart, funcls, key, num, maxsub
REAL :: integrals(1:Nfun), low(1:Ndim), upp(1:Ndim), labs, lrel, absest(1:Nfun)
REAL, allocatable :: rwk(:)
@@ -135,9 +102,6 @@ subroutine casing( teta, zeta, gBn, icasing )
BEGIN(casing)
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
- jk = globaljk ! shorthand; globaljk is a "global" variable which must be passed through to subroutine dvcfield;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
@@ -185,31 +149,31 @@ subroutine casing( teta, zeta, gBn, icasing )
select case( idcuhre ) ! "123456789012345678901234"
case(0) ;
; ; exit
- case(1) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "maxpts too smal; "
+ case(1) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "maxpts too smal; "
; ;!exit
- case(2) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "illegal key; "
+ case(2) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "illegal key; "
; ; exit
- case(3) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "illegal Ndim; "
+ case(3) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "illegal Ndim; "
; ; exit
- case(4) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "key.eq.1 & Ndim.ne.2; "
+ case(4) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "key.eq.1 & Ndim.ne.2; "
; ; exit
- case(5) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "key.eq.2 & Ndim.ne.3; "
+ case(5) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "key.eq.2 & Ndim.ne.3; "
; ; exit
- case(6) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "numfun < 1; "
+ case(6) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "numfun < 1; "
; ; exit
- case(7) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "volume is zero; "
+ case(7) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "volume is zero; "
; ; exit
- case(8) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "maxpts < 3*NUM; "
+ case(8) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "maxpts < 3*NUM; "
; ; exit
- case(9) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "maxpts < minpts; "
+ case(9) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "maxpts < minpts; "
; ; exit
- case(10) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "epsabs < 0 & epsrel < 0;"
+ case(10) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "epsabs < 0 & epsrel < 0;"
; ; exit
- case(11) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "NW is too small; "
+ case(11) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "NW is too small; "
; ; exit
- case(12) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "illegal irestart; "
+ case(12) ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "illegal irestart; "
; ; exit
- case default ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "tryin' to kill me? "
+ case default ; write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts, "tryin' to kill me? "
; ; exit
end select
@@ -227,7 +191,7 @@ subroutine casing( teta, zeta, gBn, icasing )
enddo ! end of virtual casing accuracy infinite-do-loop; 10 Apr 13;
#ifdef DEBUG
- ; ; if( Wcasing ) write(ounit,1001) cput-cpus, myid, Dxyz(1:3,jk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts
+ ; ; if( Wcasing ) write(ounit,1001) cput-cpus, myid, Dxyz(1:3,globaljk), gBn, absest(1:Nfun), idcuhre, minpts, maxpts
#endif
1001 format("casing : ",f10.2," : myid=",i3," ; [x,y,z]=["es10.2" ,"es10.2" ,"es10.2" ]; gBn="es12.4" , ",&
@@ -249,265 +213,399 @@ subroutine casing( teta, zeta, gBn, icasing )
end subroutine casing
+!> \brief Compute the field produced by the plasma currents, at an arbitrary, external location using virtual casing.
+!> \ingroup grp_free-boundary
+!>
+!> casinggrid() computes the virtual casing field using a fixed resolution grid on the plasma boundary, in contrast to casing() which uses an adaptive integration routine.
+!> Because the evaluation positions are known in advance, the most expensive terms of the computation (surfacecurrents, or in the casing() routine dvcfield) are precomputed.
+!> The integral computed by casinggrid() is
+!> \f{eqnarray}{ \int_0^{2\pi} \int_0^{2\pi} \frac{{\bf j}(\theta, \zeta) \times {\bf D}(\theta, \zeta)}{\|{\bf D}(\theta', \zeta')\|^3}+ \;\; d\theta d\zeta
+!> \f}
+!> with \f${\bf r}\f$ approximately the distance vector from an evaluation point of the current surface \f${\bf x}(\theta, \zeta)\f$ to the
+!> external point \f$(x,y,z)\f$, and the exact formula includes a regularization factor \f$\epsilon\f$ Eqn.\f$(\ref{eq:vcasing_distance})\f$.
+!>
+!> Numerically it uses a Kahan Summation algorithm to accumulate the contributions of all inner points, because truncation error quickly dominates the computation when thousands
+!> of points are used in both dimensions.
+!>
+!> @param[in] xyz \f$(x,y,z)\f$ cartesian coordinates on the computational boundary
+!> @param[in] ntz cartesian normal vector on the computational boundary
+!> @param[in] Pbxyz array of cartesian coordinates on the plasma boundary
+!> @param[in] Jxyz array of surface currents \f${\bf B}_{Plasma} \cdot {\bf e}_\theta \times {\bf e}_\zeta \;\f$ on the plasma boundary
+!> @param[in] vcstride integer stride for the fixed resolution grid. vcstride = 1 computes the integral at full resolution, vcstride = 2 computes the integral at half resolution, etc.
+!> @param[out] gBn normal field \f${\bf B}_{Plasma} \cdot n \;\f$ on the computational boundary
+subroutine casinggrid( xyz, ntz, Pbxyz, Jxyz, vcstride, gBn)
+ use constants, only : zero, one, three, pi2
+
+ use fileunits, only : ounit, vunit
+
+ use inputlist, only : vcasingeps, Nfp, vcNz, vcNt
+
+ use allglobal, only : myid, ncpu, cpus, MPI_COMM_SPEC, Dxyz, Nxyz,globaljk,pi2nfp
+
+ LOCALS
+
+ REAL, intent(in) :: xyz(1:3) ! arbitrary location; Cartesian;
+ REAL, intent(in) :: ntz(1:3) ! surface normal on the computational boundary; Cartesian;
+ REAL, intent(in) :: Pbxyz(1:vcNz*Nfp*vcNt, 1:3), Jxyz(1:vcNz*Nfp*vcNt, 1:3)
+ INTEGER, intent(in) :: vcstride
+ REAL, intent(out) :: gBn ! B.n on the computational boundary;
+
+ REAL :: rr(1:3), distance(1:3), jj(1:3), Bxyz(1:3), accumulator, firstorderfactor
+ REAL :: gBnlocal ! Local accumulator so we don't directly sum into the destination pointer (Avoid cache thrashing due to false share)
+ REAL :: c, t, y ! Helper variables for the Kahan summation
+ INTEGER :: plasmaNtz, jk
+!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+ ! plasmaNtz = SIZE(Pbxyz, 1)
+ plasmaNtz = vcNz*Nfp*vcNt
+ gBnlocal = zero
+ c = zero
+ ! loop over the high resolution plasma boundary (inner boundary for virtual casing)
+ do jk = 1, plasmaNtz, vcstride ;
+ ! distance vector between point on computational boundary and point on plasma boundary
+ rr = xyz - Pbxyz(jk, 1:3)
+ jj = Jxyz(jk, 1:3)
+
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+ distance(2) = sum( rr * rr ) + vcasingeps**2 ! 04 May 17;
+ distance(1) = sqrt( distance(2) ) ;
+ distance(3) = distance(1) * distance(2) ! powers of distance; 24 Nov 16;
+
+ firstorderfactor = ( one + three * vcasingeps**2 / distance(2) ) / distance(3) ! 04 May 17;
+
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+ Bxyz(1:3) = (/ jj(2) * rr(3) - jj(3) * rr(2), &
+ jj(3) * rr(1) - jj(1) * rr(3), &
+ jj(1) * rr(2) - jj(2) * rr(1) /)
+
+ ! Accumulate B.r/r^3 contributions
+ ! Kahan Summation algorithm
+ y = sum( Bxyz * ntz ) * firstorderfactor - c ! Correct the next value with the compensation
+ t = gBnlocal + y ! Add the corrected value to the sum
+ c = (t - gBnlocal) - y ! Update compensation
+ gBnlocal = t ! Update the running total
+ enddo
+ gBn = gBnlocal * pi2 * pi2 / (plasmaNtz / vcstride)
+ return
+
+end subroutine casinggrid
+
+
+!> \brief Compute the position and surface current on the plasma boundary for the virtual casing
+!> \ingroup grp_free-boundary
+!> **Calculation of integrand**
+!>
+!>
+!> - An adaptive integration is used to compute the integrals.
+!> Consequently, the magnetic field tangential to the plasma boundary is required at an arbitrary point.
+!> This is computed, as always, from \f${\bf B} = \nabla \times {\bf A}\f$, and this provides \f${\bf B} = B^\theta {\bf e}_\theta + B^\zeta {\bf e}_\zeta\f$.
+!> Recall that \f$B^s=0\f$ by construction on the plasma boundary.
+!>
+!>
+!> - Then, the metric elements \f$g_{\theta\theta}\f$, \f$g_{\theta\zeta}\f$ and \f$g_{\zeta\zeta}\f$ are computed.
+!> These are used to "lower" the components of the magnetic field, \f${\bf B} = B_\theta \nabla \theta + B_\zeta \nabla \zeta\f$.
+!>
+!>
+!> - The distance between the "evaluate" point, \f$(X,Y,Z)\f$, and the given point on the surface, \f$(x,y,z)\f$ is computed.
+!> - If the computational boundary becomes too close to the plasma boundary, the distance is small and this causes problems for the numerics.
+!> I have tried to regularize this problem by introducing \f$\epsilon \equiv \f$\c inputvar \c vcasingeps.
+!> Let the "distance" be
+!> \f{eqnarray}{ D \equiv \sqrt{(X-x)^2 + (Y-y)^2 + (Z-Z)^2} + \epsilon^2.
+!> \label{eq:vcasing_distance} \f}
+!> - On taking the limit that \f$\epsilon \rightarrow 0\f$, the virtual casing integrand is
+!> \f{eqnarray}{ \verb+vcintegrand+ \equiv ( B_x n_x + B_y n_y + B_z n_z ) ( 1 + 3 \epsilon^2 / D^2 ) / D^3, \label{eq:integrand_casing}
+!> \f}
+!> where the normal vector is \f${\bf n} \equiv n_x {\bf i} + n_y {\bf j} + n_z {\bf k}\f$.
+!> The normal vector, \c Nxyz , to the computational boundary (which does not change) is computed in preset().
+!> \todo This needs to be revised.
+!>
+!>
+!>
+!>
+!> @param[in] teta \f$\theta\f$
+!> @param[in] zeta \f$\zeta\f$
+!> @param[out] pxyz \f$(x,y,z)\f$ cartesian coordinates on the plasma boundary
+!> @param[out] jj Surface current \f${\bf B}_{Plasma} \cdot {\bf e}_\theta \times {\bf e}_\zeta \;\f$ on the plasma boundary
+subroutine surfacecurrent( teta, zeta, pxyz, jj)
+
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+ use constants, only : zero, half, one
+
+ use inputlist, only : Nvol, Igeometry, Lrad
+
+ use allglobal, only : myid, ncpu, cpus, MPI_COMM_SPEC, &
+ Mvol, &
+ mn, im, in, &
+ iRbc, iZbs, iRbs, iZbc, &
+ Ate, Aze, Ato, Azo, &
+ TT, &
+ YESstellsym, &
+ first_free_bound
+
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+ LOCALS
+
+ REAL , intent(in) :: teta, zeta ! theta and zeta coordinates on the plasma boundary
+ REAL , intent(out) :: pxyz(1:3) ! Position on the plasma boundary
+ REAL , intent(out) :: jj(1:3) ! Cartesian surface current components on the plasma boundary
+
+ INTEGER :: ii, mi, ni, ll, ideriv
+ REAL :: dR(0:3), dZ(0:3), gBut, gBuz, gtt, gtz, gzz, sqrtg, Blt, Blz, czeta, szeta, arg, carg, sarg
+
+ REAL :: XXt, XXz, YYt, YYz, ZZt, ZZz, ds
+
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+ dR(0:3) = zero ; ideriv = 0 ; gBut = zero ; gBuz = zero ! initialize summation of coordinates and tangential field;
+ dZ(0:3) = zero
+
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+ select case( Igeometry )
+
+ case( 1 ) ! Igeometry = 1 ; 09 Mar 17;
+
+ if( YESstellsym ) then
+
+ do ii = 1, mn ; mi = im(ii) ; ni = in(ii) ! loop over Fourier modes; construct surface current; slow transform required as position is arbitrary;
+
+ arg = mi * teta - ni * zeta ; carg = cos(arg) ; sarg = sin(arg)
+
+ dR(0) = dR(0) + ( iRbc(ii,Nvol) ) * carg
+ dR(1) = dR(1) + ( ( iRbc(ii,Mvol) - iRbc(ii,Nvol) ) * carg ) * half
+ dR(2) = dR(2) + mi * ( - ( iRbc(ii,Nvol) ) * sarg )
+ dR(3) = dR(3) - ni * ( - ( iRbc(ii,Nvol) ) * sarg )
+
+ !dZ(0) = dZ(0) + ( iZbs(ii,Nvol) ) * sarg
+ !dZ(1) = dZ(1) + ( ( iZbs(ii,Mvol) - iZbs(ii,Nvol) ) * sarg ) * half
+ !dZ(2) = dZ(2) + mi * ( ( iZbs(ii,Nvol) ) * carg )
+ !dZ(3) = dZ(3) - ni * ( ( iZbs(ii,Nvol) ) * carg )
+
+ do ll = 0, Lrad(Mvol)
+ gBut = gBut - ( Aze(Mvol,ideriv,ii)%s(ll) * carg ) * TT(ll,0,1) ! contravariant; Jacobian comes later;
+ gBuz = gBuz + ( Ate(Mvol,ideriv,ii)%s(ll) * carg ) * TT(ll,0,1)
+ enddo
+
+ enddo ! end of do ii = 1, mn ;
+
+ else ! NOTstellsym ; 08 Feb 16;
+
+ do ii = 1, mn ; mi = im(ii) ; ni = in(ii) ! loop over Fourier modes; construct surface current; slow transform required as position is arbitrary;
+
+ arg = mi * teta - ni * zeta ; carg = cos(arg) ; sarg = sin(arg)
+
+ dR(0) = dR(0) + iRbc(ii,Nvol) * carg + iRbs(ii,Nvol) * sarg
+ dR(1) = dR(1) + ( ( iRbc(ii,Mvol) - iRbc(ii,Nvol) ) * carg + ( iRbs(ii,Mvol) - iRbs(ii,Nvol) ) * sarg ) * half
+ dR(2) = dR(2) + mi * ( - iRbc(ii,Nvol) * sarg + iRbs(ii,Nvol) * carg )
+ dR(3) = dR(3) - ni * ( - iRbc(ii,Nvol) * sarg + iRbs(ii,Nvol) * carg )
+
+ !dZ(0) = dZ(0) + iZbc(ii,Nvol) * carg + iZbs(ii,Nvol) * sarg
+ !dZ(1) = dZ(1) + ( ( iZbc(ii,Mvol) - iZbc(ii,Nvol) ) * carg + ( iZbs(ii,Mvol) - iZbs(ii,Nvol) ) * sarg ) * half
+ !dZ(2) = dZ(2) + mi * ( - iZbc(ii,Nvol) * sarg + iZbs(ii,Nvol) * carg )
+ !dZ(3) = dZ(3) - ni * ( - iZbc(ii,Nvol) * sarg + iZbs(ii,Nvol) * carg )
+
+ do ll = 0, Lrad(Mvol)
+ gBut = gBut - ( Aze(Mvol,ideriv,ii)%s(ll) * carg + Azo(Mvol,ideriv,ii)%s(ll) * sarg ) * TT(ll,0,1) ! contravariant; Jacobian comes later;
+ gBuz = gBuz + ( Ate(Mvol,ideriv,ii)%s(ll) * carg + Ato(Mvol,ideriv,ii)%s(ll) * sarg ) * TT(ll,0,1)
+ enddo
+
+ enddo ! end of do ii = 1, mn ;
+
+ endif ! end of if( YESstellsym ) ; 08 Feb 16;
+
+ case( 2 ) ! Igeometry = 2 ; 09 Mar 17;
+
+ FATAL( casing, .true., virtual casing under construction for cylindrical geometry )
+
+ case( 3 ) ! Igeometry = 3 ; 09 Mar 17;
+
+ if( YESstellsym ) then
+
+ do ii = 1, mn ; mi = im(ii) ; ni = in(ii) ! loop over Fourier modes; construct surface current; slow transform required as position is arbitrary;
+
+ arg = mi * teta - ni * zeta ; carg = cos(arg) ; sarg = sin(arg)
+ dR(0) = dR(0) + ( iRbc(ii,Nvol) ) * carg
+ dR(1) = dR(1) + ( ( iRbc(ii,Mvol) - iRbc(ii,Nvol) ) * carg ) * half
+ dR(2) = dR(2) + mi * ( - ( iRbc(ii,Nvol) ) * sarg )
+ dR(3) = dR(3) - ni * ( - ( iRbc(ii,Nvol) ) * sarg )
+
+ dZ(0) = dZ(0) + ( iZbs(ii,Nvol) ) * sarg
+ dZ(1) = dZ(1) + ( ( iZbs(ii,Mvol) - iZbs(ii,Nvol) ) * sarg ) * half
+ dZ(2) = dZ(2) + mi * ( ( iZbs(ii,Nvol) ) * carg )
+ dZ(3) = dZ(3) - ni * ( ( iZbs(ii,Nvol) ) * carg )
+
+ if (first_free_bound) then
+ do ll = 0, Lrad(Nvol) ! 1 is for outside thr volume
+ gBut = gBut - ( Aze(Nvol,ideriv,ii)%s(ll) * carg ) * TT(ll,1,1) ! contravariant; Jacobian comes later;
+ gBuz = gBuz + ( Ate(Nvol,ideriv,ii)%s(ll) * carg ) * TT(ll,1,1)
+ enddo
+ else
+ do ll = 0, Lrad(Mvol)
+ gBut = gBut - ( Aze(Mvol,ideriv,ii)%s(ll) * carg ) * TT(ll,0,1) ! contravariant; Jacobian comes later;
+ gBuz = gBuz + ( Ate(Mvol,ideriv,ii)%s(ll) * carg ) * TT(ll,0,1)
+ enddo
+ endif
+
+ enddo ! end of do ii = 1, mn ;
+
+ else ! NOTstellsym ; 08 Feb 16;
+
+ do ii = 1, mn ; mi = im(ii) ; ni = in(ii) ! loop over Fourier modes; construct surface current; slow transform required as position is arbitrary;
+
+ arg = mi * teta - ni * zeta ; carg = cos(arg) ; sarg = sin(arg)
+
+ dR(0) = dR(0) + iRbc(ii,Nvol) * carg + iRbs(ii,Nvol) * sarg
+ dR(1) = dR(1) + ( ( iRbc(ii,Mvol) - iRbc(ii,Nvol) ) * carg + ( iRbs(ii,Mvol) - iRbs(ii,Nvol) ) * sarg ) * half
+ dR(2) = dR(2) + mi * ( - iRbc(ii,Nvol) * sarg + iRbs(ii,Nvol) * carg )
+ dR(3) = dR(3) - ni * ( - iRbc(ii,Nvol) * sarg + iRbs(ii,Nvol) * carg )
+
+ dZ(0) = dZ(0) + iZbc(ii,Nvol) * carg + iZbs(ii,Nvol) * sarg
+ dZ(1) = dZ(1) + ( ( iZbc(ii,Mvol) - iZbc(ii,Nvol) ) * carg + ( iZbs(ii,Mvol) - iZbs(ii,Nvol) ) * sarg ) * half
+ dZ(2) = dZ(2) + mi * ( - iZbc(ii,Nvol) * sarg + iZbs(ii,Nvol) * carg )
+ dZ(3) = dZ(3) - ni * ( - iZbc(ii,Nvol) * sarg + iZbs(ii,Nvol) * carg )
+
+ if (first_free_bound) then
+ do ll = 0, Lrad(Mvol) ! omit the possible current sheet due to a jump in tangential field at the plasma boundary; Zhu 20190603;
+ gBut = gBut - ( Aze(Mvol,ideriv,ii)%s(ll) * carg + Azo(Mvol,ideriv,ii)%s(ll) * sarg ) * TT(ll,1,1) ! contravariant; Jacobian comes later;
+ gBuz = gBuz + ( Ate(Mvol,ideriv,ii)%s(ll) * carg + Ato(Mvol,ideriv,ii)%s(ll) * sarg ) * TT(ll,1,1)
+ enddo
+ else
+ do ll = 0, Lrad(Mvol)
+ gBut = gBut - ( Aze(Mvol,ideriv,ii)%s(ll) * carg + Azo(Mvol,ideriv,ii)%s(ll) * sarg ) * TT(ll,0,1) ! contravariant; Jacobian comes later;
+ gBuz = gBuz + ( Ate(Mvol,ideriv,ii)%s(ll) * carg + Ato(Mvol,ideriv,ii)%s(ll) * sarg ) * TT(ll,0,1)
+ enddo
+ endif
+
+ enddo ! end of do ii = 1, mn ;
+
+ endif ! end of if( YESstellsym ) ; 08 Feb 16;
+
+ end select ! end of select case( Igeometry ) ; 09 Mar 17;
+
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+ select case( Igeometry )
+
+ case( 1 ) ! Igeometry = 1 ; 09 Mar 17;
+
+ gtt = one + dR(2)*dR(2)
+ gtz = dR(2)*dR(3)
+ gzz = one + dR(3)*dR(3)
+
+ sqrtg = dR(1)
+
+ case( 2 ) ! Igeometry = 2 ; 09 Mar 17;
+
+ FATAL( casing, .true., virtual casing under construction for cylindrical geometry )
+
+ case( 3 ) ! Igeometry = 3 ; 09 Mar 17;
+
+ gtt = dR(2)*dR(2) + dZ(2)*dZ(2)
+ gtz = dR(2)*dR(3) + dZ(2)*dZ(3)
+ gzz = dR(3)*dR(3) + dZ(3)*dZ(3) + dR(0)*dR(0)
+
+ sqrtg = dR(0) * ( dZ(1) * dR(2) - dR(1) * dZ(2) )
+
+ end select
+
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+ Blt = ( gBut * gtt + gBuz * gtz ) / sqrtg
+ Blz = ( gBut * gtz + gBuz * gzz ) / sqrtg
+
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+ select case( Igeometry )
+
+ case( 1 ) ! Igeometry = 1 ; 09 Mar 17;
+
+ pxyz(1:3) = (/ teta , zeta, dR(0) /)
+ XXt = one ; XXz = zero
+ YYt = zero ; YYz = one
+ ZZt = dR(2) ; ZZz = dR(3)
+
+ case( 2 ) ! Igeometry = 2 ; 09 Mar 17;
+
+ FATAL( casing, .true., virtual casing under construction for cylindrical geometry )
+
+ case( 3 ) ! Igeometry = 3 ; toroidal geometry;
+
+ czeta = cos( zeta ) ; szeta = sin( zeta )
+
+ pxyz(1:3) = (/ dR(0) * czeta, dR(0) * szeta , dZ(0) /)
+ XXt = dR(2) * czeta ; XXz = dR(3) * czeta - dR(0) * szeta ! 10 Apr 13;
+ YYt = dR(2) * szeta ; YYz = dR(3) * szeta + dR(0) * czeta
+ ZZt = dZ(2) ; ZZz = dZ(3)
+
+ end select
+
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+ jj(1:3) = (/ Blz * XXt - Blt * XXz, &
+ Blz * YYt - Blt * YYz, &
+ Blz * ZZt - Blt * ZZz /)
+
+ return
+
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+ end subroutine surfacecurrent
+
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!> @brief Differential virtual casing integrand
!> \ingroup grp_free-boundary
!>
-!> Differential virtual casing integrand
+!> Differential virtual casing integrand with the calling convention required by NAG for the adaptive integration routine.
+!> This function wraps surfacecurrent and uses a global variable globaljk to determine the position on the computational boundary (not thread safe!).
!>
!> @param[in] Ndim number of parameters (==2)
!> @param[in] tz \f$\theta\f$ and \f$\zeta\f$
!> @param[in] Nfun number of function values (==3)
!> @param[out] vcintegrand cartesian components of magnetic field
-subroutine dvcfield( Ndim, tz, Nfun, vcintegrand ) ! differential virtual-casing field; format is fixed by NAG requirements;
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
- use constants, only : zero, half, one, three, four
+subroutine dvcfield( Ndim, tz, Nfun, vcintegrand )
- use numerical, only : small
+ use constants, only : zero, half, one, three
- use fileunits, only : ounit, vunit
-
- use inputlist, only : Wcasing, Nvol, Igeometry, Lrad, vcasingeps
-
- use cputiming, only :
+ use inputlist, only : vcasingeps
use allglobal, only : myid, ncpu, cpus, MPI_COMM_SPEC, &
- pi2nfp, &
- Mvol, &
- mn, im, in, &
- iRbc, iZbs, iRbs, iZbc, &
- Ate, Aze, Ato, Azo, &
- TT, &
- YESstellsym, NOTstellsym, &
- globaljk, Dxyz, Nxyz, &
- first_free_bound
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+ globaljk, Dxyz, Nxyz, first_free_bound
+
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
LOCALS
-
+
INTEGER , intent(in) :: Ndim, Nfun
REAL , intent(in) :: tz(1:Ndim)
REAL , intent(out) :: vcintegrand(1:Nfun) ! integrand; components of magnetic field due to plasma currents in Cartesian coordinates;
-
- INTEGER :: ii, mi, ni, ll, ideriv, jk
- REAL :: dR(0:3), dZ(0:3), gBut, gBuz, gtt, gtz, gzz, sqrtg, Blt, Blz, czeta, szeta, arg, carg, sarg, XX, YY, ZZ, teta, zeta
- REAL :: jj(1:3), rr(1:3), distance(1:3), firstorderfactor
-
- REAL :: XXt, XXz, YYt, YYz, ZZt, ZZz, ds, Bxyz(1:3)
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+ REAL :: jj(1:3), & ! Cartesian surface current components on the plasma boundary
+ pxyz(1:3), & ! Position on the plasma boundary
+ rr(1:3), & ! distance vector between evaluation point on the computational boundary and the current surface point
+ distance(1:3), & ! powers of distance |(x-x')|, |(x-x')|^2, |(x-x')|^3 including a small eps to avoid division by zero
+ Bxyz(1:3), firstorderfactor
#ifdef DEBUG
FATAL( casing, Ndim.ne. 2, incorrect )
FATAL( casing, Nfun.ne. 1, incorrect )
#endif
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
- dR(0:3) = zero ; ideriv = 0 ; gBut = zero ; gBuz = zero ! initialize summation of coordinates and tangential field;
- dZ(0:3) = zero
-
- teta = tz(1) ; zeta = tz(2) ! shorthand; 09 Mar 17;
-
- jk = globaljk
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
- select case( Igeometry )
-
- case( 1 ) ! Igeometry = 1 ; 09 Mar 17;
-
- if( YESstellsym ) then
-
- do ii = 1, mn ; mi = im(ii) ; ni = in(ii) ! loop over Fourier modes; construct surface current; slow transform required as position is arbitrary;
-
- arg = mi * teta - ni * zeta ; carg = cos(arg) ; sarg = sin(arg)
-
- dR(0) = dR(0) + ( iRbc(ii,Nvol) ) * carg
- dR(1) = dR(1) + ( ( iRbc(ii,Mvol) - iRbc(ii,Nvol) ) * carg ) * half
- dR(2) = dR(2) + mi * ( - ( iRbc(ii,Nvol) ) * sarg )
- dR(3) = dR(3) - ni * ( - ( iRbc(ii,Nvol) ) * sarg )
-
- !dZ(0) = dZ(0) + ( iZbs(ii,Nvol) ) * sarg
- !dZ(1) = dZ(1) + ( ( iZbs(ii,Mvol) - iZbs(ii,Nvol) ) * sarg ) * half
- !dZ(2) = dZ(2) + mi * ( ( iZbs(ii,Nvol) ) * carg )
- !dZ(3) = dZ(3) - ni * ( ( iZbs(ii,Nvol) ) * carg )
-
- do ll = 0, Lrad(Mvol)
- gBut = gBut - ( Aze(Mvol,ideriv,ii)%s(ll) * carg ) * TT(ll,0,1) ! contravariant; Jacobian comes later;
- gBuz = gBuz + ( Ate(Mvol,ideriv,ii)%s(ll) * carg ) * TT(ll,0,1)
- enddo
-
- enddo ! end of do ii = 1, mn ;
-
- else ! NOTstellsym ; 08 Feb 16;
-
- do ii = 1, mn ; mi = im(ii) ; ni = in(ii) ! loop over Fourier modes; construct surface current; slow transform required as position is arbitrary;
-
- arg = mi * teta - ni * zeta ; carg = cos(arg) ; sarg = sin(arg)
-
- dR(0) = dR(0) + iRbc(ii,Nvol) * carg + iRbs(ii,Nvol) * sarg
- dR(1) = dR(1) + ( ( iRbc(ii,Mvol) - iRbc(ii,Nvol) ) * carg + ( iRbs(ii,Mvol) - iRbs(ii,Nvol) ) * sarg ) * half
- dR(2) = dR(2) + mi * ( - iRbc(ii,Nvol) * sarg + iRbs(ii,Nvol) * carg )
- dR(3) = dR(3) - ni * ( - iRbc(ii,Nvol) * sarg + iRbs(ii,Nvol) * carg )
-
- !dZ(0) = dZ(0) + iZbc(ii,Nvol) * carg + iZbs(ii,Nvol) * sarg
- !dZ(1) = dZ(1) + ( ( iZbc(ii,Mvol) - iZbc(ii,Nvol) ) * carg + ( iZbs(ii,Mvol) - iZbs(ii,Nvol) ) * sarg ) * half
- !dZ(2) = dZ(2) + mi * ( - iZbc(ii,Nvol) * sarg + iZbs(ii,Nvol) * carg )
- !dZ(3) = dZ(3) - ni * ( - iZbc(ii,Nvol) * sarg + iZbs(ii,Nvol) * carg )
-
- do ll = 0, Lrad(Mvol)
- gBut = gBut - ( Aze(Mvol,ideriv,ii)%s(ll) * carg + Azo(Mvol,ideriv,ii)%s(ll) * sarg ) * TT(ll,0,1) ! contravariant; Jacobian comes later;
- gBuz = gBuz + ( Ate(Mvol,ideriv,ii)%s(ll) * carg + Ato(Mvol,ideriv,ii)%s(ll) * sarg ) * TT(ll,0,1)
- enddo
-
- enddo ! end of do ii = 1, mn ;
-
- endif ! end of if( YESstellsym ) ; 08 Feb 16;
-
- case( 2 ) ! Igeometry = 2 ; 09 Mar 17;
-
- FATAL( casing, .true., virtual casing under construction for cylindrical geometry )
-
- case( 3 ) ! Igeometry = 3 ; 09 Mar 17;
-
- if( YESstellsym ) then
-
- do ii = 1, mn ; mi = im(ii) ; ni = in(ii) ! loop over Fourier modes; construct surface current; slow transform required as position is arbitrary;
-
- arg = mi * teta - ni * zeta ; carg = cos(arg) ; sarg = sin(arg)
- dR(0) = dR(0) + ( iRbc(ii,Nvol) ) * carg
- dR(1) = dR(1) + ( ( iRbc(ii,Mvol) - iRbc(ii,Nvol) ) * carg ) * half
- dR(2) = dR(2) + mi * ( - ( iRbc(ii,Nvol) ) * sarg )
- dR(3) = dR(3) - ni * ( - ( iRbc(ii,Nvol) ) * sarg )
-
- dZ(0) = dZ(0) + ( iZbs(ii,Nvol) ) * sarg
- dZ(1) = dZ(1) + ( ( iZbs(ii,Mvol) - iZbs(ii,Nvol) ) * sarg ) * half
- dZ(2) = dZ(2) + mi * ( ( iZbs(ii,Nvol) ) * carg )
- dZ(3) = dZ(3) - ni * ( ( iZbs(ii,Nvol) ) * carg )
-
- if (first_free_bound) then
- do ll = 0, Lrad(Nvol) ! 1 is for outside thr volume
- gBut = gBut - ( Aze(Nvol,ideriv,ii)%s(ll) * carg ) * TT(ll,1,1) ! contravariant; Jacobian comes later;
- gBuz = gBuz + ( Ate(Nvol,ideriv,ii)%s(ll) * carg ) * TT(ll,1,1)
- enddo
- else
- do ll = 0, Lrad(Mvol)
- gBut = gBut - ( Aze(Mvol,ideriv,ii)%s(ll) * carg ) * TT(ll,0,1) ! contravariant; Jacobian comes later;
- gBuz = gBuz + ( Ate(Mvol,ideriv,ii)%s(ll) * carg ) * TT(ll,0,1)
- enddo
- endif
-
- enddo ! end of do ii = 1, mn ;
-
- else ! NOTstellsym ; 08 Feb 16;
-
- do ii = 1, mn ; mi = im(ii) ; ni = in(ii) ! loop over Fourier modes; construct surface current; slow transform required as position is arbitrary;
-
- arg = mi * teta - ni * zeta ; carg = cos(arg) ; sarg = sin(arg)
-
- dR(0) = dR(0) + iRbc(ii,Nvol) * carg + iRbs(ii,Nvol) * sarg
- dR(1) = dR(1) + ( ( iRbc(ii,Mvol) - iRbc(ii,Nvol) ) * carg + ( iRbs(ii,Mvol) - iRbs(ii,Nvol) ) * sarg ) * half
- dR(2) = dR(2) + mi * ( - iRbc(ii,Nvol) * sarg + iRbs(ii,Nvol) * carg )
- dR(3) = dR(3) - ni * ( - iRbc(ii,Nvol) * sarg + iRbs(ii,Nvol) * carg )
-
- dZ(0) = dZ(0) + iZbc(ii,Nvol) * carg + iZbs(ii,Nvol) * sarg
- dZ(1) = dZ(1) + ( ( iZbc(ii,Mvol) - iZbc(ii,Nvol) ) * carg + ( iZbs(ii,Mvol) - iZbs(ii,Nvol) ) * sarg ) * half
- dZ(2) = dZ(2) + mi * ( - iZbc(ii,Nvol) * sarg + iZbs(ii,Nvol) * carg )
- dZ(3) = dZ(3) - ni * ( - iZbc(ii,Nvol) * sarg + iZbs(ii,Nvol) * carg )
-
- if (first_free_bound) then
- do ll = 0, Lrad(Mvol) ! omit the possible current sheet due to a jump in tangential field at the plasma boundary; Zhu 20190603;
- gBut = gBut - ( Aze(Mvol,ideriv,ii)%s(ll) * carg + Azo(Mvol,ideriv,ii)%s(ll) * sarg ) * TT(ll,1,1) ! contravariant; Jacobian comes later;
- gBuz = gBuz + ( Ate(Mvol,ideriv,ii)%s(ll) * carg + Ato(Mvol,ideriv,ii)%s(ll) * sarg ) * TT(ll,1,1)
- enddo
- else
- do ll = 0, Lrad(Mvol)
- gBut = gBut - ( Aze(Mvol,ideriv,ii)%s(ll) * carg + Azo(Mvol,ideriv,ii)%s(ll) * sarg ) * TT(ll,0,1) ! contravariant; Jacobian comes later;
- gBuz = gBuz + ( Ate(Mvol,ideriv,ii)%s(ll) * carg + Ato(Mvol,ideriv,ii)%s(ll) * sarg ) * TT(ll,0,1)
- enddo
- endif
-
- enddo ! end of do ii = 1, mn ;
-
- endif ! end of if( YESstellsym ) ; 08 Feb 16;
-
- end select ! end of select case( Igeometry ) ; 09 Mar 17;
-
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
- select case( Igeometry )
-
- case( 1 ) ! Igeometry = 1 ; 09 Mar 17;
-
- gtt = one + dR(2)*dR(2)
- gtz = dR(2)*dR(3)
- gzz = one + dR(3)*dR(3)
-
- sqrtg = dR(1)
-
- case( 2 ) ! Igeometry = 2 ; 09 Mar 17;
+ call surfacecurrent(tz(1), tz(2), pxyz, jj)
- FATAL( casing, .true., virtual casing under construction for cylindrical geometry )
-
- case( 3 ) ! Igeometry = 3 ; 09 Mar 17;
-
- gtt = dR(2)*dR(2) + dZ(2)*dZ(2)
- gtz = dR(2)*dR(3) + dZ(2)*dZ(3)
- gzz = dR(3)*dR(3) + dZ(3)*dZ(3) + dR(0)*dR(0)
-
- sqrtg = dR(0) * ( dZ(1) * dR(2) - dR(1) * dZ(2) )
-
- end select
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
- Blt = ( gBut * gtt + gBuz * gtz ) / sqrtg
- Blz = ( gBut * gtz + gBuz * gzz ) / sqrtg
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
- select case( Igeometry )
-
- case( 1 ) ! Igeometry = 1 ; 09 Mar 17;
-
- XX = teta ; XXt = one ; XXz = zero
- YY = zeta ; YYt = zero ; YYz = one
- ZZ = dR(0) ; ZZt = dR(2) ; ZZz = dR(3)
-
- case( 2 ) ! Igeometry = 2 ; 09 Mar 17;
-
- FATAL( casing, .true., virtual casing under construction for cylindrical geometry )
-
- case( 3 ) ! Igeometry = 3 ; toroidal geometry;
-
- czeta = cos( zeta ) ; szeta = sin( zeta )
-
- XX = dR(0) * czeta ; XXt = dR(2) * czeta ; XXz = dR(3) * czeta - dR(0) * szeta ! 10 Apr 13;
- YY = dR(0) * szeta ; YYt = dR(2) * szeta ; YYz = dR(3) * szeta + dR(0) * czeta
- ZZ = dZ(0) ; ZZt = dZ(2) ; ZZz = dZ(3)
-
- end select
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
- rr(1:3) = (/ Dxyz(1,jk) - XX, &
- Dxyz(2,jk) - YY, &
- Dxyz(3,jk) - ZZ /)
-
- jj(1:3) = (/ Blz * XXt - Blt * XXz, &
- Blz * YYt - Blt * YYz, &
- Blz * ZZt - Blt * ZZz /)
+ ! distance vector between point on computational boundary and point on plasma boundary
+ rr(1:3) = Dxyz(1:3,globaljk) - pxyz(1:3)
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
distance(2) = sum( rr(1:3) * rr(1:3) ) + vcasingeps**2 ! 04 May 17;
-
- distance(1) = sqrt( distance(2) ) ; distance(3) = distance(1) * distance(2) ! powers of distance; 24 Nov 16;
+ distance(1) = sqrt( distance(2) ) ;
+ distance(3) = distance(1) * distance(2) ! powers of distance; 24 Nov 16;
firstorderfactor = ( one + three * vcasingeps**2 / distance(2) ) / distance(3) ! 04 May 17;
@@ -517,28 +615,23 @@ subroutine dvcfield( Ndim, tz, Nfun, vcintegrand ) ! differential virtual-casing
jj(3) * rr(1) - jj(1) * rr(3), &
jj(1) * rr(2) - jj(2) * rr(1) /)
- vcintegrand(1) = sum( Bxyz(1:3) * Nxyz(1:3,jk) ) * firstorderfactor
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+ vcintegrand(1) = sum( Bxyz(1:3) * Nxyz(1:3,globaljk) ) * firstorderfactor
-! vcintegrand( 4) = - three * vcintegrand(1) * rr(1) / distance(2) ! dBxdx; ! need to divide by distance(3) ; 14 Apr 17;
-! vcintegrand( 5) = - three * vcintegrand(1) * rr(2) / distance(2) - jj(3) ! dBxdy;
-! vcintegrand( 6) = - three * vcintegrand(1) * rr(3) / distance(2) + jj(2) ! dBxdz;
-!
-! vcintegrand( 7) = - three * vcintegrand(2) * rr(1) / distance(2) + jj(3) ! dBydx;
-! vcintegrand( 8) = - three * vcintegrand(2) * rr(2) / distance(2) ! dBydy;
-! vcintegrand( 9) = - three * vcintegrand(2) * rr(3) / distance(2) - jj(1) ! dBydz;
-!
-! vcintegrand(10) = - three * vcintegrand(3) * rr(1) / distance(2) - jj(2) ! dBzdx;
-! vcintegrand(11) = - three * vcintegrand(3) * rr(2) / distance(2) + jj(1) ! dBzdy;
-! vcintegrand(12) = - three * vcintegrand(3) * rr(3) / distance(2) ! dBzdz;
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+ ! vcintegrand( 4) = - three * vcintegrand(1) * rr(1) / distance(2) ! dBxdx; ! need to divide by distance(3) ; 14 Apr 17;
+ ! vcintegrand( 5) = - three * vcintegrand(1) * rr(2) / distance(2) - jj(3) ! dBxdy;
+ ! vcintegrand( 6) = - three * vcintegrand(1) * rr(3) / distance(2) + jj(2) ! dBxdz;
+ !
+ ! vcintegrand( 7) = - three * vcintegrand(2) * rr(1) / distance(2) + jj(3) ! dBydx;
+ ! vcintegrand( 8) = - three * vcintegrand(2) * rr(2) / distance(2) ! dBydy;
+ ! vcintegrand( 9) = - three * vcintegrand(2) * rr(3) / distance(2) - jj(1) ! dBydz;
+ !
+ ! vcintegrand(10) = - three * vcintegrand(3) * rr(1) / distance(2) - jj(2) ! dBzdx;
+ ! vcintegrand(11) = - three * vcintegrand(3) * rr(2) / distance(2) + jj(1) ! dBzdy;
+ ! vcintegrand(12) = - three * vcintegrand(3) * rr(3) / distance(2) ! dBzdz;
+ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
return
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
end subroutine dvcfield
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
diff --git a/src/global.f90 b/src/global.f90
index 902cd112..6c5b1ea9 100644
--- a/src/global.f90
+++ b/src/global.f90
@@ -836,19 +836,18 @@ module allglobal
!> \addtogroup grp_misc Miscellaneous
!> The following are miscellaneous flags required for the virtual casing field, external (vacuum) field integration, ...
!> @{
- INTEGER :: globaljk !< labels position
- REAL, allocatable :: Dxyz(:,:) !< computational boundary; position
- REAL, allocatable :: Nxyz(:,:) !< computational boundary; normal
- REAL, allocatable :: Jxyz(:,:) !< plasma boundary; surface current
-
- REAL :: tetazeta(1:2) !< what is this?
+ INTEGER :: globaljk !< labels position
+ REAL, allocatable :: Dxyz(:,:) !< computational boundary; position
+ REAL, allocatable :: Nxyz(:,:) !< computational boundary; normal
+ REAL, allocatable :: Jxyz(:,:) !< plasma boundary; surface current
+ REAL, allocatable :: Pbxyz(:,:) !< plasma boundary; position
+ INTEGER :: prevcstride!< previous virtual casing stride for casinggrid (only relevant when Lvcgrid=1)
REAL :: virtualcasingfactor = -one / (four*pi) !< this agrees with diagno
-
+
INTEGER :: IBerror !< for computing error in magnetic field
INTEGER :: nfreeboundaryiterations !< number of free-boundary iterations already performed
-
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
INTEGER, parameter :: Node = 2 !< best to make this global for consistency between calling and called routines
@@ -1345,7 +1344,7 @@ subroutine check_inputs()
1020 format("readin : ",f10.2," : Linitialize=",i3," ;LautoinitBn=",i3," ; Lzerovac=",i2," ; Ndiscrete="i2" ;")
1021 format("readin : ", 10x ," : Nquad="i4" ; iMpol="i4" ; iNtor="i4" ;")
1022 format("readin : ", 10x ," : Lsparse="i2" ; Lsvdiota="i2" ; imethod="i2" ; iorder="i2" ; iprecon="i2" ; iotatol="es13.5" ;")
-1023 format("readin : ", 10x ," : Lextrap="i2" ; Mregular="i3" ; Lrzaxis="i2" ; Ntoraxis="i2" ;")
+1023 format("readin : ", 10x ," : Lextrap="i2" ; Mregular="i3" ; Lrzaxis="i2" ; Ntoraxis="i2" ; Lvcgrid="i2" ;")
FATAL( readin, Ndiscrete.le.0, error )
@@ -1386,13 +1385,13 @@ subroutine check_inputs()
write(ounit,1041) escale, opsilon, pcondense, epsilon, wpoloidal, upsilon
write(ounit,1042) forcetol, c05xmax, c05xtol, c05factor, LreadGF
write(ounit,1043) mfreeits, gBntol, gBnbld
- write(ounit,1044) vcasingeps, vcasingtol, vcasingits, vcasingper
+ write(ounit,1044) vcasingeps, vcasingtol, vcasingits, vcasingper, vcNt, vcNz
1040 format("readin : ",f10.2," : Lfindzero="i2" ;")
1041 format("readin : ", 10x ," : escale="es13.5" ; opsilon="es13.5" ; pcondense="f7.3" ; epsilon="es13.5" ; wpoloidal="f7.4" ; upsilon="es13.5" ;")
1042 format("readin : ", 10x ," : forcetol="es13.5" ; c05xmax="es13.5" ; c05xtol="es13.5" ; c05factor="es13.5" ; LreadGF="L2" ; ")
1043 format("readin : ", 10x ," : mfreeits="i4" ; gBntol="es13.5" ; gBnbld="es13.5" ;")
-1044 format("readin : ", 10x ," : vcasingeps="es13.5" ; vcasingtol="es13.5" ; vcasingits="i6" ; vcasingper="i6" ;")
+1044 format("readin : ", 10x ," : vcasingeps="es13.5" ; vcasingtol="es13.5" ; vcasingits="i6" ; vcasingper="i6" ; vcNt="i6" ; vcNz="i6" ;")
FATAL( readin, escale .lt.zero , error )
FATAL( readin, pcondense .lt.one , error )
@@ -1513,6 +1512,7 @@ subroutine broadcast_inputs
IlBCAST( Mregular , 1, 0 )
IlBCAST( Lrzaxis , 1, 0 )
IlBCAST( Ntoraxis , 1, 0 )
+ IlBCAST( Lvcgrid , 1, 0 )
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
@@ -1540,6 +1540,8 @@ subroutine broadcast_inputs
RlBCAST( vcasingtol, 1 , 0 )
IlBCAST( vcasingits, 1 , 0 )
IlBCAST( vcasingper, 1 , 0 )
+ IlBCAST( vcNt, 1 , 0 )
+ IlBCAST( vcNz, 1 , 0 )
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
@@ -1824,6 +1826,7 @@ subroutine wrtend
write(iunit,'(" Mregular = ",i9 )') Mregular
write(iunit,'(" Lrzaxis = ",i9 )') Lrzaxis
write(iunit,'(" Ntoraxis = ",i9 )') Ntoraxis
+ write(iunit,'(" Lvcgrid = ",i9 )') Lvcgrid
write(iunit,'("/")')
if( Wwrtend ) then ; cput = GETTIME ; write(ounit,'("wrtend : ",f10.2," : myid=",i3," ; writing locallist ;")') cput-cpus, myid
@@ -1865,6 +1868,8 @@ subroutine wrtend
write(iunit,'(" vcasingtol = ",es23.15 )') vcasingtol
write(iunit,'(" vcasingits = ",i9 )') vcasingits
write(iunit,'(" vcasingper = ",i9 )') vcasingper
+ write(iunit,'(" vcNt = ",i9 )') vcNt
+ write(iunit,'(" vcNz = ",i9 )') vcNz
write(iunit,'("/")')
if( Wwrtend ) then ; cput = GETTIME ; write(ounit,'("wrtend : ",f10.2," : myid=",i3," ; writing diagnosticslist ;")') cput-cpus, myid
diff --git a/src/inputlist.f90 b/src/inputlist.f90
index 30c0b04b..24d960d4 100644
--- a/src/inputlist.f90
+++ b/src/inputlist.f90
@@ -315,6 +315,7 @@ module inputlist
!<
INTEGER :: Ntoraxis = 3 !< the number of \f$n\f$ harmonics used in the Jacobian \f$m=1\f$ harmonic elimination method;
!< only relevant if \c Lrzaxis.ge.1 .
+ INTEGER :: Lvcgrid = 0 !< Which method to use for the virtual casing integral. 0 = adaptive integration routine with guaranteed accuracy, 1 = fixed resolution grid
!> @}
!> \addtogroup grp_global_local locallist
@@ -469,8 +470,10 @@ module inputlist
REAL :: vcasingtol = 1.e-08 !< accuracy on virtual casing integral; see bnorml(), casing()
INTEGER :: vcasingits = 8 !< minimum number of calls to adaptive virtual casing routine; see casing()
INTEGER :: vcasingper = 1 !< periods of integragion in adaptive virtual casing routine; see casing()
+ INTEGER :: vcNt = 256 !< theta resolution of the real space grid on which the surface current is computed during virtual casing [0, 2pi]
+ INTEGER :: vcNz = 256 !< zeta resolution of the real space grid on which the surface current is computed during virtual casing [0, 2pi], over one field period
INTEGER :: mcasingcal = 8 !< minimum number of calls to adaptive virtual casing routine; see casing(); redundant;
-!> @}
+ !> @}
!> \addtogroup grp_global_diagnostics diagnosticslist
@@ -675,7 +678,8 @@ module inputlist
Lextrap ,&
Mregular ,&
Lrzaxis ,&
- Ntoraxis
+ Ntoraxis ,&
+ Lvcgrid
namelist/locallist/&
LBeltrami ,&
@@ -711,6 +715,8 @@ module inputlist
vcasingtol ,&
vcasingits ,&
vcasingper ,&
+ vcNt ,&
+ vcNz ,&
mcasingcal
namelist/diagnosticslist/&
@@ -889,6 +895,7 @@ subroutine initialize_inputs
Mregular = -1
Lrzaxis = 1
Ntoraxis = 3
+ Lvcgrid = 0
! locallist
@@ -925,6 +932,8 @@ subroutine initialize_inputs
vcasingtol = 1.e-08
vcasingits = 8
vcasingper = 1
+ vcNt = 256
+ vcNz = 256
mcasingcal = 8
! diagnosticslist
diff --git a/src/preset.f90 b/src/preset.f90
index e3a6102a..4ecd2e42 100644
--- a/src/preset.f90
+++ b/src/preset.f90
@@ -1735,21 +1735,31 @@ subroutine preset
SALLOCATE( dlambdaout, (1:lmns,1:Mvol,0:1), zero )
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
if (Lfreebound > 0) then ! Only do for free-boundary; 7 Nov 18;
+ if (Lvcgrid .eq. 1) then
+ if ( vcNt.lt.Nt ) then
+ FATAL( bnorml, .true., The plasma boundary resolution for virtual casing vcNt must be greater than or equal to Nt )
+ endif
+ if ( vcNz.lt.Nz ) then
+ FATAL( bnorml, .true., The plasma boundary resolution for virtual casing vcNz must be greater than or equal to Nz )
+ endif
+ ! Factor of nfp because plasma boundary for virtual casing has to cover all field periods, but resolution is given for one field period.
+ SALLOCATE( Jxyz, (1:vcNt*vcNz*Nfp,1:3), zero ) ! Cartesian components of virtual casing surface current; needs to be recalculated at each iteration;
+ SALLOCATE( Pbxyz, (1:vcNt*vcNz*Nfp,1:3), zero ) ! Cartesian points on the plasma boundary; needs to be recalculated at each iteration;
+ prevcstride = huge(prevcstride)
+ endif
SALLOCATE( Dxyz, (1:3,1:Ntz), zero ) ! Cartesian components of computational boundary; position; 14 Apr 17;
SALLOCATE( Nxyz, (1:3,1:Ntz), zero ) ! Cartesian components of computational boundary; normal ; 14 Apr 17;
- SALLOCATE( Jxyz, (1:Ntz,1:3), zero ) ! Cartesian components of virtual casing surface current; needs to be recalculated at each iteration;
-
lvol = Mvol ; lss = one ; Lcurvature = 1 ; Lcoordinatesingularity = .false. ! will only require normal field on outer interface = computational boundary;
WCALL( preset, coords,( lvol, lss, Lcurvature, Ntz, mn ) ) ! will need Rij, Zij; THE COMPUTATIONAL BOUNDARY DOES NOT CHANGE;
do kk = 0, Nz-1 ; zeta = kk * pi2nfp / Nz
- if( Igeometry.eq.3 ) then ; cszeta(0:1) = (/ cos(zeta), sin(zeta) /)
+ if( Igeometry.eq.3 ) then ;
+ cszeta(0:1) = (/ cos(zeta), sin(zeta) /)
endif
do jj = 0, Nt-1 ; teta = jj * pi2 / Nt ; jk = 1 + jj + kk*Nt
diff --git a/src/sphdf5.f90 b/src/sphdf5.f90
index f6c833bf..ca8e0303 100644
--- a/src/sphdf5.f90
+++ b/src/sphdf5.f90
@@ -295,6 +295,7 @@ subroutine mirror_input_to_outfile
HWRITEIV( grpInputNumerics, 1, Mregular , (/ Mregular /))
HWRITEIV( grpInputNumerics, 1, Lrzaxis , (/ Lrzaxis /))
HWRITEIV( grpInputNumerics, 1, Ntoraxis , (/ Ntoraxis /))
+ HWRITEIV( grpInputNumerics, 1, Lvcgrid , (/ Lvcgrid /))
HCLOSEGRP( grpInputNumerics, __FILE__, __LINE__)
@@ -345,6 +346,8 @@ subroutine mirror_input_to_outfile
HWRITEIV( grpInputGlobal, 1, vcasingits , (/ vcasingits /))
HWRITEIV( grpInputGlobal, 1, vcasingper , (/ vcasingper /))
HWRITEIV( grpInputGlobal, 1, mcasingcal , (/ mcasingcal /)) ! redundant;
+ HWRITEIV( grpInputGlobal, 1, vcNt , (/ vcNt /))
+ HWRITEIV( grpInputGlobal, 1, vcNz , (/ vcNz /))
HCLOSEGRP( grpInputGlobal )