Merge branch 'master' of github.com:jonathanschilling/educational_VMEC

.gitmodules

@@ -7,3 +7,6 @@
 [submodule "src/vac3"]
 	path = src/vac3
 	url = [email protected]:jonathanschilling/vac3.git
+[submodule "abscab-fortran"]
+	path = abscab-fortran
+	url = https://github.com/jonathanschilling/abscab-fortran

CMakeLists.txt

@@ -51,6 +51,7 @@ endif ()
 set (vmec_sources "")
 add_subdirectory(src)
 add_subdirectory(json-fortran)
+add_subdirectory(abscab-fortran)
 
 # if you have access to the vac2 and vac3 versions of NESTOR (contact me), they are included here
 #add_subdirectory(src/vac2)

src/NESTOR/analyt.f90

@@ -115,17 +115,20 @@ SUBROUTINE analyt(grpmn, bvec, ivacskip, lasym, m_map, n_map, grpmn_m_map, grpmn
   LLOOP: DO l = 0, mf + nf
      fl = fl1
 
+     ! here, tlp/m are available for the current value of l
+     ! --> save into matrix for debugging
+     all_tlp(l,:) = tlp
+     all_tlm(l,:) = tlm
+
      ! COMPUTE SL+ and SL- , Eq (A17)
      ! SLP(M): SL+(-)
      IF (ivacskip .eq. 0) THEN
         slp = (r1p*fl + ra1p)*tlp + r0p*fl*tlp1 - (r1p + r0p)/sqrtc + sign1*(r0p - r1p)/sqrta
         slm = (r1m*fl + ra1m)*tlm + r0m*fl*tlm1 - (r1m + r0m)/sqrtc + sign1*(r0m - r1m)/sqrta
         slpm = slp + slm
 
-        ! here, tlp/m and slp/m are available for the current value of l
+        ! here, slp/m are available for the current value of l
         ! --> save into matrix for debugging
-        all_tlp(l,:) = tlp
-        all_tlm(l,:) = tlm
         all_slp(l,:) = slp
         all_slm(l,:) = slm
      ENDIF
@@ -169,11 +172,12 @@ SUBROUTINE analyt(grpmn, bvec, ivacskip, lasym, m_map, n_map, grpmn_m_map, grpmn
 
     call add_real_3d("all_tlp", mf+nf+1, nv, nu3, all_tlp)
     call add_real_3d("all_tlm", mf+nf+1, nv, nu3, all_tlm)
-    call add_real_2d("bvec", mf1, nf1, bvec)
+    call add_real_2d("bvec", mf1, nf1, bvec) ! (mf+1)x(2*nf+1)x(ndim: 1 or 2)
+
     if (ivacskip .eq. 0) then
       call add_real_3d("all_slp", mf+nf+1, nv, nu3, all_slp)
       call add_real_3d("all_slm", mf+nf+1, nv, nu3, all_slm)
-      call add_real_4d("grpmn", mf1, nf1, nv, nu3, grpmn)
+      call add_real_4d("grpmn", mf1, nf1, nv, nu3, grpmn) ! missing dim: (ndim: 1 or 2)
     else
       call add_null("all_slp")
       call add_null("all_slm")

src/NESTOR/belicu.f90

@@ -15,6 +15,8 @@ SUBROUTINE belicu(torcur, bx, by, bz, cos1, sin1, rp, zp)
 
   USE vacmod, vm_bz => bz
 
+  use abscab, only: magneticFieldPolygonFilament
+
   IMPLICIT NONE
 
   REAL(rprec), INTENT(in) :: torcur
@@ -23,63 +25,48 @@ SUBROUTINE belicu(torcur, bx, by, bz, cos1, sin1, rp, zp)
 
   REAL(rprec) :: current
   INTEGER :: i, j, kper, kv
-  REAL(rprec), DIMENSION(3) :: xpt, bvec, dvec, Ri_vec
 
-  ! quantities from Hanson & Hirshman, "Compact expressions for the Biot-Savart fields of a filamentary segment" (2002)
-  REAL(rprec) :: L, Ri, Rf, Ri_p_Rf, Bmag
+  real(rprec), dimension(3, nuv2) :: eval_pos, magnetic_field
 
   ! B_External due to LIne CUrrent
 
   ! net toroidal plasma current in A
   current = torcur/mu0
-
-  ! initialize target array
-  bx = 0
-  by = 0
-  bz = 0
-
-  ! first point (at index 0) is equal to last point --> closed curve
-  xpts(1, 0) = raxis_nestor(nv)*(cosper(nvper)*cosuv(nv) - sinper(nvper)*sinuv(nv))
-  xpts(2, 0) = raxis_nestor(nv)*(sinper(nvper)*cosuv(nv) + cosper(nvper)*sinuv(nv))
-  xpts(3, 0) = zaxis_nestor(nv)
 
   ! loops over source geometry
-  i = 1
+  i = 0
   DO kper = 1, nvper
      DO kv = 1, nv
+        i = i + 1
 
         ! xpts == xpt of _s_ource (current filament)
         xpts(1, i) = raxis_nestor(kv)*(cosper(kper)*cosuv(kv) - sinper(kper)*sinuv(kv))
         xpts(2, i) = raxis_nestor(kv)*(sinper(kper)*cosuv(kv) + cosper(kper)*sinuv(kv))
         xpts(3, i) = zaxis_nestor(kv)
+     end do
+  end do
 
-        ! filament geometry: from current point (R_i == xpts(:,i)) to previous point (R_f == xpts(:,i-1))
-        dvec = xpts(:,i)-xpts(:,i-1)
-        L = norm2(dvec)
+  ! last point is equal to first point --> closed curve
+  xpts(1, nvp+1) = xpts(1, 1)
+  xpts(2, nvp+1) = xpts(2, 1)
+  xpts(3, nvp+1) = xpts(3, 1)
 
-        ! loop over evaluation points
-        DO j = 1,nuv2
-           ! xpt == evaluation position
-           xpt(1) = rp(j) * cos1(j)
-           xpt(2) = rp(j) * sin1(j)
-           xpt(3) = zp(j)
+  DO j = 1, nuv2
+    ! evaluation positions
+    eval_pos(1, j) = rp(j) * cos1(j)
+    eval_pos(2, j) = rp(j) * sin1(j)
+    eval_pos(3, j) = zp(j)
+  end do
 
-           Ri_vec = xpt - xpts(:,i-1)
-           Ri = norm2(Ri_vec)
-           Rf = norm2(xpt - xpts(:,i))
-           Ri_p_Rf = Ri + Rf
-
-           ! 1.0e-7 == mu0/4 pi
-           Bmag = 1.0E-7_dp * current * 2.0_dp * Ri_p_Rf / ( Ri * Rf * (Ri_p_Rf*Ri_p_Rf - L*L) )
+  ! initialize target array
+  magnetic_field = 0.0_dp
 
-           ! cross product of L*hat(eps)==dvec with Ri_vec, scaled by Bmag
-           bx(j) = bx(j) + Bmag * (dvec(2)*Ri_vec(3) - dvec(3)*Ri_vec(2))
-           by(j) = by(j) + Bmag * (dvec(3)*Ri_vec(1) - dvec(1)*Ri_vec(3))
-           bz(j) = bz(j) + Bmag * (dvec(1)*Ri_vec(2) - dvec(2)*Ri_vec(1))
-        end do
+  ! use ABSCAB to compute the line-current-along-axis magnetic field contribution
+  call magneticFieldPolygonFilament(nvper * nv + 1, xpts, current, &
+                                    nuv2, eval_pos, magnetic_field)
 
-        i = i + 1
-     END DO
-  END DO
+  bx(:) = magnetic_field(1,:)
+  by(:) = magnetic_field(2,:)
+  bz(:) = magnetic_field(3,:)
 
 END SUBROUTINE belicu

src/NESTOR/bextern.f90

@@ -18,6 +18,7 @@ SUBROUTINE bextern(plascur, wint)
   REAL(rprec), DIMENSION(nuv2), INTENT(in) :: wint
 
   INTEGER :: i
+  logical :: dbgout_active
 
   ! exterior Neumann problem
 
@@ -30,6 +31,19 @@ SUBROUTINE bextern(plascur, wint)
   ! This sets brad, bphi and bz to the interpolated field from the mgrid.
   CALL becoil(r1b, z1b, bvac(1,1), bvac(1,2), bvac(1,3))
 
+  dbgout_active = open_dbg_context("vac1n_bextern", num_eqsolve_retries)
+  if (dbgout_active) then
+
+    ! these are only from the mgrid at this point
+    call add_real_2d("mgrid_brad", nv, nu3, brad)
+    call add_real_2d("mgrid_bphi", nv, nu3, bphi)
+    call add_real_2d("mgrid_bz",   nv, nu3, bz)
+
+    ! This should be in Amperes.
+    call add_real("axis_current", plascur/mu0)
+
+  end if ! dbgout_active
+
   ! COMPUTE B (ON PLASMA BOUNDARY) FROM NET TOROIDAL PLASMA CURRENT
   ! THE NET CURRENT IS MODELLED AS A WIRE AT THE MAGNETIC AXIS, AND THE
   ! BIOT-SAVART LAW IS USED TO COMPUTE THE FIELD AT THE PLASMA SURFACE
@@ -58,8 +72,12 @@ SUBROUTINE bextern(plascur, wint)
   ! NOTE: BEXN == NP*F = -B0 dot [Xu cross Xv] NP        (see PKM, Eq. 2.13)
   bexni(:nuv2) = wint(:nuv2)*bexn(:nuv2)*pi2*pi2
 
-  if (open_dbg_context("vac1n_bextern", num_eqsolve_retries)) then
+  if (dbgout_active) then
+
+    ! axis geometry used in belicu
+    call add_real_2d("xpts_axis", 3, nvper * nv + 1, xpts)
 
+    ! these are now mgrid + axis-current
     call add_real_2d("brad", nv, nu3, brad)
     call add_real_2d("bphi", nv, nu3, bphi)
     call add_real_2d("bz",   nv, nu3, bz)

src/NESTOR/data/vacmod.f90

@@ -229,8 +229,9 @@ subroutine allocate_nestor
   IF (i .ne. 0) STOP 'allocation error in bextern'
 
   ! from tolicu
-  ! need 0:nvp for "virtual" point at index 0 which is equal to last point for closed curve
-  ALLOCATE (xpts(3,0:nvp), stat=i)
+  ! need nvp+1 for "virtual" point at last index,
+  ! which is equal to first point for a closed curve
+  ALLOCATE (xpts(3,nvp+1), stat=i)
   IF (i .ne. 0) STOP ' allocation error in tolicu'
 
   ! from scalpot

src/bcovar.f90

@@ -394,7 +394,7 @@ SUBROUTINE bcovar (lu, lv)
     IF (MOD(iter2-iter1,ns4).eq.0) THEN
        ! only update preconditioner every ns4==25 iterations (?) (for ns4, see vmec_params)
 
-!        write(*,*) "update 1d preconditioner"
+       ! write(*,*) "update 1d preconditioner at iter2=", iter2
 
        phipog(:nrzt) = phipog(:nrzt)*wint(:nrzt) ! remember that actually phipog == 1/sqrt(g)