Merge pull request #403 from gafusion/tglf_mxhgeo

Implementation of Miller eXtended Harmonic (MXH) flux surface parameterization in TGLF

profiles_gen/locpargen/locpargen_tglf.f90

@@ -22,6 +22,31 @@ subroutine locpargen_tglf
   write(1,10) 'S_DELTA_LOC=',s_delta_loc
   write(1,10) 'ZETA_LOC=',zeta_loc
   write(1,10) 'S_ZETA_LOC=',s_zeta_loc
+  write(1,*)
+  write(1,'(a)')  '# Geometry (Advanced)'
+  write(1,10) 'SHAPE_SIN3=',shape_sin3_loc
+  write(1,10) 'SHAPE_S_SIN3=',shape_s_sin3_loc
+  write(1,10) 'SHAPE_SIN4=',shape_sin4_loc
+  write(1,10) 'SHAPE_S_SIN4=',shape_s_sin4_loc
+  write(1,10) 'SHAPE_SIN5=',shape_sin5_loc
+  write(1,10) 'SHAPE_S_SIN5=',shape_s_sin5_loc
+  write(1,10) 'SHAPE_SIN6=',shape_sin6_loc
+  write(1,10) 'SHAPE_S_SIN6=',shape_s_sin6_loc
+  write(1,10) 'SHAPE_COS0=',shape_cos0_loc
+  write(1,10) 'SHAPE_S_COS0=',shape_s_cos0_loc
+  write(1,10) 'SHAPE_COS1=',shape_cos1_loc
+  write(1,10) 'SHAPE_S_COS1=',shape_s_cos1_loc
+  write(1,10) 'SHAPE_COS2=',shape_cos2_loc
+  write(1,10) 'SHAPE_S_COS2=',shape_s_cos2_loc
+  write(1,10) 'SHAPE_COS3=',shape_cos3_loc
+  write(1,10) 'SHAPE_S_COS3=',shape_s_cos3_loc 
+  write(1,10) 'SHAPE_COS4=',shape_cos4_loc
+  write(1,10) 'SHAPE_S_COS4=',shape_s_cos4_loc 
+  write(1,10) 'SHAPE_COS5=',shape_cos5_loc
+  write(1,10) 'SHAPE_S_COS5=',shape_s_cos5_loc 
+  write(1,10) 'SHAPE_COS6=',shape_cos6_loc
+  write(1,10) 'SHAPE_S_COS6=',shape_s_cos6_loc 
+  write(1,*)
   write(1,10) 'SIGN_BT=',btccw
   write(1,10) 'SIGN_IT=',ipccw
   write(1,*)

profiles_gen/locpargen/locpargen_tglf_stack.f90

@@ -17,7 +17,18 @@ subroutine locpargen_tglf_stack
        'ZMAJ_LOC=',zmag_loc, 'DZMAJDX_LOC=',dzmag_loc, 'Q_LOC=',abs(q_loc), &
        'Q_PRIME_LOC=',(q_loc/r0)**2*s_loc, 'KAPPA_LOC=',kappa_loc, 'S_KAPPA_LOC=',s_kappa_loc, &
        'DELTA_LOC=',delta_loc, 'S_DELTA_LOC=',s_delta_loc, 'ZETA_LOC=',zeta_loc, &
-       'S_ZETA_LOC=',s_zeta_loc, 'SIGN_BT=',btccw, 'SIGN_IT=',ipccw, &
+       'S_ZETA_LOC=',s_zeta_loc, 'SHAPE_SIN3=',shape_sin3_loc, 'SHAPE_S_SIN3=',shape_s_sin3_loc, &
+       'SHAPE_SIN4=' ,shape_sin4_loc, 'SHAPE_S_SIN4=' ,shape_s_sin4_loc, &
+       'SHAPE_SIN5=' ,shape_sin5_loc, 'SHAPE_S_SIN5=' ,shape_s_sin5_loc, &
+       'SHAPE_SIN6=' ,shape_sin6_loc, 'SHAPE_S_SIN6=' ,shape_s_sin6_loc, &
+       'SHAPE_COS0=' ,shape_cos0_loc, 'SHAPE_S_COS0=' ,shape_s_cos0_loc, &
+       'SHAPE_COS1=' ,shape_cos1_loc, 'SHAPE_S_COS1=' ,shape_s_cos1_loc, &
+       'SHAPE_COS2=' ,shape_cos2_loc, 'SHAPE_S_COS2=' ,shape_s_cos2_loc, &
+       'SHAPE_COS3=' ,shape_cos3_loc, 'SHAPE_S_COS3=' ,shape_s_cos3_loc, &
+       'SHAPE_COS4=' ,shape_cos4_loc, 'SHAPE_S_COS4=' ,shape_s_cos4_loc, &
+       'SHAPE_COS5=' ,shape_cos5_loc, 'SHAPE_S_COS5=' ,shape_s_cos5_loc, &
+       'SHAPE_COS6=' ,shape_cos6_loc, 'SHAPE_S_COS6=' ,shape_s_cos6_loc, &
+       'SIGN_BT=',btccw, 'SIGN_IT=',ipccw, &
        'ZEFF=',z_eff_loc, 'XNUE=',nu_ee*a/cs_loc, &
        'DEBYE=',7.43*sqrt(1e3*temp_loc(ise)/(1e13*dens_loc(ise)))/abs(rhos_loc), &
        'BETAE=',betae_unit, 'P_PRIME_LOC=',(abs(q_loc)/r0)*(-beta_star_loc/(8*pi)), &

tglf/bin/tglf_defaults.py

@@ -162,6 +162,28 @@ def set_defaults() :
      t.add('S_DELTA_LOC','0.0')
      t.add('ZETA_LOC','0.0')
      t.add('S_ZETA_LOC','0.0')
+     t.add('SHAPE_COS0','0.0')
+     t.add('SHAPE_S_COS0','0.0')
+     t.add('SHAPE_COS1','0.0')
+     t.add('SHAPE_S_COS1','0.0')
+     t.add('SHAPE_COS2','0.0')
+     t.add('SHAPE_S_COS2','0.0')
+     t.add('SHAPE_COS3','0.0')
+     t.add('SHAPE_S_COS3','0.0')
+     t.add('SHAPE_COS4','0.0')
+     t.add('SHAPE_S_COS4','0.0')
+     t.add('SHAPE_COS5','0.0')
+     t.add('SHAPE_S_COS5','0.0')
+     t.add('SHAPE_COS6','0.0')
+     t.add('SHAPE_S_COS6','0.0')
+     t.add('SHAPE_SIN3','0.0')
+     t.add('SHAPE_S_SIN3','0.0')
+     t.add('SHAPE_SIN4','0.0')
+     t.add('SHAPE_S_SIN4','0.0')
+     t.add('SHAPE_SIN5','0.0')
+     t.add('SHAPE_S_SIN5','0.0')
+     t.add('SHAPE_SIN6','0.0')
+     t.add('SHAPE_S_SIN6','0.0')
      t.add('Q_LOC','2.0')
      t.add('Q_PRIME_LOC','16.0')
      t.add('P_PRIME_LOC','0.0')

tglf/src/tglf_driver.f90

@@ -48,7 +48,11 @@ PROGRAM tglf_driver
       elseif(igeo_tg.eq.1)then
 !        q_prime_tg = shat_tg*(q_tg/rmin_tg)**2
         CALL put_Miller_geometry(rmin_tg,rmaj_tg,zmaj_tg,drmindx_tg,drmajdx_tg,dzmajdx_tg, &
-         kappa_tg,s_kappa_tg,delta_tg,s_delta_tg,zeta_tg,s_zeta_tg,q_tg,q_prime_tg,p_prime_tg,kx0_tg)
+         kappa_tg,s_kappa_tg,delta_tg,s_delta_tg,zeta_tg,s_zeta_tg,shape_sin3_tg,shape_s_sin3_tg,& 
+         shape_sin4_tg,shape_s_sin4_tg,shape_sin5_tg,shape_s_sin5_tg,shape_sin6_tg,shape_s_sin6_tg,&
+         shape_cos0_tg,shape_s_cos0_tg,shape_cos1_tg,shape_s_cos1_tg,shape_cos2_tg,shape_s_cos2_tg,&
+         shape_cos3_tg,shape_s_cos3_tg,shape_cos4_tg,shape_s_cos4_tg,shape_cos5_tg,shape_s_cos5_tg,&
+         shape_cos6_tg,shape_s_cos6_tg,q_tg,q_prime_tg,p_prime_tg,kx0_tg)
       elseif(igeo_tg.eq.2)then
         CALL put_Fourier_geometry(q_tg,q_prime_tg,p_prime_tg,nfourier_tg,fourier_tg)
       elseif(igeo_tg.eq.3)then

