Skip to content
Snippets Groups Projects
Commit 4c3a5b4c authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

Debug of the geometry routine. Works fine for Miller and s-a. Zpinch must be tested...

parent f66b0353
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/diagnose.F90
+ 3
− 3
View file @ 4c3a5b4c
......@@ -158,9 +158,9 @@ SUBROUTINE diagnose_full(kstep)
CALL putarrnd(fidres, "/data/metric/hatZ", hatZ(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/hatB", hatB(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/hatB_NL", hatB_NL(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradxB", gradxB(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradyB", gradyB(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradzB", gradzB(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/dBdx", dBdx(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/dBdy", dBdy(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/dBdz", dBdz(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/Jacobian", Jacobian(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/1, 1, 3/))
......
src/geometry_mod.F90
+ 39
− 42
View file @ 4c3a5b4c
......@@ -55,7 +55,7 @@ implicit none
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gxx, gxy, gxz, gyy, gyz, gzz ! dimensions: z, odd/even p
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dxdr, dxdZ, Rc, phic, Zc
! derivatives of magnetic field strength
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gradxB, gradyB, gradzB
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dBdx, dBdy, dBdz
! Relative magnetic field strength
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB, hatB_NL
! Relative strength of major radius
......@@ -111,8 +111,8 @@ CONTAINS
ELSE
SELECT CASE(geom)
CASE('s-alpha')
IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry'
call eval_salphaB_geometry
IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha geometry'
call eval_salpha_geometry
CASE('Z-pinch')
IF( my_id .eq. 0 ) WRITE(*,*) 'Z-pinch geometry'
call eval_zpinch_geometry
......@@ -123,7 +123,7 @@ CONTAINS
call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta)
call get_miller(eps,major_R,major_Z,q0,shear,alpha_MHD,edge_opt,&
C_y,C_xy,dpdx_pm_geom,gxx,gyy,gzz,gxy,gxz,gyz,&
gradxB,gradyB,hatB,jacobian,gradzB,hatR,hatZ,dxdR,dxdZ,&
dBdx,dBdy,hatB,jacobian,dBdz,hatR,hatZ,dxdR,dxdZ,&
Ckxky,gradz_coeff)
CASE DEFAULT
STOP 'geometry not recognized!!'
......@@ -139,9 +139,10 @@ CONTAINS
ky = kyarray(iky)
DO ikx = ikxs, ikxe
kx = kxarray(ikx)
kparray(iky, ikx, iz, eo) = &
SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/hatB(iz,eo)
! there is a factor 1/B from the normalization; important to match GENE
! this factor comes from $b_a$ argument in the Bessel. Kperp is not used otherwise.
kparray(iky, ikx, iz, eo) = &
SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/ hatB(iz,eo)
ENDDO
ENDDO
ENDDO
......@@ -150,18 +151,24 @@ CONTAINS
G1 = gxy(iz,eo)*gxy(iz,eo)-gxx(iz,eo)*gyy(iz,eo)
G2 = gxy(iz,eo)*gxz(iz,eo)-gxx(iz,eo)*gyz(iz,eo)
G3 = gyy(iz,eo)*gxz(iz,eo)-gxy(iz,eo)*gyz(iz,eo)
Cx = (G1*gradyB(iz,eo) + G2*gradzB(iz,eo))/hatB(iz,eo)
Cy = (G3*gradzB(iz,eo) - G1*gradxB(iz,eo))/hatB(iz,eo)
! Here we divide by hatB because our equation is formulated with grad(lnB) terms (not gradB like in GENE)
Cx =-(dBdy(iz,eo) + G2/G1*dBdz(iz,eo))/hatB(iz,eo)
Cy = (dBdx(iz,eo) - G3/G1*dBdz(iz,eo))/hatB(iz,eo)
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
Ckxky(iky, ikx, iz,eo) = (Cx*kx + Cy*ky)
Ckxky(iky, ikx, iz,eo) = (Cx*kx + Cy*ky)/C_xy
ENDDO
ENDDO
! coefficient in the front of parallel derivative
gradz_coeff(iz,eo) = 1._dp / jacobian(iz,eo) / hatB(iz,eo)
gradz_coeff(iz,eo) = 1._dp /(jacobian(iz,eo)*hatB(iz,eo))
! Nonlinear term prefactor
! (according to my derivations, there should be a metric dependent factor in front of the Poisson bracket)
hatB_NL(iz,eo) = 1._dp !Jacobian(iz,eo)*(gxx(iz,eo)*gyy(iz,eo) - gxy(iz,eo)**2)/hatB(iz,eo)
ENDDO
ENDDO
......@@ -180,7 +187,7 @@ CONTAINS
!--------------------------------------------------------------------------------
!
SUBROUTINE eval_salphaB_geometry
SUBROUTINE eval_salpha_geometry
! evaluate s-alpha geometry model
implicit none
REAL(dp) :: z, kx, ky, Gx, Gy
......@@ -196,11 +203,11 @@ CONTAINS
gxz(iz,eo) = 0._dp
gyy(iz,eo) = 1._dp + (shear*z - alpha_MHD*SIN(z))**2
gyz(iz,eo) = 1._dp/eps
gzz(iz,eo) = 0._dp
gzz(iz,eo) = 1._dp/eps**2
dxdR(iz,eo)= COS(z)
dxdZ(iz,eo)= SIN(z)
! Relative strengh of radius
! Poloidal plane coordinates
hatR(iz,eo) = 1._dp + eps*COS(z)
hatZ(iz,eo) = 1._dp + eps*SIN(z)
......@@ -209,37 +216,24 @@ CONTAINS
phic(iz,eo) = z
Zc (iz,eo) = hatZ(iz,eo)
! Jacobian
Jacobian(iz,eo) = q0*hatR(iz,eo)
! Relative strengh of modulus of B
hatB(iz,eo) = 1._dp / hatR(iz,eo)
hatB(iz,eo) = 1._dp/(1._dp + eps*COS(z))
! Jacobian
Jacobian(iz,eo) = q0/hatB(iz,eo)
! Derivative of the magnetic field strenght
gradxB(iz,eo) = -COS(z) ! Gene put a factor hatB^2 or 1/hatR^2 in this
gradyB(iz,eo) = 0._dp
gradzB(iz,eo) = eps*SIN(z)/hatR(iz,eo) ! Gene put a factor hatB or 1/hatR in this
dBdx(iz,eo) = -COS(z)*hatB(iz,eo) ! LB = 1
dBdy(iz,eo) = 0._dp
dBdz(iz,eo) = eps*SIN(z)*hatB(iz,eo)**2
! Curvature operator
Gx = (gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo)) *eps*SIN(Z) ! Kx
Gy = -COS(z) + (gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)) *eps*SIN(Z) ! Ky
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - COS(z)*ky -(shear*z-alpha_MHD*SIN(z))*SIN(z)*ky)/ hatB(iz,eo)
! Ckxky(iky, ikx, iz,eo) = (Gx*kx + Gy*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
! Nonlinear term prefactor
hatB_NL(iz,eo) = 1._dp !Jacobian(iz,eo)*(gxx(iz,eo)*gyy(iz,eo) - gxy(iz,eo)**2)/hatB(iz,eo)
! Curvature factor
C_xy = 1._dp
ENDDO zloop
ENDDO parity
END SUBROUTINE eval_salphaB_geometry
END SUBROUTINE eval_salpha_geometry
!
!--------------------------------------------------------------------------------
!
......@@ -280,9 +274,9 @@ CONTAINS
hatB_NL(iz,eo) = 1._dp
! Derivative of the magnetic field strenght
gradxB(iz,eo) = -1._dp ! Gene put a factor hatB^2 or 1/hatR^2 in this
gradyB(iz,eo) = 0._dp
gradzB(iz,eo) = 0._dp ! Gene put a factor hatB or 1/hatR in this
dBdx(iz,eo) = -hatB(iz,eo) ! LB = 1
dBdy(iz,eo) = 0._dp
dBdz(iz,eo) = 0._dp ! Gene put a factor hatB or 1/hatR in this
! Curvature operator
DO iky = ikys, ikye
......@@ -431,10 +425,13 @@ CONTAINS
! DO iky = ikys,ikye
! print*, ikx_zBC_L(iky,:)
! enddo
! print*, pb_phase_L
! write(*,*) 'ikx_zBC_R :-----------'
! DO iky = ikys,ikye
! print*, ikx_zBC_R(iky,:)
! enddo
! print*, pb_phase_R
! print*, shift_y
! stop
!!!!!!! Example of maps ('x' means connected to 0 value, in the table it is -99)
! dirichlet connection map BC of the RIGHT boundary (z=pi*Npol-dz)
......@@ -481,9 +478,9 @@ END SUBROUTINE set_ikx_zBC_map
CALL allocate_array( gyy,izgs,izge, 0,1)
CALL allocate_array( gyz,izgs,izge, 0,1)
CALL allocate_array( gzz,izgs,izge, 0,1)
CALL allocate_array( gradxB,izgs,izge, 0,1)
CALL allocate_array( gradyB,izgs,izge, 0,1)
CALL allocate_array( gradzB,izgs,izge, 0,1)
CALL allocate_array( dBdx,izgs,izge, 0,1)
CALL allocate_array( dBdy,izgs,izge, 0,1)
CALL allocate_array( dBdz,izgs,izge, 0,1)
CALL allocate_array( hatB,izgs,izge, 0,1)
CALL allocate_array( hatB_NL,izgs,izge, 0,1)
CALL allocate_array( hatR,izgs,izge, 0,1)
......
src/moments_eq_rhs_mod.F90
+ 6
− 6
View file @ 4c3a5b4c
......@@ -33,7 +33,7 @@ SUBROUTINE moments_eq_rhs_e
USE model
USE prec_const
USE collision
USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatB
USE geometry, ONLY: gradz_coeff, dBdz, Ckxky, hatB_NL, hatB
USE calculus, ONLY : interp_z, grad_z, grad_z2
IMPLICIT NONE
......@@ -121,13 +121,13 @@ SUBROUTINE moments_eq_rhs_e
!! Sum of all RHS terms
moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
-imagu*Ckxky(iky,ikx,iz,eo)*hatB(iz,eo) * Tperp&
-imagu*Ckxky(iky,ikx,iz,eo) * Tperp&
! Perpendicular pressure effects
-i_ky*beta*dpdx * Tperp&
! Parallel coupling (Landau Damping)
-gradz_coeff(iz,eo) * Tpar &
! Mirror term (parallel magnetic gradient)
-gradzB(iz,eo)*gradz_coeff(iz,eo) * Tmir&
-dBdz(iz,eo)*gradz_coeff(iz,eo) * Tmir&
! Drives (density + temperature gradients)
-i_ky * (Tphi - Tpsi) &
! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose)
......@@ -182,7 +182,7 @@ SUBROUTINE moments_eq_rhs_i
USE model
USE prec_const
USE collision
USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatB
USE geometry, ONLY: gradz_coeff, dBdz, Ckxky, hatB_NL, hatB
USE calculus, ONLY : interp_z, grad_z, grad_z2
IMPLICIT NONE
......@@ -272,13 +272,13 @@ SUBROUTINE moments_eq_rhs_i
!! Sum of all RHS terms
moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
-imagu*Ckxky(iky,ikx,iz,eo)*hatB(iz,eo) * Tperp &
-imagu*Ckxky(iky,ikx,iz,eo) * Tperp &
! Perpendicular pressure effects
-i_ky*beta*dpdx * Tperp&
! Parallel coupling (Landau damping)
-gradz_coeff(iz,eo) * Tpar &
! Mirror term (parallel magnetic gradient)
-gradzB(iz,eo)*gradz_coeff(iz,eo) * Tmir &
-dBdz(iz,eo)*gradz_coeff(iz,eo) * Tmir &
! Drives (density + temperature gradients)
-i_ky * (Tphi - Tpsi) &
! Numerical hyperdiffusion (totally artificial, for stability purpose)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment