diff --git a/bash scripts/fort_00.90 b/bash scripts/fort_00.90
new file mode 100644
index 0000000000000000000000000000000000000000..183029f992e968409f6412c35caefa64e93fcf3a
--- /dev/null
+++ b/bash scripts/fort_00.90
@@ -0,0 +1,107 @@
+&BASIC
+ nrun = 100000000
+ dt = 0.002
+ tmax = 100
+ maxruntime = 356400
+ job2load = -1
+/
+&GRID
+ pmax = 2
+ jmax = 1
+ Nx = 192
+ Lx = 90
+ Ny = 48
+ Ly = 200
+ Nz = 32
+ SG = .false.
+ Nexc = 0
+/
+&GEOMETRY
+ geom = 'miller'
+ q0 =-2.15
+ shear = 3.62
+ eps = 0.28
+ kappa = 1.53
+ s_kappa = 0.77
+ delta = 0.23
+ s_delta = 1.05
+ zeta =-0.01
+ s_zeta =-0.17
+ parallel_bc = 'dirichlet'
+ shift_y = 0
+ Npol = 1
+ PB_PHASE= .false.
+/
+&OUTPUT_PAR
+ dtsave_0d = 0.5
+ dtsave_1d = -1
+ dtsave_2d = -1
+ dtsave_3d = 5.0
+ dtsave_5d = 10
+ write_doubleprecision = .false.
+ write_gamma = .true.
+ write_hf = .true.
+ write_phi = .true.
+ write_Na00 = .true.
+ write_Napj = .true.
+ write_dens = .true.
+ write_temp = .true.
+/
+&MODEL_PAR
+LINEARITY = 'nonlinear'
+RM_LD_T_EQ= .false.
+ Na = 2
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 2.0
+ HYP_V = 'hypcoll'
+ mu_p = 0
+ mu_j = 0
+ nu = 0.5
+ k_gB = 1
+ k_cB = 1
+ lambdaD = 0
+ beta = 0.0034
+ ADIAB_E = .false.
+ ADIAB_I = .false.
+ MHD_PD = .true.
+/
+&CLOSURE_PAR
+ hierarchy_closure='truncation'
+ dmax =-1
+ nonlinear_closure='truncation'
+ nmax =0
+/
+&SPECIES
+ name_ = 'ions'
+ tau_ = 1
+ sigma_ = 1
+ q_ = 1
+ K_N_ = 1.33
+ K_T_ = 8.25
+/
+&SPECIES
+ name_ = 'electrons'
+ tau_ = 1
+ sigma_ = 0.04!0.023338
+ q_ = -1
+ K_N_ = 1.33
+ K_T_ = 12
+/
+&COLLISION_PAR
+ collision_model = 'DG'
+ GK_CO = .true.
+ INTERSPECIES = .true.
+ mat_file = '/Users/ahoffmann/gyacomo/iCa/null'
+ collision_kcut = 1
+/
+&INITIAL_CON
+ INIT_OPT = 'blob'
+ init_background = 0
+ init_noiselvl = 1e-05
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
\ No newline at end of file
diff --git a/bash scripts/submit_00.cmd b/bash scripts/submit_00.cmd
index 480ee2b3ffacc6b1d4373c1b1aa964c17e70e725..0d1f9436f14f046c8a0198ab88b73282e77d02ee 100644
--- a/bash scripts/submit_00.cmd
+++ b/bash scripts/submit_00.cmd
@@ -10,5 +10,4 @@
#SBATCH --account=FUA36_TSVVT422
#SBATCH --partition=skl_fua_dbg
#SBATCH --qos=normal
-srun --cpu-bind=cores ./gyacomo 2 24 1 0
-
+srun --cpu-bind=cores ./gyacomo 2 24 1 0
\ No newline at end of file
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 5675ced44c526cf619d6366e212b517f3c18cff1..cf3adae4259d686581839b48b0c03b3cabde3af2 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -18,6 +18,8 @@ implicit none
REAL(xp), PUBLIC, PROTECTED :: s_delta = 0._xp ! '' sdelta = r/sqrt(1-delta2) ddelta/dr
REAL(xp), PUBLIC, PROTECTED :: zeta = 0._xp ! squareness
REAL(xp), PUBLIC, PROTECTED :: s_zeta = 0._xp ! '' szeta = r dzeta/dr
+ REAL(xp), PUBLIC, PROTECTED :: theta1 = 0._xp ! for Miller global
+ REAL(xp), PUBLIC, PROTECTED :: theta2 = 0._xp !
! to apply shift in the parallel z-BC if shearless
REAL(xp), PUBLIC, PROTECTED :: shift_y = 0._xp ! for Arno <3
REAL(xp), PUBLIC, PROTECTED :: Npol = 1._xp ! number of poloidal turns (real for 3D zpinch studies)
@@ -80,7 +82,8 @@ CONTAINS
! Read the input parameters
IMPLICIT NONE
NAMELIST /GEOMETRY/ geom, q0, shear, eps,&
- kappa, s_kappa,delta, s_delta, zeta, s_zeta,& ! For miller
+ kappa, s_kappa,delta, s_delta, zeta, s_zeta,& ! For Miller
+ theta1, theta2,& ! for Miller global
parallel_bc, shift_y, Npol, PB_PHASE
READ(lu_in,geometry)
PB_PHASE = .false.
@@ -139,15 +142,15 @@ CONTAINS
kappa = 1._xp
C_y = 1._xp
Cyq0_x0 = 1._xp
- CASE('miller','Miller')
+ CASE('miller','Miller','Miller_global','miller_global')
CALL speak('Miller geometry')
IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for Miller geometry"
IF(MODULO(FLOOR(Npol),2) .EQ. 0) THEN
write(*,*) "Npol must be odd for Miller, (Npol = ",Npol,")"
ERROR STOP
ENDIF
- call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta)
- call get_miller(eps,major_R,major_Z,q0,shear,FLOOR(Npol),alpha_MHD,edge_opt,&
+ call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta,theta1,theta2)
+ call get_miller(geom,eps,major_R,major_Z,q0,shear,FLOOR(Npol),alpha_MHD,edge_opt,&
C_y,C_xy,Cyq0_x0,xpdx_pm_geom,gxx,gxy,gxz,gyy,gyz,gzz,&
dBdx,dBdy,dBdz,hatB,jacobian,hatR,hatZ,dxdR,dxdZ)
CASE('circular','circ')
diff --git a/src/miller_mod.F90 b/src/miller_mod.F90
index 5f8528439cf351faff7586a7cce0a089272b0179..705675f876b2f74334e8b0a149603c1b4178f2a4 100644
--- a/src/miller_mod.F90
+++ b/src/miller_mod.F90
@@ -21,11 +21,12 @@ MODULE miller
real(xp) :: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta
real(xp) :: thetaShift
real(xp) :: thetak, thetad
+ real(xp) :: asurf, theta1, theta2, theta3, delta2, delta3, Raxis, Zaxis
CONTAINS
!>Set defaults for miller parameters
- subroutine set_miller_parameters(kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_)
- real(xp), INTENT(IN) :: kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_
+ subroutine set_miller_parameters(kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_,theta1_,theta2_)
+ real(xp), INTENT(IN) :: kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_,theta1_,theta2_
rho = -1._xp
kappa = kappa_
s_kappa = s_kappa_
@@ -36,15 +37,24 @@ CONTAINS
drR = 0._xp
drZ = 0._xp
thetak = 0._xp
- thetad = 0._xp
+ thetad = 0._xp
+ asurf = 0.54_xp
+ theta1 = theta1_
+ theta2 = theta2_
+ theta3 = 0._xp
+ delta2 = 1._xp
+ delta3 = 1._xp
+ Raxis = 1._xp
+ Zaxis = 0._xp
end subroutine set_miller_parameters
!>Get Miller metric, magnetic field, jacobian etc.
- subroutine get_miller(trpeps,major_R,major_Z,q0,shat,Npol,amhd,edge_opt,&
+ subroutine get_miller(geom,trpeps,major_R,major_Z,q0,shat,Npol,amhd,edge_opt,&
C_y,C_xy,Cyq0_x0,xpdx_pm_geom,gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,dBdx_,dBdy_,dBdz_,&
Bfield_,jacobian_,R_hat_,Z_hat_,dxdR_,dxdZ_)
!!!!!!!!!!!!!!!! GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ CHARACTER(len=16), INTENT(IN) :: geom
real(xp), INTENT(INOUT) :: trpeps ! eps in gyacomo (inverse aspect ratio)
real(xp), INTENT(INOUT) :: major_R ! major radius
real(xp), INTENT(INOUT) :: major_Z ! major Z
@@ -73,8 +83,8 @@ CONTAINS
real(xp), dimension(500*(Npol+2)):: dtdtZcirc, dtdtZelong, dtdtZelongTilt, dtdtZtri, dtdtZtriTilt, dtRcirc, dtRelong
real(xp), dimension(500*(Npol+2)):: dtRelongTilt, dtRtri, dtRtriTilt, dtZcirc, dtZelong, dtZelongTilt, dtZtri, dtZtriTilt
real(xp), dimension(500*(Npol+2)):: Rcirc, Relong, RelongTilt, Rtri, RtriTilt, Zcirc, Zelong, ZelongTilt, Ztri, ZtriTilt
- ! real(xp), dimension(500*(Npol+2)):: drrShape, drrAng, drxAng, dryAng, dtdtrShape, dtdtrAng, dtdtxAng
- ! real(xp), dimension(500*(Npol+2)):: dtdtyAng, dtrShape, dtrAng, dtxAng, dtyAng, rShape, rAng, xAng, yAng
+ real(xp), dimension(500*(Npol+2)):: drrShape, drrAng, drxAng, dryAng, dtdtrShape, dtdtrAng, dtdtxAng
+ real(xp), dimension(500*(Npol+2)):: dtdtyAng, dtrShape, dtrAng, dtxAng, dtyAng, rShape, rAng, xAng, yAng
real(xp), dimension(500*(Npol+2)):: theta, thAdj, J_r, B, Bp, D0, D1, D2, D3, nu, chi, psi1, nu1
real(xp), dimension(500*(Npol+2)):: tmp_reverse, theta_reverse, tmp_arr
@@ -91,9 +101,9 @@ CONTAINS
real(xp), dimension(1:total_nz+ngz):: gxx_out,gxy_out,gxz_out,gyy_out,gyz_out,gzz_out,Bfield_out,jacobian_out, dBdx_out, dBdz_out, chi_out
real(xp), dimension(1:total_nz+ngz):: R_out, Z_out, dxdR_out, dxdZ_out
real(xp):: d_inv, drPsi, dxPsi, dq_dx, dq_dpsi, R0, Z0, B0, F, D0_full, D1_full, D2_full, D3_full
- !real(xp) :: Lnorm, Psi0 ! currently module-wide defined anyway
+ real(xp) :: Lnorm, Psi0 ! currently module-wide defined anyway
real(xp):: pprime, ffprime, D0_mid, D1_mid, D2_mid, D3_mid, dx_drho, pi, mu_0, dzprimedz
- ! real(xp):: rho_a, psiN, drpsiN, CN2, CN3, Rcenter, Zcenter, drRcenter, drZcenter
+ real(xp):: rho_a, psiN, drpsiN, CN2, CN3, Rcenter, Zcenter, drRcenter, drZcenter
logical:: bMaxShift
integer:: i, k, iBmax
@@ -124,76 +134,156 @@ CONTAINS
do while (bMaxShift)
!flux surface parametrization
thAdj = theta + thetaShift
- if (zeta/=0._xp .or. s_zeta/=0._xp) then
- R = R0 + rho*Cos(thAdj + d_inv*Sin(thAdj))
- Z = Z0 + kappa*rho*Sin(thAdj + zeta*Sin(2*thAdj))
-
- R_rho = drR + Cos(thAdj + d_inv*Sin(thAdj)) - s_delta*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
- Z_rho = drZ + kappa*s_zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
- + kappa*Sin(thAdj + zeta*Sin(2*thAdj)) + kappa*s_kappa*Sin(thAdj + zeta*Sin(2*thAdj))
-
- R_theta = -(rho*(1 + d_inv*Cos(thAdj))*Sin(thAdj + d_inv*Sin(thAdj)))
- Z_theta = kappa*rho*(1 + 2*zeta*Cos(2*thAdj))*Cos(thAdj + zeta*Sin(2*thAdj))
-
- R_theta_theta = -(rho*(1 + d_inv*Cos(thAdj))**2*Cos(thAdj + d_inv*Sin(thAdj))) &
- + d_inv*rho*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
- Z_theta_theta = -4*kappa*rho*zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
- - kappa*rho*(1 + 2*zeta*Cos(2*thAdj))**2*Sin(thAdj + zeta*Sin(2*thAdj))
- else
- Rcirc = rho*Cos(thAdj - thetad + thetak)
- Zcirc = rho*Sin(thAdj - thetad + thetak)
- Relong = Rcirc
- Zelong = Zcirc + (-1 + kappa)*rho*Sin(thAdj - thetad + thetak)
- RelongTilt = Relong*Cos(thetad - thetak) - Zelong*Sin(thetad - thetak)
- ZelongTilt = Zelong*Cos(thetad - thetak) + Relong*Sin(thetad - thetak)
- Rtri = RelongTilt - rho*Cos(thAdj) + rho*Cos(thAdj + delta*Sin(thAdj))
- Ztri = ZelongTilt
- RtriTilt = Rtri*Cos(thetad) + Ztri*Sin(thetad)
- ZtriTilt = Ztri*Cos(thetad) - Rtri*Sin(thetad)
- R = R0 + RtriTilt
- Z = Z0 + ZtriTilt
-
- drRcirc = Cos(thAdj - thetad + thetak)
- drZcirc = Sin(thAdj - thetad + thetak)
- drRelong = drRcirc
- drZelong = drZcirc - (1 - kappa - kappa*s_kappa)*Sin(thAdj - thetad + thetak)
- drRelongTilt = drRelong*Cos(thetad - thetak) - drZelong*Sin(thetad - thetak)
- drZelongTilt = drZelong*Cos(thetad - thetak) + drRelong*Sin(thetad - thetak)
- drRtri = drRelongTilt - Cos(thAdj) + Cos(thAdj + delta*Sin(thAdj)) &
- - s_delta*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj))
- drZtri = drZelongTilt
- drRtriTilt = drRtri*Cos(thetad) + drZtri*Sin(thetad)
- drZtriTilt = drZtri*Cos(thetad) - drRtri*Sin(thetad)
- R_rho = drR + drRtriTilt
- Z_rho = drZ + drZtriTilt
-
- dtRcirc = -(rho*Sin(thAdj - thetad + thetak))
- dtZcirc = rho*Cos(thAdj - thetad + thetak)
- dtRelong = dtRcirc
- dtZelong = dtZcirc + (-1 + kappa)*rho*Cos(thAdj - thetad + thetak)
- dtRelongTilt = dtRelong*Cos(thetad - thetak) - dtZelong*Sin(thetad - thetak)
- dtZelongTilt = dtZelong*Cos(thetad - thetak) + dtRelong*Sin(thetad - thetak)
- dtRtri = dtRelongTilt + rho*Sin(thAdj) - rho*(1 + delta*Cos(thAdj))*Sin(thAdj + delta*Sin(thAdj))
- dtZtri = dtZelongTilt
- dtRtriTilt = dtRtri*Cos(thetad) + dtZtri*Sin(thetad)
- dtZtriTilt = dtZtri*Cos(thetad) - dtRtri*Sin(thetad)
- R_theta = dtRtriTilt
- Z_theta = dtZtriTilt
-
- dtdtRcirc = -(rho*Cos(thAdj - thetad + thetak))
- dtdtZcirc = -(rho*Sin(thAdj - thetad + thetak))
- dtdtRelong = dtdtRcirc
- dtdtZelong = dtdtZcirc - (-1 + kappa)*rho*Sin(thAdj - thetad + thetak)
- dtdtRelongTilt = dtdtRelong*Cos(thetad - thetak) - dtdtZelong*Sin(thetad - thetak)
- dtdtZelongTilt = dtdtZelong*Cos(thetad - thetak) + dtdtRelong*Sin(thetad - thetak)
- dtdtRtri = dtdtRelongTilt + rho*Cos(thAdj) - rho*(1 + delta*Cos(thAdj))**2*Cos(thAdj + delta*Sin(thAdj)) &
- + delta*rho*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj))
- dtdtZtri = dtdtZelongTilt
- dtdtRtriTilt = dtdtRtri*Cos(thetad) + dtdtZtri*Sin(thetad)
- dtdtZtriTilt = dtdtZtri*Cos(thetad) - dtdtRtri*Sin(thetad)
- R_theta_theta = dtdtRtriTilt
- Z_theta_theta = dtdtZtriTilt
- endif
+ SELECT CASE (geom)
+ CASE('Miller','miller')
+ if (zeta/=0._xp .or. s_zeta/=0._xp) then
+ R = R0 + rho*Cos(thAdj + d_inv*Sin(thAdj))
+ Z = Z0 + kappa*rho*Sin(thAdj + zeta*Sin(2*thAdj))
+
+ R_rho = drR + Cos(thAdj + d_inv*Sin(thAdj)) - s_delta*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
+ Z_rho = drZ + kappa*s_zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
+ + kappa*Sin(thAdj + zeta*Sin(2*thAdj)) + kappa*s_kappa*Sin(thAdj + zeta*Sin(2*thAdj))
+
+ R_theta = -(rho*(1 + d_inv*Cos(thAdj))*Sin(thAdj + d_inv*Sin(thAdj)))
+ Z_theta = kappa*rho*(1 + 2*zeta*Cos(2*thAdj))*Cos(thAdj + zeta*Sin(2*thAdj))
+
+ R_theta_theta = -(rho*(1 + d_inv*Cos(thAdj))**2*Cos(thAdj + d_inv*Sin(thAdj))) &
+ + d_inv*rho*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
+ Z_theta_theta = -4*kappa*rho*zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
+ - kappa*rho*(1 + 2*zeta*Cos(2*thAdj))**2*Sin(thAdj + zeta*Sin(2*thAdj))
+ else
+ Rcirc = rho*Cos(thAdj - thetad + thetak)
+ Zcirc = rho*Sin(thAdj - thetad + thetak)
+ Relong = Rcirc
+ Zelong = Zcirc + (-1 + kappa)*rho*Sin(thAdj - thetad + thetak)
+ RelongTilt = Relong*Cos(thetad - thetak) - Zelong*Sin(thetad - thetak)
+ ZelongTilt = Zelong*Cos(thetad - thetak) + Relong*Sin(thetad - thetak)
+ Rtri = RelongTilt - rho*Cos(thAdj) + rho*Cos(thAdj + delta*Sin(thAdj))
+ Ztri = ZelongTilt
+ RtriTilt = Rtri*Cos(thetad) + Ztri*Sin(thetad)
+ ZtriTilt = Ztri*Cos(thetad) - Rtri*Sin(thetad)
+ R = R0 + RtriTilt
+ Z = Z0 + ZtriTilt
+
+ drRcirc = Cos(thAdj - thetad + thetak)
+ drZcirc = Sin(thAdj - thetad + thetak)
+ drRelong = drRcirc
+ drZelong = drZcirc - (1 - kappa - kappa*s_kappa)*Sin(thAdj - thetad + thetak)
+ drRelongTilt = drRelong*Cos(thetad - thetak) - drZelong*Sin(thetad - thetak)
+ drZelongTilt = drZelong*Cos(thetad - thetak) + drRelong*Sin(thetad - thetak)
+ drRtri = drRelongTilt - Cos(thAdj) + Cos(thAdj + delta*Sin(thAdj)) &
+ - s_delta*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj))
+ drZtri = drZelongTilt
+ drRtriTilt = drRtri*Cos(thetad) + drZtri*Sin(thetad)
+ drZtriTilt = drZtri*Cos(thetad) - drRtri*Sin(thetad)
+ R_rho = drR + drRtriTilt
+ Z_rho = drZ + drZtriTilt
+
+ dtRcirc = -(rho*Sin(thAdj - thetad + thetak))
+ dtZcirc = rho*Cos(thAdj - thetad + thetak)
+ dtRelong = dtRcirc
+ dtZelong = dtZcirc + (-1 + kappa)*rho*Cos(thAdj - thetad + thetak)
+ dtRelongTilt = dtRelong*Cos(thetad - thetak) - dtZelong*Sin(thetad - thetak)
+ dtZelongTilt = dtZelong*Cos(thetad - thetak) + dtRelong*Sin(thetad - thetak)
+ dtRtri = dtRelongTilt + rho*Sin(thAdj) - rho*(1 + delta*Cos(thAdj))*Sin(thAdj + delta*Sin(thAdj))
+ dtZtri = dtZelongTilt
+ dtRtriTilt = dtRtri*Cos(thetad) + dtZtri*Sin(thetad)
+ dtZtriTilt = dtZtri*Cos(thetad) - dtRtri*Sin(thetad)
+ R_theta = dtRtriTilt
+ Z_theta = dtZtriTilt
+
+ dtdtRcirc = -(rho*Cos(thAdj - thetad + thetak))
+ dtdtZcirc = -(rho*Sin(thAdj - thetad + thetak))
+ dtdtRelong = dtdtRcirc
+ dtdtZelong = dtdtZcirc - (-1 + kappa)*rho*Sin(thAdj - thetad + thetak)
+ dtdtRelongTilt = dtdtRelong*Cos(thetad - thetak) - dtdtZelong*Sin(thetad - thetak)
+ dtdtZelongTilt = dtdtZelong*Cos(thetad - thetak) + dtdtRelong*Sin(thetad - thetak)
+ dtdtRtri = dtdtRelongTilt + rho*Cos(thAdj) - rho*(1 + delta*Cos(thAdj))**2*Cos(thAdj + delta*Sin(thAdj)) &
+ + delta*rho*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj))
+ dtdtZtri = dtdtZelongTilt
+ dtdtRtriTilt = dtdtRtri*Cos(thetad) + dtdtZtri*Sin(thetad)
+ dtdtZtriTilt = dtdtZtri*Cos(thetad) - dtdtRtri*Sin(thetad)
+ R_theta_theta = dtdtRtriTilt
+ Z_theta_theta = dtdtZtriTilt
+ endif
+ CASE('Miller_global','miller_global')
+ rho_a = rho/aSurf
+
+ CN2 = (-1._xp + Delta2**2)/(1._xp + Delta2**2)
+ CN3 = (-1._xp + Delta3**2)/(1._xp + Delta3**3)
+
+ psiN = rho_a**2 + CN2*rho_a**2 + CN3*rho_a**3
+ drpsiN = 2._xp*rho_a + 2._xp*CN2*rho_a + 3._xp*CN3*rho_a**2
+
+ xAng = 2*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**3 - 27._xp*CN3**2*Cos(3._xp*(thAdj + theta3))**2*psiN
+ drxAng = -27._xp*CN3**2*Cos(3._xp*(thAdj + theta3))**2*drpsiN
+ dtxAng = -12._xp*CN2*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**2*Sin(2._xp*(thAdj + theta2)) &
+ + 162._xp*CN3**2*Cos(3._xp*(thAdj + theta3))*psiN*Sin(3._xp*(thAdj + theta3))
+ dtdtxAng = 486._xp*CN3**2*Cos(6._xp*(thAdj + theta3))*psiN + 24._xp*CN2*(1._xp + CN2*Cos(2._xp*(thAdj + theta2))) &
+ *(-(Cos(2._xp*(thAdj + theta2))*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))) + 2._xp*CN2*Sin(2._xp*(thAdj + theta2))**2)
+
+ yAng = 3._xp*Sqrt(3._xp)*CN3*Cos(3._xp*(thAdj + theta3))*Sqrt(psiN)*Sqrt(2._xp*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**3 + xAng)
+ dryAng = (yAng*(drpsiN/psiN + drxAng/(2._xp*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**3 + xAng)))/2._xp
+ dtyAng = yAng*(-3._xp*Tan(3._xp*(thAdj + theta3)) &
+ + (-12._xp*CN2*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**2*Sin(2._xp*(thAdj + theta2)) + dtxAng) &
+ /(2._xp*(2*(1._xp + CN2*Cos(2*(thAdj + theta2)))**3 + xAng)))
+ dtdtyAng = dtyAng**2/yAng + yAng*(-9._xp/Cos(3._xp*(thAdj + theta3))**2 &
+ - (-12._xp*CN2*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**2*Sin(2._xp*(thAdj + theta2)) + dtxAng)**2 &
+ /(2._xp*(2._xp*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**3 + xAng)**2) &
+ + (-24._xp*CN2*Cos(2._xp*(thAdj + theta2))*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**2 &
+ + 48._xp*CN2**2*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))*Sin(2._xp*(thAdj + theta2))**2 + dtdtxAng) &
+ /(2._xp*(2._xp*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**3 + xAng)))
+
+ rAng = atan2(yAng,xAng)/3.
+ drrAng = (-(yAng*drxAng) + xAng*dryAng)/(3._xp*(xAng**2 + yAng**2))
+ dtrAng = (-(yAng*dtxAng) + xAng*dtyAng)/(3._xp*(xAng**2 + yAng**2))
+ dtdtrAng = (-6._xp*dtrAng**2*yAng)/xAng &
+ + ((2._xp*yAng*dtxAng**2)/xAng - 2._xp*dtxAng*dtyAng - yAng*dtdtxAng + xAng*dtdtyAng)/(3._xp*(xAng**2 + yAng**2))
+
+ rShape = (aSurf*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))/Cos(3._xp*(thAdj + theta3))*&
+ &(-1._xp + Cos(rAng) + Sqrt(3._xp)*Sin(rAng)))/(3._xp*CN3)
+ drrShape = (aSurf*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))/Cos(3._xp*(thAdj + theta3))*&
+ &(Sqrt(3._xp)*Cos(rAng) - Sin(rAng))*drrAng)/(3._xp*CN3)
+ dtrShape = rShape*((-2._xp*CN2*Sin(2._xp*(thAdj + theta2)))/(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))&
+ &+ 3._xp*Tan(3._xp*(thAdj + theta3)) &
+ &+ ((Sqrt(3._xp)*Cos(rAng) - Sin(rAng))*dtrAng)/(-1._xp + Cos(rAng) + Sqrt(3._xp)*Sin(rAng)))
+ dtdtrShape = dtrShape**2/rShape + rShape*(9._xp - (4._xp*CN2*Cos(2._xp*(thAdj + theta2)))/(1._xp + CN2*Cos(2._xp*(thAdj + theta2))) &
+ - (4._xp*CN2**2*Sin (2._xp*(thAdj + theta2))**2)/(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**2 + 9._xp*Tan (3._xp*(thAdj + theta3))**2 &
+ + ((-4._xp + Cos(rAng) + Sqrt(3._xp)*Sin(rAng))*dtrAng**2)/(-1._xp + Cos(rAng) + Sqrt(3._xp)*Sin(rAng))**2 &
+ + ((Sqrt(3._xp)*Cos(rAng) - Sin(rAng))*dtdtrAng)/(-1._xp + Cos(rAng) + Sqrt(3._xp)*Sin(rAng)))
+
+ do i=1._xp,np
+ if (abs(CN3*cos(3._xp*(thAdj(i)+theta3)))<1e-8) then
+ rShape(i) = aSurf*Sqrt(psiN/(1._xp + CN2*Cos(2*(thAdj(i) + theta2))))
+ drrShape(i) = (rShape(i)*drpsiN)/(2._xp*psiN)
+ dtrShape(i) = (CN2*rShape(i)*Sin(2._xp*(thAdj(i) + theta2)))/(1._xp + CN2*Cos(2._xp*(thAdj(i) + theta2)))
+ dtdtrShape(i) = dtrShape(i)**2/rShape(i) &
+ + (2*(CN2**2 + CN2*Cos(2._xp*(thAdj(i) + theta2)))*rShape(i))/(1._xp + CN2*Cos(2._xp*(thAdj(i) + theta2)))**2
+ endif
+
+ if (abs(1._xp + CN2*Cos(2._xp*(thAdj(i) + theta2)))<1e-8) then
+ rShape(i) = aSurf*Sqrt(psiN)
+ drrShape(i) = 1._xp
+ dtrShape(i) = 0._xp
+ dtdtrShape(i) = 0._xp
+ endif
+ enddo
+
+ Rcenter = Raxis - (-R0 + Raxis)*rho_a**2
+ Zcenter = Zaxis - rho_a**2*(-Z0 + Zaxis)
+ drRcenter = -2._xp*(-R0 + Raxis)*rho_a
+ drZcenter = -2._xp*rho_a*(-Z0 + Zaxis)
+
+ R = Rcenter + Cos(thAdj)*rShape
+ Z = rShape*Sin(thAdj) + Zcenter
+ R_rho = (drRcenter + Cos(thAdj)*drrShape)/aSurf ! Adjust for rho deriv, not rho_a deriv
+ Z_rho = (drZcenter + Sin(thAdj)*drrShape)/aSurf ! Adjust for rho deriv, not rho_a deriv
+ R_theta = -(rShape*Sin(thAdj)) + Cos(thAdj)*dtrShape
+ Z_theta = Cos(thAdj)*rShape + Sin(thAdj)*dtrShape
+ R_theta_theta = -(Cos(thAdj)*rShape) - 2._xp*Sin(thAdj)*dtrShape + Cos(thAdj)*dtdtrShape
+ Z_theta_theta = -(rShape*Sin(thAdj)) + 2._xp*Cos(thAdj)*dtrShape + Sin(thAdj)*dtdtrShape
+
+ END SELECT
!dl/dtheta
dlp=(R_theta**2+Z_theta**2)**0.5_xp
diff --git a/testcases/CBC_ExBshear/fort_01.90 b/testcases/CBC_ExBshear/fort_01.90
index bf5ea9976e90d08811b3feb582aa81292f0f3262..45b8353393cd502bd52ad6871475bbc040c62096 100644
--- a/testcases/CBC_ExBshear/fort_01.90
+++ b/testcases/CBC_ExBshear/fort_01.90
@@ -57,10 +57,10 @@
HYP_V = 'hypcoll'
mu_p = 0.0
mu_j = 0.0
- nu = 1.0
+ nu = 1.00
beta = 0.0
ADIAB_E = .t.
- ExBrate = 0.01
+ ExBrate = 0.0
/
&CLOSURE
hierarchy_closure='truncation'
diff --git a/testcases/CBC_ExBshear/std/fort_00.90 b/testcases/CBC_ExBshear/std/fort_00.90
new file mode 100644
index 0000000000000000000000000000000000000000..b26e8952adf9141cf362ec36b9b315131b390ebe
--- /dev/null
+++ b/testcases/CBC_ExBshear/std/fort_00.90
@@ -0,0 +1,97 @@
+&BASIC
+ nrun = 99999999
+ dt = 0.01
+ tmax = 50
+ maxruntime = 72000
+ job2load = -1
+/
+&GRID
+ pmax = 2
+ jmax = 1
+ Nx = 64
+ Lx = 120
+ Ny = 48
+ Ly = 120
+ Nz = 16
+ SG = .f.
+ Nexc = 0
+/
+&GEOMETRY
+ geom = 's-alpha'
+ !geom = 'miller'
+ q0 = 1.4
+ shear = 0.8
+ eps = 0.18
+ kappa = 1.0
+ s_kappa= 0.0
+ delta = 0.0
+ s_delta= 0.0
+ zeta = 0.0
+ s_zeta = 0.0
+ parallel_bc = 'dirichlet'
+ shift_y= 0.0
+/
+&DIAGNOSTICS
+ dtsave_0d = 1
+ dtsave_1d = -1
+ dtsave_2d = -1
+ dtsave_3d = 1
+ dtsave_5d = 20
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .t.
+ write_Napj = .t.
+ write_dens = .t.
+ write_fvel = .t.
+ write_temp = .t.
+/
+&MODEL
+ LINEARITY = 'nonlinear'
+ Na = 1 ! number of species
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 1.0
+ HYP_V = 'hypcoll'
+ mu_p = 0.0
+ mu_j = 0.0
+ nu = 1.00
+ beta = 0.0
+ ADIAB_E = .t.
+ ExBrate = 0.0
+/
+&CLOSURE
+ hierarchy_closure='truncation'
+ dmax = -1
+ nonlinear_closure='truncation'
+ nmax = 0
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_ = 1.0
+ sigma_= 1.0
+ q_ = 1.0
+ k_N_ = 2.22
+ k_T_ = 6.96
+/
+
+&COLLISION
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ GK_CO = .f.
+ INTERSPECIES = .true.
+ mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL
+ INIT_OPT = 'blob'
+ ACT_ON_MODES = 'donothing'
+ init_background = 0.0
+ init_noiselvl = 0.005
+ iseed = 42
+/
+&TIME_INTEGRATION
+ !numerical_scheme = 'RK4'
+ numerical_scheme = 'SSP_RK3'
+/
diff --git a/testcases/CBC_ExBshear/std/fort_01.90 b/testcases/CBC_ExBshear/std/fort_01.90
new file mode 100644
index 0000000000000000000000000000000000000000..bf5ea9976e90d08811b3feb582aa81292f0f3262
--- /dev/null
+++ b/testcases/CBC_ExBshear/std/fort_01.90
@@ -0,0 +1,97 @@
+&BASIC
+ nrun = 99999999
+ dt = 0.01
+ tmax = 100
+ maxruntime = 72000
+ job2load = 0
+/
+&GRID
+ pmax = 2
+ jmax = 1
+ Nx = 64
+ Lx = 120
+ Ny = 48
+ Ly = 120
+ Nz = 16
+ SG = .f.
+ Nexc = 0
+/
+&GEOMETRY
+ geom = 's-alpha'
+ !geom = 'miller'
+ q0 = 1.4
+ shear = 0.8
+ eps = 0.18
+ kappa = 1.0
+ s_kappa= 0.0
+ delta = 0.0
+ s_delta= 0.0
+ zeta = 0.0
+ s_zeta = 0.0
+ parallel_bc = 'dirichlet'
+ shift_y= 0.0
+/
+&DIAGNOSTICS
+ dtsave_0d = 1
+ dtsave_1d = -1
+ dtsave_2d = -1
+ dtsave_3d = 1
+ dtsave_5d = 20
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .t.
+ write_Napj = .t.
+ write_dens = .t.
+ write_fvel = .t.
+ write_temp = .t.
+/
+&MODEL
+ LINEARITY = 'nonlinear'
+ Na = 1 ! number of species
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 1.0
+ HYP_V = 'hypcoll'
+ mu_p = 0.0
+ mu_j = 0.0
+ nu = 1.0
+ beta = 0.0
+ ADIAB_E = .t.
+ ExBrate = 0.01
+/
+&CLOSURE
+ hierarchy_closure='truncation'
+ dmax = -1
+ nonlinear_closure='truncation'
+ nmax = 0
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_ = 1.0
+ sigma_= 1.0
+ q_ = 1.0
+ k_N_ = 2.22
+ k_T_ = 6.96
+/
+
+&COLLISION
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ GK_CO = .f.
+ INTERSPECIES = .true.
+ mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL
+ INIT_OPT = 'blob'
+ ACT_ON_MODES = 'donothing'
+ init_background = 0.0
+ init_noiselvl = 0.005
+ iseed = 42
+/
+&TIME_INTEGRATION
+ !numerical_scheme = 'RK4'
+ numerical_scheme = 'SSP_RK3'
+/
diff --git a/testcases/CBC_ExBshear/fort_02.90 b/testcases/CBC_ExBshear/std/fort_02.90
similarity index 100%
rename from testcases/CBC_ExBshear/fort_02.90
rename to testcases/CBC_ExBshear/std/fort_02.90
diff --git a/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_fit_to_LLNL_GK.txt b/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_fit_to_LLNL_GK.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8e72be2dba002b521dd0788b41f3de7457e114aa
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_fit_to_LLNL_GK.txt
@@ -0,0 +1,18 @@
+5.988372093023255, -0.020876826722338038
+6.162790697674419, 0.4384133611691041
+6.308139534883722, 0.8559498956158684
+6.482558139534884, 1.1482254697286027
+6.627906976744185, 1.4822546972860131
+6.7441860465116275, 1.7954070981210872
+6.947674418604652, 2.0668058455114835
+7.093023255813954, 2.3382045929018798
+7.151162790697676, 2.44258872651357
+7.267441860465118, 2.713987473903968
+7.38372093023256, 2.8810020876826723
+7.500000000000002, 3.0480167014613784
+7.587209302325581, 3.1941544885177464
+7.674418604651164, 3.3402922755741127
+7.790697674418606, 3.549060542797495
+7.936046511627907, 3.757828810020877
+8.13953488372093, 4.050104384133613
+8.313953488372093, 4.237995824634655
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index c3418eb3a76d4cd59388c70b0438aebb07546daf..09b285b87e7526eed4415cc0159c3430f000b7ea 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -16,7 +16,7 @@ PARTITION = '';
% resdir = 'LM_DIIID_rho95/5x3x512x92x32';
% resdir = 'LM_DIIID_rho95/3x2x512x92x32';
% resdir = '../testcases/cyclone_example';
-% resdir = '../testcases/CBC_ExBshear';
+resdir = '../testcases/CBC_ExBshear/std';
% resdir = '../results/paper_3/HM_DTT_rho98/3x2x128x64x64';
%%
J0 = 00; J1 = 10;
@@ -49,7 +49,7 @@ if 0
data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
options.INTERP = 1;
options.POLARPLOT = 0;
-options.BWR = 1; % bluewhitered plot or gray
+options.BWR = 0; % bluewhitered plot or gray
options.CLIMAUTO = 1; % adjust the colormap auto
% options.NAME = '\phi';
% options.NAME = '\phi^{NZ}';
@@ -60,7 +60,7 @@ options.NAME = 'N_i^{00}';
% options.NAME = 'Q_x';
% options.NAME = 'n_i';
% options.NAME = 'n_i-n_e';
-options.PLAN = 'kxky';
+options.PLAN = 'xy';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
@@ -85,11 +85,12 @@ options.POLARPLOT = 0;
options.AXISEQUAL = 0;
options.NORMALIZE = 0;
options.LOGSCALE = 1;
-% options.NAME = 'N_i^{00}';
-options.NAME = '\phi';
+options.CLIMAUTO = 1;
+options.NAME = 'N_i^{00}';
+% options.NAME = '\phi';
options.PLAN = 'kxky';
options.COMP = 'avg';
-options.TIME = [0 10 50 100 500];
+options.TIME = [0 100];
% options.TIME = data.Ts3D(1:2:end);
options.RESOLUTION = 256;
fig = photomaton(data,options);
diff --git a/wk/lin_run_script.m b/wk/lin_run_script.m
index eddd7273611d1ccf7f977f4915f18d6420c6d9fd..13da10b74bff69506a55c95ff094a89455ed646b 100644
--- a/wk/lin_run_script.m
+++ b/wk/lin_run_script.m
@@ -25,23 +25,26 @@ EXECNAME = 'gyacomo23_dp'; % double precision
%% Setup parameters
% run lin_DTT_HM_rho85
% run lin_DTT_HM_rho98
-run lin_DIIID_LM_rho90
+% run lin_DIIID_LM_rho90
% run lin_DIIID_LM_rho95
% run lin_JET_rho97
% run lin_Entropy
-% run lin_ITG
+run lin_ITG
% run lin_KBM
%% Change parameters
-EXBRATE = 0.0; % Background ExB shear flow
-NY = 2;
-NX = 4;
-% PMAX = 4;
-% JMAX = PMAX/2;
-ky = 0.8; LY = 2*pi/ky;
-DT = 1e-4;
-TAU = 2.1;
+% EXBRATE = 0.0; % Background ExB shear flow
+NU = 0.1;
+CO = 'SG';
+GKCO = 1;
+NY = 16;
+NX = 2;
+PMAX = 4;
+JMAX = PMAX/2;
+ky = 0.1; LY = 2*pi/ky;
+% DT = 1e-4;
+% TAU = 2.1;
% SIGMA_E = 0.02;
-% TMAX = 50;
+TMAX = 25;
% DTSAVE0D = 200*DT;
% DTSAVE3D = TMAX/50;
%%-------------------------------------------------------------------------
@@ -52,7 +55,8 @@ setup
if RUN
MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
% RUN =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
- RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
+ % RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
+ RUN =['time ',mpirun,' -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
% RUN =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;'];
% RUN =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
% RUN = ['./../../../bin/gyacomo23_sp 0;'];
diff --git a/wk/parameters/lin_DIIID_LM_rho90.m b/wk/parameters/lin_DIIID_LM_rho90.m
index c8ff82e27f3ceebee4d0ddb69aea3da9f0e873ad..9995af00da479279d9fd14fab07431aa0c7780e0 100644
--- a/wk/parameters/lin_DIIID_LM_rho90.m
+++ b/wk/parameters/lin_DIIID_LM_rho90.m
@@ -20,8 +20,8 @@ P = 2;
J = P/2;%P/2;
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
-NX = 2; % real space x-gridpoints
-NY = 4; % real space y-gridpoints
+NX = 4; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
LY = 2*pi/0.3; % Size of the squared frequency domain in y direction
NZ = 32; % number of perpendicular planes (parallel grid)
diff --git a/wk/parameters/lin_DIIID_LM_rho95.m b/wk/parameters/lin_DIIID_LM_rho95.m
index 77a21328f20e3fe8936554d2f2a8d04ebd6e77ba..2bd5ba9b2f383bfb8a885ecb30666f6b79194743 100644
--- a/wk/parameters/lin_DIIID_LM_rho95.m
+++ b/wk/parameters/lin_DIIID_LM_rho95.m
@@ -1,7 +1,7 @@
%% Reference values
% See Neiser et al. 2019 Gyrokinetic GENE simulations of DIII-D near-edge L-mode plasmas
%% Set simulation parameters
-SIMID = 'lin_DTT_HM_rho90'; % Name of the simulation
+SIMID = 'lin_DIIID_HM_rho90'; % Name of the simulation
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
NU = 1.0; %(0.00235 in GENE)
@@ -20,7 +20,7 @@ P = 2;
J = P/2;%P/2;
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
-NX = 16; % real space x-gridpoints
+NX = 6; % real space x-gridpoints
NY = 2; % real space y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
LY = 2*pi/0.3; % Size of the squared frequency domain in y direction
@@ -46,7 +46,7 @@ NPOL = 1; % Number of poloidal turns
PB_PHASE = 0;
%% TIME PARAMETERS
TMAX = 15; % Maximal time unit
-DT = 1e-3; % Time step
+DT = 1e-4; % Time step
DTSAVE0D = 0.5; % Sampling time for 0D arrays
DTSAVE2D = -1; % Sampling time for 2D arrays
DTSAVE3D = 0.5; % Sampling time for 3D arrays