tglf/src/tglf_driver_mpi.f90

@@ -52,7 +52,11 @@ PROGRAM tglf_driver_mpi
       elseif(igeo_tg.eq.1)then
 !        q_prime_tg = shat_tg*(q_tg/rmin_tg)**2
         CALL put_Miller_geometry(rmin_tg,rmaj_tg,zmaj_tg,drmindx_tg,drmajdx_tg,dzmajdx_tg, &
-         kappa_tg,s_kappa_tg,delta_tg,s_delta_tg,zeta_tg,s_zeta_tg,q_tg,q_prime_tg,p_prime_tg,kx0_tg)
+         kappa_tg,s_kappa_tg,delta_tg,s_delta_tg,zeta_tg,s_zeta_tg,shape_sin3_tg,shape_s_sin3_tg,& 
+         shape_sin4_tg,shape_s_sin4_tg,shape_sin5_tg,shape_s_sin5_tg,shape_sin6_tg,shape_s_sin6_tg,&
+         shape_cos0_tg,shape_s_cos0_tg,shape_cos1_tg,shape_s_cos1_tg,shape_cos2_tg,shape_s_cos2_tg,&
+         shape_cos3_tg,shape_s_cos3_tg,shape_cos4_tg,shape_s_cos4_tg,shape_cos5_tg,shape_s_cos5_tg,&
+         shape_cos6_tg,shape_s_cos6_tg,q_tg,q_prime_tg,p_prime_tg,kx0_tg)
       elseif(igeo_tg.eq.2)then
         CALL put_Fourier_geometry(q_tg,q_prime_tg,p_prime_tg,nfourier_tg,fourier_tg)
       elseif(igeo_tg.eq.3)then

tglf/src/tglf_geometry.f90

@@ -1324,6 +1324,7 @@ END SUBROUTINE mercier_write
 ! squarness shear s_zeta_loc = rmin*d(zeta)/dx. Included elevation Zmaj_loc.
 ! Also changed definition of s_delta = rmin*d(delta)/dx from Waltz-Miller convention to GYRO's.
 ! July 26, 2021 set elevation Zmaj_loc =0.0 and DZMAJDX_LOC=0.0 since these break the up/down symmetry of Miller by contributing to Grad_r
+! August 1, 2024 Sophia Guizzo([email protected]) revised with MXH parameterization (Arbon et al PPCF 2020)
 !---------------------------------------------------------------
 
 SUBROUTINE miller_geo
@@ -1374,10 +1375,36 @@ SUBROUTINE miller_geo
   ! compute the arclength around the flux surface
   !
   theta = 0.0
