Merge pull request #12 from GKV-developers/maeyama_shearflow_rotating

Rotating flux tube model

README_for_namelist.txt

@@ -1,12 +1,13 @@
-NOTE for gkvp_f0.30                                     S. Maeyama   March 2013
-Updated for gkvp_f0.40                                   M. Nakata    June 2014
-Updated for gkvp_f0.45                                   M. Nakata    July 2015
-Updated for gkvp_f0.46                                  S. Maeyama     May 2016
-Updated for gkvp_f0.47                                  S. Maeyama     Nov 2016
-Updated for gkvp_f0.48                                  S. Maeyama     Dec 2016
-Updated for gkvp_f0.50                                  S. Maeyama     Sep 2017
-Updated for gkvp_f0.55                                   M. Nakata     Dec 2019
+Updated for gkvp_f0.62                                  S. Maeyama   March 2023
 Updated for gkvp_f0.61                                  S. Maeyama   March 2021
+Updated for gkvp_f0.55                                   M. Nakata     Dec 2019
+Updated for gkvp_f0.50                                  S. Maeyama     Sep 2017
+Updated for gkvp_f0.48                                  S. Maeyama     Dec 2016
+Updated for gkvp_f0.47                                  S. Maeyama     Nov 2016
+Updated for gkvp_f0.46                                  S. Maeyama     May 2016
+Updated for gkvp_f0.45                                   M. Nakata    July 2015
+Updated for gkvp_f0.40                                   M. Nakata    June 2014
+NOTE for gkvp_f0.30                                     S. Maeyama   March 2013
 
 %%% How to run the code %%%
 
@@ -66,6 +67,7 @@ equib_type: "analytic" -  Analytic helical field with the metrics in cylinder
             "vmec"     -  Tokamak/stellarator field from the VMEC code
             "eqdsk"    -  Tokamak field (MEUDAS/TOPICS or G-EQDSK) via IGS code
             "slab"     -  Shearless slab geometry
+            "ring"     -  Ring dipole geometry
 
 inum: current shot number
 
@@ -156,6 +158,15 @@ del_c: mode connection phase in fluxtube model
 
 eps_r ~~ malpha : geometrical parameters such as safety factor, B-shear, etc.  
 
+&ring: parameters for ring dipole geometry
+       !  There is a ring current at R=a. The field line passing through (R,Z)=(R0,0) is picked up as a flux-tube domain.
+       !  The reference length is set to be R0 (not the ring current at R=a).
+       !  The reference magnetic field strength is B0 at (R,Z)=(R0,0).
+
+ring_a: = a / R0, which specify a flux tube of the ring dipole.
+
+kxmin: Minimum wavenumber in kx, valid only when equib_type == "ring"
+
 &vmecp -- &bozxf : parameters for vmec equilibrium
 
 &igsp -- &igsf : parameters for tokamak (g-eqdsk) equilibrium

Version_memo.txt

@@ -1,3 +1,12 @@
+gkvp_f0.62                                              S. Maeyama     Mar 2023
+1) equib_type = "ring" is added for ring dipole geometry.
+
+2) Rotating flux-tube model is implemented to treat equilibrium shearflows,
+   available for torus: equib_type = "s-alpha", "s-alpha-shift", "analytic",
+   "circMHD", "vmec", "igs". (But not available for "slab", "ring")
+
+
+
 gkvp_f0.61                                              S. Maeyama     Mar 2021
 1) equib_type = "s-alpha-shift" is added. s-alpha model with Shafranov shift.
 
@@ -7,7 +16,6 @@ gkvp_f0.61                                              S. Maeyama     Mar 2021
 
 
 
-
 gkvp_f0.60                                              S. Maeyama     Feb 2021
 1) NetCDF4+parallel HDF5 is added for optional output of GKV.
 

lib/gkvp_math_portable.f90

@@ -5,6 +5,8 @@ MODULE GKV_math
 !
 !    Update history of gkvp_set.f90
 !    --------------
+!      gkvp_f0.62 (S. Maeyama, Mar 2023)
+!        - Elliptic integrals math_eli1, math_eli2 are added.
 !      gkvp_f0.61 (S. Maeyama, Mar 2021)
 !        - random_seed is added for reproducibility.
 !
@@ -17,7 +19,9 @@ MODULE GKV_math
   private
 
   public   math_j0, math_j0zero, math_j1, math_j2, &
-           math_i0, math_g0,     math_random
+           math_i0, math_g0,     math_random, &
+           math_eli1, math_eli2
+
 
   integer, parameter :: ifnc = 99 ! unit number preserved 
                                   ! for this module
@@ -158,6 +162,47 @@ SUBROUTINE math_g0( x, g0 )
 
   END SUBROUTINE math_g0
 
+  SUBROUTINE math_eli1( x, eli1 )
+!
+!     Complete elliptic integral of the first kind K
+!
+    real(kind=DP), intent(in)  :: x
+    real(kind=DP), intent(out) :: eli1
+!
+    real(kind=DP) :: sqrtx
+!
+      sqrtx = sqrt(x)
+      if( 0._DP <= sqrtx  .and.  sqrtx < 1._DP ) then
+        eli1 = ellipk(sqrtx)
+      else
+        print *, "### math_eli1:  x is out of range!"
+      end if
+
+    return
+
+  END SUBROUTINE math_eli1
+
+  SUBROUTINE math_eli2( x, eli2 )
+!
+!     0th-order modified Bessel function
+!
+    real(kind=DP), intent(in)  :: x
+    real(kind=DP), intent(out) :: eli2
+!
+    real(kind=DP) :: sqrtx
+!
+      sqrtx = sqrt(x)
+      if( 0._DP <= sqrtx  .and.  sqrtx < 1._DP ) then
+        eli2 = ellipe(sqrtx)
+      else
+        print *, "### math_eli2:  x is out of range!"
+      end if
+
+    return
+
+  END SUBROUTINE math_eli2
+
+
 
   SUBROUTINE math_random( rr )
 !
@@ -825,6 +870,160 @@ function dbesi1(x)
       dbesi1 = y
       end function dbesi1
 !
+FUNCTION ellipk(k)
+!
+! Complete elliptic integral of the 1st kind
+!   K(k) = \int_0^(pi/2) (1-k**2*sin(q)**2)**(-1/2) dq
+! where 0<=k<1.
+!
+! Algorithm based on Jacobi's nome q expansion [Refs1-3]
+! Implemented by S. Maeyama (July,2021)
+!
+! Ref1. T.Fukushima, "Fast computation of complete elliptic integrals and 
+!       Jacobian elliptic functions", Celest. Mech. Dyn. Astr. 105,
+!       305-328 (2009). DOI 10.1007/s10569-009-9228-z
+! Ref2. J.Yamauchi, T.Uno, S.Ichimatsu, "Denshikeisanki notameno
+!       suuchikeisanhou III" (Baifukan, 1971), p.258. (Japanese)
+! Ref3. http://www.totoha.net/fc2_mirror2/Complete_Elliptical_Integral.pdf
+!
+  real(kind=8) :: ellipk
+  real(kind=8), intent(in) :: k
+
+  real(kind=8), parameter :: halfpi=2.d0*atan(1.d0)
+  real(kind=8), parameter :: piinv=1.d0/(2.d0*halfpi)
+  real(kind=8) :: kp2,kp,sqrtkp,eps,eps4,q,q4,q9,q16,q25,t3
+
+  if (k<0.7d0) then
+
+    kp2=1.d0-k*k
+    kp=sqrt(kp2)
+    sqrtkp=sqrt(kp)
+    eps=0.5d0*(1.d0-sqrtkp)/(1.d0+sqrtkp)
+    eps4=eps*eps*eps*eps
+    q=eps*(1.d0+eps4*(2.d0+eps4*(15.d0  &
+     +eps4*(150.d0+eps4*(1707.d0+eps4*(20910.d0+eps4*(268616.d0)))))))
+    q4=q*q*q*q
+    q9=q4*q4*q
+    q16=q4*q4*q4*q4
+    q25=q16*q9
+    t3=1.d0+2.d0*(q+q4+q9+q16+q25)
+    ellipk=halfpi*t3*t3
+
+  else
+
+    kp2=k*k ! modify for large k
+    kp=sqrt(kp2)
+    sqrtkp=sqrt(kp)
+    eps=0.5d0*(1.d0-sqrtkp)/(1.d0+sqrtkp)
+    eps4=eps*eps*eps*eps
+    q=eps*(1.d0+eps4*(2.d0+eps4*(15.d0  &
+     +eps4*(150.d0+eps4*(1707.d0+eps4*(20910.d0+eps4*(268616.d0)))))))
+    q4=q*q*q*q
+    q9=q4*q4*q
+    q16=q4*q4*q4*q4
+    q25=q16*q9
+    t3=1.d0+2.d0*(q+q4+q9+q16+q25)
+    ellipk=halfpi*t3*t3
+    ellipk=-piinv*ellipk*log(q) ! modify for large k
+
+  end if
+
+    return
 
+END FUNCTION ellipk
+
+
+FUNCTION ellipe(k)
+!                                                                                                                                                                 
+! Complete elliptic integral of the 2nd kind
+!   E(k) = \int_0^(pi/2) (1-k**2*sin(q)**")**(1/2) dq
+! where 0<=k<1.
+!
+! Algorithm based on Jacobi's nome q expansion [Refs1-3]
+! Implemented by S. Maeyama (July,2021)
+!
+! Ref1. T.Fukushima, "Fast computation of complete elliptic integrals and 
+!       Jacobian elliptic functions", Celest. Mech. Dyn. Astr. 105,
+!       305-328 (2009). DOI 10.1007/s10569-009-9228-z
+! Ref2. J.Yamauchi, T.Uno, S.Ichimatsu, "Denshikeisanki notameno
+!       suuchikeisanhou III" (Baifukan, 1971), p.258. (Japanese)
+! Ref3. http://www.totoha.net/fc2_mirror2/Complete_Elliptical_Integral.pdf
+!
+  real(kind=8) :: ellipe
+  real(kind=8), intent(in) :: k
+
+  real(kind=8), parameter :: piover6=2.d0*atan(1.d0)/3.d0
+  real(kind=8), parameter :: halfpi=2.d0*atan(1.d0)
+  real(kind=8), parameter :: piinv=1.d0/(2.d0*halfpi)
+  real(kind=8) :: kp2,kp,sqrtkp,eps,eps4,q,q2,q4,q6,q8,q9,&
+                  q10,q12,q14,q16,q18,q20,q22,q25,t3,t3sq,fac,wq
+  real(kind=8) :: ellipkp, ellipk, ellipep
+
+  if (k<0.7d0) then
+
+    kp2=1.d0-k*k
+    kp=sqrt(kp2)
+    sqrtkp=sqrt(kp)
+    eps=0.5d0*(1.d0-sqrtkp)/(1.d0+sqrtkp)
+    eps4=eps*eps*eps*eps
+    q=eps*(1.d0+eps4*(2.d0+eps4*(15.d0  &
+     +eps4*(150.d0+eps4*(1707.d0+eps4*(20910.d0+eps4*(268616.d0)))))))
+    q2=q*q
+    q4=q2*q2
+    q6=q4*q2
+    q8=q4*q4
+    q9=q8*q
+    q10=q4*q6
+    q12=q6*q6
+    q14=q6*q8
+    q16=q8*q8
+    q18=q9*q9
+    q20=q10*q10
+    q22=q20*q2
+    q25=q16*q9
+    t3=1.d0+2.d0*(q+q4+q9+q16+q25) !&
+    t3sq=t3*t3
+    fac=piover6/t3sq
+    wq=q2+3.d0*q4+4.d0*q6+7.d0*q8+6.d0*q10+12.d0*q12+8.d0*q14  &
+      +15.d0*q16+13.d0*q18+18.d0*q20+12.d0*q22
+    ellipe=fac*(1.d0+(1.d0+kp2)*t3sq*t3sq-24.d0*wq)
+
+  else
+
+    kp2=k*k ! modify for large k
+    kp=sqrt(kp2)
+    sqrtkp=sqrt(kp)
+    eps=0.5d0*(1.d0-sqrtkp)/(1.d0+sqrtkp)
+    eps4=eps*eps*eps*eps
+    q=eps*(1.d0+eps4*(2.d0+eps4*(15.d0  &
+     +eps4*(150.d0+eps4*(1707.d0+eps4*(20910.d0+eps4*(268616.d0)))))))
+    q2=q*q
+    q4=q2*q2
+    q6=q4*q2
+    q8=q4*q4
+    q9=q8*q
+    q10=q4*q6
+    q12=q6*q6
+    q14=q6*q8
+    q16=q8*q8
+    q18=q9*q9
+    q20=q10*q10
+    q22=q20*q2
+    q25=q16*q9
+    t3=1.d0+2.d0*(q+q4+q9+q16+q25)
+    t3sq=t3*t3
+    fac=piover6/t3sq
+    wq=q2+3.d0*q4+4.d0*q6+7.d0*q8+6.d0*q10+12.d0*q12+8.d0*q14  &
+      +15.d0*q16+13.d0*q18+18.d0*q20+12.d0*q22
+    ellipep=fac*(1.d0+(1.d0+kp2)*t3sq*t3sq-24.d0*wq) ! modify for large k
+    ellipkp=halfpi*t3sq                              ! modify for large k
+    ellipk=-piinv*ellipkp*log(q)                     ! modify for large k
+    ellipe=ellipk+(halfpi-ellipep*ellipk)/ellipkp    ! modify for large k
+
+  end if
+
+    return
+
+END FUNCTION ellipe
 
 END MODULE GKV_math

run/Makefile

@@ -50,13 +50,15 @@ gkvp:	$(SRC)gkvp_header.f90\
 	$(SRC)gkvp_tips.f90\
 	$(SRC)gkvp_vmecbzx.f90\
 	$(SRC)gkvp_igs.f90\
+	$(SRC)gkvp_ring.f90\
 	$(SRC)gkvp_bndry.f90\
 	$(SRC)gkvp_colli.f90\
 	$(SRC)$(FFT).f90\
 	$(SRC)gkvp_fld.f90\
 	$(SRC)gkvp_colliimp.f90\
 	$(SRC)gkvp_freq.f90\
 	$(SRC)gkvp_zfilter.f90\
+	$(SRC)gkvp_geom.f90\
 	$(SRC)gkvp_exb.f90\
 	$(SRC)gkvp_trans.f90\
 	$(SRC)gkvp_advnc.f90\
@@ -75,13 +77,15 @@ gkvp:	$(SRC)gkvp_header.f90\
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_tips.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_vmecbzx.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_igs.f90
+	$(FC) $(FFLAGS) -c $(SRC)gkvp_ring.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_bndry.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_colli.f90
 	$(FC) $(FFLAGS) -c $(SRC)$(FFT).f90 $(INC)
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_fld.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_colliimp.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_freq.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_zfilter.f90
+	$(FC) $(FFLAGS) -c $(SRC)gkvp_geom.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_exb.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_trans.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_advnc.f90
@@ -101,13 +105,15 @@ gkvp:	$(SRC)gkvp_header.f90\
 			gkvp_tips.o\
 			gkvp_vmecbzx.o\
 			gkvp_igs.o\
+			gkvp_ring.o\
 			gkvp_bndry.o\
 			gkvp_colli.o\
 			$(FFT).o\
 			gkvp_fld.o\
 			gkvp_colliimp.o\
 			gkvp_freq.o\
 			gkvp_zfilter.o\
+			gkvp_geom.o\
 			gkvp_exb.o\
 			gkvp_trans.o\
 			gkvp_advnc.o\

run/backup/Makefile_flow

@@ -50,13 +50,15 @@ gkvp:	$(SRC)gkvp_header.f90\
 	$(SRC)gkvp_tips.f90\
 	$(SRC)gkvp_vmecbzx.f90\
 	$(SRC)gkvp_igs.f90\
+	$(SRC)gkvp_ring.f90\
 	$(SRC)gkvp_bndry.f90\
 	$(SRC)gkvp_colli.f90\
 	$(SRC)$(FFT).f90\
 	$(SRC)gkvp_fld.f90\
 	$(SRC)gkvp_colliimp.f90\
 	$(SRC)gkvp_freq.f90\
 	$(SRC)gkvp_zfilter.f90\
+	$(SRC)gkvp_geom.f90\
 	$(SRC)gkvp_exb.f90\
 	$(SRC)gkvp_trans.f90\
 	$(SRC)gkvp_advnc.f90\
@@ -75,13 +77,15 @@ gkvp:	$(SRC)gkvp_header.f90\
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_tips.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_vmecbzx.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_igs.f90
+	$(FC) $(FFLAGS) -c $(SRC)gkvp_ring.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_bndry.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_colli.f90
 	$(FC) $(FFLAGS) -c $(SRC)$(FFT).f90 $(INC)
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_fld.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_colliimp.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_freq.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_zfilter.f90
+	$(FC) $(FFLAGS) -c $(SRC)gkvp_geom.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_exb.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_trans.f90
 	$(FC) $(FFLAGS) -c $(SRC)gkvp_advnc.f90
@@ -101,13 +105,15 @@ gkvp:	$(SRC)gkvp_header.f90\
 			gkvp_tips.o\
 			gkvp_vmecbzx.o\
 			gkvp_igs.o\
+			gkvp_ring.o\
 			gkvp_bndry.o\
 			gkvp_colli.o\
 			$(FFT).o\
 			gkvp_fld.o\
 			gkvp_colliimp.o\
 			gkvp_freq.o\
 			gkvp_zfilter.o\
+			gkvp_geom.o\
 			gkvp_exb.o\
 			gkvp_trans.o\
 			gkvp_advnc.o\