-  arg_r = theta+x_delta*sin(theta)
-  darg_r = 1.0+x_delta*cos(theta)
-  arg_z = theta + zeta_loc*sin(2.0*theta)
-  darg_z = 1.0 + zeta_loc*2.0*cos(2.0*theta)
+  arg_r = theta + shape_cos0_loc + &
+  shape_cos1_loc*cos(theta) + &
+  shape_cos2_loc*cos(2*theta) + &
+  shape_cos3_loc*cos(3*theta) + &
+  shape_cos4_loc*cos(4*theta) + &
+  shape_cos5_loc*cos(5*theta) + &
+  shape_cos6_loc*cos(6*theta) + &
+  x_delta*sin(theta)- &
+  zeta_loc*sin(2*theta) + &
+  shape_sin3_loc*sin(3*theta) +&
+  shape_sin4_loc*sin(4*theta) + &
+  shape_sin5_loc*sin(5*theta) +&
+  shape_sin6_loc*sin(6*theta)
+
+  darg_r = 1 - shape_cos1_loc*sin(theta) - &
+  2*shape_cos2_loc*sin(2*theta) - &
+  3*shape_cos3_loc*sin(3*theta) - &
+  4*shape_cos4_loc*sin(4*theta) - &
+  5*shape_cos5_loc*sin(5*theta) - &
+  6*shape_cos6_loc*sin(6*theta) + &
+  x_delta*cos(theta)-&
+  2*zeta_loc*cos(2*theta) + &
+  3*shape_sin3_loc*cos(3*theta) + &
+  4*shape_sin4_loc*cos(4*theta) + &
+  5*shape_sin5_loc*cos(5*theta) + &
+  6*shape_sin6_loc*cos(6*theta)
+
+  arg_z = theta
+  darg_z = 1.0
+
   r_t = -rmin_loc*sin(arg_r)*darg_r
   z_t = kappa_loc*rmin_loc*cos(arg_z)*darg_z
   l_t = SQRT(r_t**2+z_t**2)
@@ -1387,6 +1414,7 @@ SUBROUTINE miller_geo
   l_t1 = l_t
   scale_max=l_t
   arclength = 0.0
+
   do while(theta.lt.pi_2)
      theta = theta + dtheta
      if(theta.gt.pi_2)then
@@ -1395,15 +1423,39 @@ SUBROUTINE miller_geo
         theta = pi_2
      endif
      ! write(*,*)"theta = ",theta,"dtheta=",dtheta
-     arg_r = theta+x_delta*sin(theta)
+     arg_r = theta + shape_cos0_loc + &
+     shape_cos1_loc*cos(theta) + &
+     shape_cos2_loc*cos(2*theta) + &
+     shape_cos3_loc*cos(3*theta) + &
+     shape_cos4_loc*cos(4*theta) + &
+     shape_cos5_loc*cos(5*theta) + &
+     shape_cos6_loc*cos(6*theta) + &
+     x_delta*sin(theta)- &
+     zeta_loc*sin(2*theta) + &
+     shape_sin3_loc*sin(3*theta) +&
+     shape_sin4_loc*sin(4*theta) + &
+     shape_sin5_loc*sin(5*theta) +&
+     shape_sin6_loc*sin(6*theta)
+
      ! d(arg_r)/dtheta
-     darg_r = 1.0+x_delta*cos(theta)
+     darg_r = 1 - shape_cos1_loc*sin(theta) - &
+     2*shape_cos2_loc*sin(2*theta) - &
+     3*shape_cos3_loc*sin(3*theta) - &
+     4*shape_cos4_loc*sin(4*theta) - &
+     5*shape_cos5_loc*sin(5*theta) - &
+     6*shape_cos6_loc*sin(6*theta) + &
+     x_delta*cos(theta)-&
+     2*zeta_loc*cos(2*theta) + &
+     3*shape_sin3_loc*cos(3*theta) + &
+     4*shape_sin4_loc*cos(4*theta) + &
+     5*shape_sin5_loc*cos(5*theta) + &
+     6*shape_sin6_loc*cos(6*theta)
      ! dR/dtheta
      r_t = -rmin_loc*sin(arg_r)*darg_r
      !
-     arg_z = theta + zeta_loc*sin(2.0*theta)
+     arg_z = theta
      ! d(arg_z)/dtheta
-     darg_z = 1.0 + zeta_loc*2.0*cos(2.0*theta)
+     darg_z = 1.0 
      ! dZ/dtheta
      z_t = kappa_loc*rmin_loc*cos(arg_z)*darg_z
      ! dl/dtheta
@@ -1431,22 +1483,68 @@ SUBROUTINE miller_geo
   t_s(ms)=-pi_2
   ! make a first guess based on theta=0.0
   theta=0.0
-  arg_r = theta+x_delta*sin(theta)
-  darg_r = 1.0+x_delta*cos(theta)
-  arg_z = theta + zeta_loc*sin(2.0*theta)
-  darg_z = 1.0 + zeta_loc*2.0*cos(2.0*theta)
+  arg_r = theta + shape_cos0_loc + &
+  shape_cos1_loc*cos(theta) + &
+  shape_cos2_loc*cos(2*theta) + &
+  shape_cos3_loc*cos(3*theta) + &
+  shape_cos4_loc*cos(4*theta) + &
+  shape_cos5_loc*cos(5*theta) + &
+  shape_cos6_loc*cos(6*theta) + &
+  x_delta*sin(theta)- &
+  zeta_loc*sin(2*theta) + &
+  shape_sin3_loc*sin(3*theta) +&
+  shape_sin4_loc*sin(4*theta) + &
+  shape_sin5_loc*sin(5*theta) +&
+  shape_sin6_loc*sin(6*theta)
+  darg_r = 1 - shape_cos1_loc*sin(theta) - &
+  2*shape_cos2_loc*sin(2*theta) - &
+  3*shape_cos3_loc*sin(3*theta) - &
+  4*shape_cos4_loc*sin(4*theta) - &
+  5*shape_cos5_loc*sin(5*theta) - &
+  6*shape_cos6_loc*sin(6*theta) + &
+  x_delta*cos(theta)-&
+  2*zeta_loc*cos(2*theta) + &
+  3*shape_sin3_loc*cos(3*theta) + &
+  4*shape_sin4_loc*cos(4*theta) + &
+  5*shape_sin5_loc*cos(5*theta) + &
+  6*shape_sin6_loc*cos(6*theta)
+  arg_z = theta
+  darg_z = 1.0 
   r_t = -rmin_loc*sin(arg_r)*darg_r
   z_t = kappa_loc*rmin_loc*cos(arg_z)*darg_z
   l_t = SQRT(r_t**2+z_t**2)
   dtheta = -ds/l_t
   theta=dtheta
   l_t1=l_t
   !
-  do m=1,ms/2
-     arg_r = theta+x_delta*sin(theta)
-     darg_r = 1.0+x_delta*cos(theta)
-     arg_z = theta + zeta_loc*sin(2.0*theta)
-     darg_z = 1.0 + zeta_loc*2.0*cos(2.0*theta)
+  do m=1,ms
+     arg_r = theta + shape_cos0_loc + &
+     shape_cos1_loc*cos(theta) + &
+     shape_cos2_loc*cos(2*theta) + &
+     shape_cos3_loc*cos(3*theta) + &
+     shape_cos4_loc*cos(4*theta) + &
+     shape_cos5_loc*cos(5*theta) + &
+     shape_cos6_loc*cos(6*theta) + &
+     x_delta*sin(theta)- &
+     zeta_loc*sin(2*theta) + &
+     shape_sin3_loc*sin(3*theta) +&
+     shape_sin4_loc*sin(4*theta) + &
+     shape_sin5_loc*sin(5*theta) +&
+     shape_sin6_loc*sin(6*theta)
+     darg_r = 1 - shape_cos1_loc*sin(theta) - &
+     2*shape_cos2_loc*sin(2*theta) - &
+     3*shape_cos3_loc*sin(3*theta) - &
+     4*shape_cos4_loc*sin(4*theta) - &
+     5*shape_cos5_loc*sin(5*theta) - &
+     6*shape_cos6_loc*sin(6*theta) + &
+     x_delta*cos(theta)-&
+     2*zeta_loc*cos(2*theta) + &
+     3*shape_sin3_loc*cos(3*theta) + &
+     4*shape_sin4_loc*cos(4*theta) + &
+     5*shape_sin5_loc*cos(5*theta) + &
+     6*shape_sin6_loc*cos(6*theta)
+     arg_z = theta 
+     darg_z = 1.0 
      r_t = -rmin_loc*sin(arg_r)*darg_r
      z_t = kappa_loc*rmin_loc*cos(arg_z)*darg_z
      l_t = SQRT(r_t**2+z_t**2)
@@ -1456,13 +1554,12 @@ SUBROUTINE miller_geo
      l_t1=l_t
   enddo
   ! distribute endpoint error over interior points
-  dtheta = (t_s(ms/2)-(-pi))/REAL(ms/2)
+  dtheta = (t_s(ms)-(-pi_2))/REAL(ms)
   ! write(*,*)"enpoint error =",dtheta
   ! dtheta=0.0
   ! t_s(ms/2)=-pi
-  do m=1,ms/2
+  do m=1,ms
      t_s(m) = t_s(m)-REAL(m)*dtheta
-     t_s(ms-m)=-pi_2 - t_s(m)
   enddo
   ! write(*,*)"t_s(ms/2)+pi=",t_s(ms/2)+pi
   !
@@ -1491,10 +1588,33 @@ SUBROUTINE miller_geo
   do m=0,ms
 
      theta = t_s(m)
-     arg_r = theta + x_delta*sin(theta)
-     darg_r = 1.0 + x_delta*cos(theta)
-     arg_z = theta + zeta_loc*sin(2.0*theta)
-     darg_z = 1.0 + zeta_loc*2.0*cos(2.0*theta)
+     arg_r = theta + shape_cos0_loc + &
+     shape_cos1_loc*cos(theta) + &
+     shape_cos2_loc*cos(2*theta) + &
+     shape_cos3_loc*cos(3*theta) + &
+     shape_cos4_loc*cos(4*theta) + &
+     shape_cos5_loc*cos(5*theta) + &
+     shape_cos6_loc*cos(6*theta) + &
+     x_delta*sin(theta)- &
+     zeta_loc*sin(2*theta) + &
+     shape_sin3_loc*sin(3*theta) +&
+     shape_sin4_loc*sin(4*theta) + &
+     shape_sin5_loc*sin(5*theta) +&
+     shape_sin6_loc*sin(6*theta)
+     darg_r = 1 - shape_cos1_loc*sin(theta) - &
+     2*shape_cos2_loc*sin(2*theta) - &
+     3*shape_cos3_loc*sin(3*theta) - &
+     4*shape_cos4_loc*sin(4*theta) - &
+     5*shape_cos5_loc*sin(5*theta) - &
+     6*shape_cos6_loc*sin(6*theta) + &
+     x_delta*cos(theta)-&
+     2*zeta_loc*cos(2*theta) + &
+     3*shape_sin3_loc*cos(3*theta) + &
+     4*shape_sin4_loc*cos(4*theta) + &
+     5*shape_sin5_loc*cos(5*theta) + &
+     6*shape_sin6_loc*cos(6*theta)
+     arg_z = theta 
+     darg_z = 1.0 
 
      ! R(theta)
      ! Z(theta)
@@ -1514,9 +1634,21 @@ SUBROUTINE miller_geo
      ! dR/dr
      ! dZ/dr
      R_r = drmajdx_loc + drmindx_loc*cos(arg_r) &
-          -sin(arg_r)*s_delta_loc*sin(theta)/sqrt(1.0 - delta_loc**2)
-     Z_r = dzmajdx_loc + kappa_loc*sin(arg_z)*(drmindx_loc +s_kappa_loc) &
-          +kappa_loc*cos(arg_z)*s_zeta_loc*sin(2.0*theta)
+          -sin(arg_r)*(shape_s_cos0_loc + &
+          shape_s_cos1_loc*cos(theta) + &
+          shape_s_cos2_loc*cos(2*theta) + &
+          shape_s_cos3_loc*cos(3*theta) + &
+          shape_s_cos4_loc*cos(4*theta) + &
+          shape_s_cos5_loc*cos(5*theta) + &
+          shape_s_cos6_loc*cos(6*theta) + &
+          s_delta_loc*sin(theta)/cos(x_delta) - &
+          s_zeta_loc*sin(2*theta) + &
+          shape_s_sin3_loc*sin(3*theta) + &
+          shape_s_sin4_loc*sin(4*theta) + &
+          shape_s_sin5_loc*sin(5*theta) + &
+          shape_s_sin6_loc*sin(6*theta))
+
+     Z_r = dzmajdx_loc + kappa_loc*sin(arg_z)*(drmindx_loc +s_kappa_loc) 
      ! Jacobian
      det = R_r*z_t - R_t*Z_r
      ! grad_r