From 286b37a44f5e4ac56e3a33b7565b29dda5066796 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Wed, 27 Oct 2021 15:29:28 +0200
Subject: [PATCH] Update scripts
---
matlab/post_processing.m | 4 +-
src/array_mod.F90 | 5 +
src/collision_mod.F90 | 8 +-
src/compute_Sapj.F90 | 8 +-
src/geometry_mod.F90 | 175 ++++++++++++++-------------------
src/inital.F90 | 4 +
src/memory.F90 | 4 +
src/processing_mod.F90 | 42 +++++++-
src/srcinfo.h | 4 +-
src/srcinfo/srcinfo.h | 4 +-
src/stepon.F90 | 11 ++-
wk/analysis_3D.m | 25 ++---
wk/linear_1D_entropy_mode.m | 20 ++--
wk/molix_helaz_benchmark.m | 21 ++++
wk/molix_plot_phi.m | 127 ++++++++++++++++++++++++
wk/plot_cosol_mat.m | 2 +-
wk/plot_phi_ballooning.m | 21 ++--
wk/shearless_linear_fluxtube.m | 12 +--
18 files changed, 335 insertions(+), 162 deletions(-)
create mode 100644 wk/molix_helaz_benchmark.m
create mode 100644 wk/molix_plot_phi.m
diff --git a/matlab/post_processing.m b/matlab/post_processing.m
index 14e8b673..434be93a 100644
--- a/matlab/post_processing.m
+++ b/matlab/post_processing.m
@@ -165,7 +165,7 @@ g_I = zeros(Nkx,Nky,Nz);
for ikx = 1:Nkx
for iky = 1:Nky
for iz = 1:Nz
- [g_I(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin),squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
+ [g_I(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin)',squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
end
end
end
@@ -186,7 +186,7 @@ g_II = zeros(Nkx,Nky);
for ikx = 1:Nkx
for iky = 1
for iz = 1:Nz
- [g_II(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin),squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
+ [g_II(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin)',squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
end
end
end
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index ccdb267a..92dd2d63 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -6,6 +6,11 @@ MODULE array
! Arrays to store the rhs, for time integration (ip,ij,ikx,iky,iz,updatetlevel)
COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: moments_rhs_e
COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: moments_rhs_i
+
+ ! Arrays of non-adiabatique moments
+ COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: nadiab_moments_e
+ COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: nadiab_moments_i
+
! Non linear term array (ip,ij,ikx,iky,iz)
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Sepj ! electron
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Sipj ! ion
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index 960d7356..eaa0b706 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -548,7 +548,7 @@ CONTAINS
DO ikp = ikps_C,ikpe_C ! Loop over everz kperp values
! we put zeros if kp>2/3kpmax because thoses frequencies are filtered through AA
- IF(kp_grid_mat(ikp) .GT. two_third_kpmax .AND. NON_LIN) THEN
+ IF( (kp_grid_mat(ikp) .GT. two_third_kpmax) .AND. NON_LIN) THEN
CiepjT_kp(:,:,ikp) = 0._dp
CiepjF_kp(:,:,ikp) = 0._dp
CeipjT_kp(:,:,ikp) = 0._dp
@@ -696,6 +696,8 @@ CONTAINS
ikp_next = ikp_mat !index of k1 (larger than kperp_sim thanks to previous loop)
ikp_prev = ikp_mat - 1 !index of k0 (smaller neighbour to interpolate inbetween)
! write(*,*) kp_grid_mat(ikp_prev), '<', kperp_sim, '<', kp_grid_mat(ikp_next)
+ if ( (kp_grid_mat(ikp_prev) .GT. kperp_sim) .OR. (kp_grid_mat(ikp_next) .LT. kperp_sim) )&
+ write(*,*) 'Warning, linear interp of collision matrix failed!!'
! 0->1 variable for linear interp, i.e. zero2one = (k-k0)/(k1-k0)
zerotoone = (kperp_sim - kp_grid_mat(ikp_prev))/(kp_grid_mat(ikp_next) - kp_grid_mat(ikp_prev))
@@ -718,8 +720,8 @@ CONTAINS
CiepjF(:,:,1,1) = CiepjF_kp(:,:,1)
ENDIF
! Deallocate auxiliary variables
- DEALLOCATE (Ceepj__kp ); DEALLOCATE (CeipjT_kp); DEALLOCATE (CeipjF_kp)
- DEALLOCATE (Ciipj__kp ); DEALLOCATE (CiepjT_kp); DEALLOCATE (CiepjF_kp)
+ DEALLOCATE (Ceepj__kp); DEALLOCATE (CeipjT_kp); DEALLOCATE (CeipjF_kp)
+ DEALLOCATE (Ciipj__kp); DEALLOCATE (CiepjT_kp); DEALLOCATE (CiepjF_kp)
IF( CO_AA_ONLY ) THEN
CeipjF = 0._dp;
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index ca6120cd..967e1039 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -43,9 +43,7 @@ zloop: DO iz = izs,ize
real_data_c = 0._dp ! initialize sum over real nonlinear term
! Set non linear sum truncation
- IF (NL_CLOS .EQ. -3) THEN
- return
- ELSEIF (NL_CLOS .EQ. -2) THEN
+ IF (NL_CLOS .EQ. -2) THEN
nmax = Jmaxe
ELSEIF (NL_CLOS .EQ. -1) THEN
nmax = Jmaxe-(ij-1)
@@ -141,9 +139,7 @@ zloop: DO iz = izs,ize
real_data_c = 0._dp ! initialize sum over real nonlinear term
! Set non linear sum truncation
- IF (NL_CLOS .EQ. -3) THEN
- return
- ELSEIF (NL_CLOS .EQ. -2) THEN
+ IF (NL_CLOS .EQ. -2) THEN
nmax = Jmaxi
ELSEIF (NL_CLOS .EQ. -1) THEN
nmax = Jmaxi-(ij-1)
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 72e13ac7..165bab1c 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -18,127 +18,98 @@ contains
! evalute metrix, elements, jacobian and gradient
implicit none
!
- IF( my_id .eq. 0 ) WRITE(*,*) 'circular geometry'
- ! circular model
- call eval_circular_geometry
+ IF( (Ny .EQ. 1) .AND. (Nz .EQ. 1)) THEN !1D perp linear run
+ IF( my_id .eq. 0 ) WRITE(*,*) '1D perpendicular geometry'
+ call eval_1D_geometry
+ ELSE
+ IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry'
+ call eval_salphaB_geometry
+ ENDIF
!
END SUBROUTINE eval_magnetic_geometry
!
!--------------------------------------------------------------------------------
!
- subroutine eval_circular_geometry
- ! evaluate circular geometry model
- ! Ref: Lapilonne et al., PoP, 2009
- implicit none
- REAL(dp) :: z, kx, ky
+ subroutine eval_salphaB_geometry
+ ! evaluate s-alpha geometry model
+ implicit none
+ REAL(dp) :: z, kx, ky
- zloop: DO iz = izs,ize
- z = zarray(iz)
+ zloop: DO iz = izs,ize
+ z = zarray(iz)
- ! Metric elements
+ ! metric
gxx(iz) = 1._dp
- gxy(iz) = shear*z -eps*SIN(z)
- gyy(iz) = 1._dp +(shear*z)**2 -2._dp*eps*COS(z) -2._dp*shear*eps*z*SIN(z)
+ gxy(iz) = shear*z
+ gyy(iz) = 1._dp + (shear*z)**2
- ! Relative strengh of radius
- hatR(iz) = 1._dp + eps*cos(z)
+ ! Relative strengh of radius
+ hatR(iz) = 1._dp + eps*COS(z)
- hatB(iz) = SQRT(gxx(iz)*gyy(iz) - gxy(iz)**2)
+ ! Jacobian
+ Jacobian(iz) = q0*hatR(iz)
- ! Jacobian
- Jacobian(iz) = q0*hatR(iz)**2
+ ! Relative strengh of modulus of B
+ hatB(iz) = 1._dp / hatR(iz)
- ! Derivative of the magnetic field strenght
- gradxB(iz) = -(COS(z) + eps*SIN(z)**2)/hatB(iz)/hatB(iz)
- gradzB(iz) = eps*SIN(z)
- ! coefficient in the front of parallel derivative
- gradz_coeff(iz) = 1._dp / Jacobian(iz) / hatB(iz)
+ ! Derivative of the magnetic field strenght
+ gradxB(iz) = -COS(z)
+ gradzB(iz) = eps * SIN(z) / hatR(iz)
- ! Curvature operator
+ ! Curvature operator
DO iky = ikys, ikye
ky = kyarray(iky)
- DO ikx= ikxs, ikxe
- kx = kxarray(ikx)
- Ckxky(ikx,iky,iz) = -SIN(z)*kx -(COS(z)+shear*z*SIN(z))*ky
- ENDDO
+ DO ikx= ikxs, ikxe
+ kx = kxarray(ikx)
+ Ckxky(ikx, iky, iz) = (-SIN(z)*kx - (COS(z) + shear* z* SIN(z))*ky) * hatB(iz) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ ENDDO
ENDDO
- IF (Nky .EQ. 1) THEN ! linear 1D run we switch kx and ky for parallel opt
- DO iky = ikys, ikye
- ky = kyarray(iky)
+ ! coefficient in the front of parallel derivative
+ gradz_coeff(iz) = 1._dp / Jacobian(iz) / hatB(iz)
+ ENDDO zloop
+
+ END SUBROUTINE eval_salphaB_geometry
+ !
+ !--------------------------------------------------------------------------------
+ !
+
+ subroutine eval_1D_geometry
+ ! evaluate 1D perp geometry model
+ implicit none
+ REAL(dp) :: z, kx, ky
+
+ zloop: DO iz = izs,ize
+ z = zarray(iz)
+
+ ! metric
+ gxx(iz) = 1._dp
+ gxy(iz) = 0._dp
+ gyy(iz) = 1._dp
+
+ ! Relative strengh of radius
+ hatR(iz) = 1._dp
+
+ ! Jacobian
+ Jacobian(iz) = 1._dp
+
+ ! Relative strengh of modulus of B
+ hatB(iz) = 1._dp
+
+ ! Curvature operator
+ DO iky = ikys, ikye
+ ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
- Ckxky(ikx,iky,iz) = -SIN(z)*ky -(COS(z)+shear*z*SIN(z))*kx
+ Ckxky(ikx, iky, iz) = -kx ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
- ENDDO
- ENDIF
+ ENDDO
+ ! coefficient in the front of parallel derivative
+ gradz_coeff(iz) = 1._dp
ENDDO zloop
- END SUBROUTINE eval_circular_geometry
- !
- !--------------------------------------------------------------------------------
- ! UNUSED
- ! subroutine eval_salpha_geometry
- ! ! evaluate s-alpha geometry model
- ! implicit none
- ! REAL(dp) :: z, kx, ky
- !
- ! zloop: DO iz = izs,ize
- ! z = zarray(iz)
- ! gxx(iz) = 1._dp
- ! gxy(iz) = shear*z
- ! gyy(iz) = 1._dp + (shear*z)**2
- ! gyz(iz) = 1._dp/eps
- ! gxz(iz) = 0._dp
- !
- ! ! Relative strengh of radius
- ! hatR(iz) = 1._dp + eps*cos(z)
- !
- ! ! Jacobian
- ! Jacobian(iz) = q0*hatR(iz)
- !
- ! ! Relative strengh of modulus of B
- ! hatB(iz) = 1._dp / hatR( iz)
- !
- ! ! Derivative of the magnetic field strenght
- ! gradxB(iz) = - cos( z) / hatR(iz)
- ! gradzB(iz) = eps * sin(z)
- !
- ! ! Gemoetrical coefficients for the curvature operator
- ! ! Note: Gamma2 and Gamma3 are obtained directly form Gamma1 in the expression of the curvature operator implemented here
- ! !
- ! Gamma1(iz) = gxy(iz) * gxy(iz) - gxx(iz) * gyy(iz)
- ! Gamma2(iz) = gxz(iz) * gxy(iz) - gxx(iz) * gyz(iz)
- ! Gamma3(iz) = gxz(iz) * gyy(iz) - gxy(iz) * gyz(iz)
- !
- ! ! Curvature operator
- ! DO iky = ikys, ikye
- ! ky = kyarray(iky)
- ! DO ikx= ikxl, ikxr
- ! kx = kxarray(ikx,iky)
- ! Cxy(ikx, iky, iz) = (-sin(z)*kx - (cos(z) + shear* z* sin(z))*ky) * hatB(iz) ! .. 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) = 1._dp / Jacobian(iz) / hatB(iz)
- ! ENDDO
- !
- ! ! Evaluate perpendicular wavenumber
- ! ! k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy}
- ! ! normalized to rhos_
- ! DO iky = ikys, ikye
- ! ky = kyarray(iky)
- ! DO ikx = ikxl, ikxr
- ! kx = kxarray(ikx,iky)
- ! kperp_array(ikx, iky, iz) = sqrt( gxx(iz)*kx**2 + 2._dp* gxy(iz) * kx*ky + gyy(iz)* ky**2) !! / hatB( iz) ! there is a factor 1/B from the normalization; important to match GENE
- ! ENDDO
- ! ENDDO
- ! ENDDO zloop
- ! !
- ! IF( me .eq. 0 ) PRINT*, 'max kperp = ', maxval( kperp_array)
- !
- ! END SUBROUTINE eval_salpha_geometry
- !
- !--------------------------------------------------------------------------------
- !
+
+ END SUBROUTINE eval_1D_geometry
+ !
+ !--------------------------------------------------------------------------------
+ !
end module geometry
diff --git a/src/inital.F90 b/src/inital.F90
index 1f08a444..3bac2e89 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -14,6 +14,7 @@ SUBROUTINE inital
USE ghosts
USE restarts
USE numerics, ONLY: wipe_turbulence, wipe_zonalflow
+ USE processing, ONLY: compute_nadiab_moments
IMPLICIT NONE
CALL set_updatetlevel(1)
@@ -63,6 +64,9 @@ SUBROUTINE inital
IF (my_id .EQ. 0) WRITE(*,*) 'Ghosts communication'
CALL update_ghosts
+ IF (my_id .EQ. 0) WRITE(*,*) 'Computing non adiab moments'
+ CALL compute_nadiab_moments
+
!!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!!
IF ( NON_LIN ) THEN;
IF (my_id .EQ. 0) WRITE(*,*) 'Init Sapj'
diff --git a/src/memory.F90 b/src/memory.F90
index 31d098e0..2d1dee0a 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -40,6 +40,10 @@ SUBROUTINE memory
CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
+ ! Non linear terms and dnjs table
+ CALL allocate_array( nadiab_moments_e, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( nadiab_moments_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize)
+
! Collision term
CALL allocate_array( TColl_e, ips_e,ipe_e, ijs_e,ije_e , ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( TColl_i, ips_i,ipe_i, ijs_i,ije_i , ikxs,ikxe, ikys,ikye, izs,ize)
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 09f164bd..73583aa4 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -9,7 +9,7 @@ MODULE processing
REAL(dp), PUBLIC, PROTECTED :: pflux_ri, gflux_ri
REAL(dp), PUBLIC, PROTECTED :: hflux_x
- PUBLIC :: compute_density, compute_temperature
+ PUBLIC :: compute_density, compute_temperature, compute_nadiab_moments
PUBLIC :: compute_radial_ion_transport, compute_radial_heatflux
CONTAINS
@@ -126,6 +126,46 @@ SUBROUTINE compute_radial_heatflux
ENDIF
END SUBROUTINE compute_radial_heatflux
+SUBROUTINE compute_nadiab_moments
+ ! evaluate the non-adiabatique ion moments
+ !
+ ! n_{pi} = N^{pj} + kernel(j) /tau_i phi delta_p0
+ !
+ USE fields, ONLY : moments_i, moments_e, phi
+ USE array, ONLY : kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i
+ USE time_integration, ONLY : updatetlevel
+ USE model, ONLY : qe_taue, qi_taui
+ implicit none
+
+ ! Add non-adiabatique term
+ DO ip=ipsg_e,ipeg_e
+ IF(parray_e(ip) .EQ. 0) THEN
+ DO ij=ijsg_e,ijeg_e
+ nadiab_moments_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize)&
+ = moments_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) &
+ + qe_taue*kernel_e(ij,ikxs:ikxe,ikys:ikye,izs:ize)*phi(ikx,iky,iz)
+ ENDDO
+ ELSE
+ nadiab_moments_e(ip,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize) &
+ = moments_e(ip,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel)
+ ENDIF
+ ENDDO
+ ! Add non-adiabatique term
+ DO ip=ipsg_i,ipeg_i
+ IF(parray_i(ip) .EQ. 0) THEN
+ DO ij=ijsg_i,ijeg_i
+ nadiab_moments_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize)&
+ = moments_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) &
+ + qi_taui*kernel_i(ij,ikxs:ikxe,ikys:ikye,izs:ize)*phi(ikx,iky,iz)
+ ENDDO
+ ELSE
+ nadiab_moments_i(ip,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize) &
+ = moments_i(ip,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel)
+ ENDIF
+ ENDDO
+ !
+END SUBROUTINE compute_nadiab_moments
+
! Compute the 2D particle density for electron and ions (sum over Laguerre)
SUBROUTINE compute_density
USE fields, ONLY : moments_i, moments_e, phi
diff --git a/src/srcinfo.h b/src/srcinfo.h
index 04893679..746ccf4f 100644
--- a/src/srcinfo.h
+++ b/src/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='d118c2b-dirty')
+parameter (VERSION='f7c491b-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Thu Oct 7 17:18:16 CEST 2021')
+parameter (EXECDATE='Thu Oct 21 10:50:35 CEST 2021')
parameter (HOST ='spcpc606')
diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h
index 04893679..746ccf4f 100644
--- a/src/srcinfo/srcinfo.h
+++ b/src/srcinfo/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='d118c2b-dirty')
+parameter (VERSION='f7c491b-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Thu Oct 7 17:18:16 CEST 2021')
+parameter (EXECDATE='Thu Oct 21 10:50:35 CEST 2021')
parameter (HOST ='spcpc606')
diff --git a/src/stepon.F90 b/src/stepon.F90
index b7a36785..26b3df80 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -13,6 +13,7 @@ SUBROUTINE stepon
use prec_const
USE time_integration
USE numerics, ONLY: wipe_zonalflow, wipe_turbulence
+ USE processing, ONLY: compute_nadiab_moments
USE utility, ONLY: checkfield
IMPLICIT NONE
@@ -23,7 +24,7 @@ SUBROUTINE stepon
DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
!----- BEFORE: All fields are updated for step = n
! Compute right hand side from current fields
- ! N_rhs(N_n,phi_n, S_n, Tcoll_n)
+ ! N_rhs(N_n, nadia_n, phi_n, S_n, Tcoll_n)
CALL moments_eq_rhs_e
CALL moments_eq_rhs_i
! ---- step n -> n+1 transition
@@ -40,6 +41,8 @@ SUBROUTINE stepon
CALL compute_TColl
! Update electrostatic potential phi_n = phi(N_n+1)
CALL poisson
+ ! Update non adiabatic moments n -> n+1
+ CALL compute_nadiab_moments
! Update nonlinear term S_n -> S_n+1(phi_n+1,N_n+1)
IF ( NON_LIN ) CALL compute_Sapj
! Cancel zonal modes artificially
@@ -61,7 +64,7 @@ SUBROUTINE stepon
CALL cpu_time(t0_checkfield)
IF(NON_LIN) CALL anti_aliasing ! ensure 0 mode for 2/3 rule
- IF(NON_LIN) CALL enforce_symetry ! Enforcing symmetry on kx = 0
+ IF(NON_LIN) CALL enforce_symmetry ! Enforcing symmetry on kx = 0
mlend=.FALSE.
IF(.NOT.nlend) THEN
@@ -109,7 +112,7 @@ SUBROUTINE stepon
END DO
END SUBROUTINE anti_aliasing
- SUBROUTINE enforce_symetry ! Force X(k) = X(N-k)* complex conjugate symmetry
+ SUBROUTINE enforce_symmetry ! Force X(k) = X(N-k)* complex conjugate symmetry
IF ( contains_kx0 ) THEN
! Electron moments
DO ip=ips_e,ipe_e
@@ -142,6 +145,6 @@ SUBROUTINE stepon
! must be real at origin
phi(ikx_0,iky_0,:) = REAL(phi(ikx_0,iky_0,:))
ENDIF
- END SUBROUTINE enforce_symetry
+ END SUBROUTINE enforce_symmetry
END SUBROUTINE stepon
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 5f32c1b5..cf0111b1 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -12,16 +12,16 @@ outfile ='simulation_A/1024x256_3x2_L_120_kN_1.6667_nu_1e-01_DGGK';
CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
system(CMD);
else% Marconi results
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/500x500_L_120_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_0e+00/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/500x500_L_120_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_0e+00/out.txt';
BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
end
%% Load the results
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 02; JOBNUMMAX = 02;
+JOBNUMMIN = 03; JOBNUMMAX = 03;
% JOBNUMMIN = 07; JOBNUMMAX = 20; % For CO damping sim A
compile_results %Compile the results from first output found to JOBNUMMAX if existing
@@ -107,8 +107,8 @@ FIELD = squeeze(plt(FIELD));
GIFNAME = [NAME,sprintf('_%.2d',JOBNUM),'_',PARAMS];
% Create movie (gif or mp4)
-% create_gif
-create_mov
+create_gif
+% create_mov
end
if 0
@@ -130,7 +130,7 @@ FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}'
% FIELD = dens_e-Z_n_e-(Z_phi-phi); FNAME = 'Non_adiab_part'; FIELDLTX = 'n_e^{NZ}-\phi^{NZ}';
% Chose when to plot it
-tf = 200:200:1000;
+tf = 1500:500:3000;
% Sliced
ix = 1; iy = 1; iz = 1;
@@ -172,7 +172,7 @@ FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
% FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e';
% Chose when to plot it
-tf = 200:400:1400;
+tf = 1500:500:3000;
% tf = 8000;
% Sliced
@@ -212,18 +212,19 @@ plot_param_evol
end
if 0
-%% Radial shear profile
-tf = 2000+[900:20:1100];
+%% Radial shear profile (with moving average)
+tf = 1000+[0:100:1000];
ymax = 0;
figure
for i_ = 1:numel(tf)
[~,it] = min(abs(Ts3D-tf(i_)));
-data = squeeze((mean(dxphi(:,:,1,it),2)));
-plot(x,data,'Displayname',sprintf('$t c_s/R=%.0f$',Ts3D(it))); hold on;
+data = squeeze((mean(dx2phi(:,:,1,it),2)));
+step = 50;
+plot(movmean(x,step),movmean(data,step),'Displayname',sprintf('$t c_s/R=%.0f$',Ts3D(it))); hold on;
ymax = max([ymax abs(min(data)) abs(max(data))]);
end
xlim([min(x), max(x)]); ylim(1.2*[-1 1]*ymax);
-xlabel('$x/\rho_s$'); ylabel('$v_{E\times B,x}$'); grid on
+xlabel('$x/\rho_s$'); ylabel('$s_{E\times B,x}$'); grid on
end
if 1
diff --git a/wk/linear_1D_entropy_mode.m b/wk/linear_1D_entropy_mode.m
index 04081247..e1cd8d83 100644
--- a/wk/linear_1D_entropy_mode.m
+++ b/wk/linear_1D_entropy_mode.m
@@ -6,23 +6,23 @@ default_plots_options
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.01; % Collision frequency
+NU = 0.1; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_N = 1/0.6; % Density gradient drive
K_T = 0.0; % Temperature '''
K_E = 0.0; % Electrostat '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
%% GRID PARAMETERS
-Nx = 82; % real space x-gridpoints
+Nx = 100; % real space x-gridpoints
Ny = 1; % '' y-gridpoints
Lx = 150; % Size of the squared frequency domain
-Ly = 0; % Size of the squared frequency domain
+Ly = 1; % Size of the squared frequency domain
Nz = 1; % number of perpendicular planes (parallel grid)
q0 = 1.0; % safety factor
shear = 0.0; % magnetic shear
eps = 0.0; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 300; % Maximal time unit
+TMAX = 100; % Maximal time unit
DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
@@ -35,7 +35,7 @@ SIMID = 'test_4.1'; % Name of the simulation
NON_LIN = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
% (0:L.Bernstein, 1:Dougherty, 2:Sugama, 3:Pitch angle, 4:Full Couloumb ; +/- for GK/DK)
-CO = 3;
+CO = 2;
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
@@ -88,9 +88,9 @@ for i = 1:Nparam
system(['rm fort*.90']);
% Run linear simulation
if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3.82 1 6 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ./../../../bin/helaz 1 2 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
end
% Load and process results
%%
@@ -111,10 +111,10 @@ for i = 1:Nparam
end
gamma_Ni00(i,:) = real(gamma_Ni00(i,:));% .* (gamma_Ni00(i,:)>=0.0));
gamma_Nipj(i,:) = real(gamma_Nipj(i,:));% .* (gamma_Nipj(i,:)>=0.0));
- if 1
+ if 0
%% Fit verification
figure;
- for i = 1:4:Nx/2+1
+ for i = 1:1:Nx/2+1
X_ = Ts3D(:); Y_ = squeeze(abs(Ni00(i,1,1,:)));
semilogy(X_,Y_,'DisplayName',['k=',num2str(kx(i))]); hold on;
end
diff --git a/wk/molix_helaz_benchmark.m b/wk/molix_helaz_benchmark.m
new file mode 100644
index 00000000..0863988a
--- /dev/null
+++ b/wk/molix_helaz_benchmark.m
@@ -0,0 +1,21 @@
+%%
+
+cd ..
+system('make');
+cd wk
+
+shearless_linear_fluxtube
+
+
+cd ../../molix
+system('make');
+system('./bin/molix');
+cd ../HeLaZ/wk
+%%
+time_2_plot = 5.0;
+[Y_,X_]=molix_plot_phi('../../molix/field.dat.h5',time_2_plot);
+
+plot_phi_ballooning; hold on
+plot(X_/pi,real(Y_),'ob');
+plot(X_/pi,imag(Y_),'or');
+plot(X_/pi,abs(Y_) ,'ok');
\ No newline at end of file
diff --git a/wk/molix_plot_phi.m b/wk/molix_plot_phi.m
new file mode 100644
index 00000000..b94bd844
--- /dev/null
+++ b/wk/molix_plot_phi.m
@@ -0,0 +1,127 @@
+function [phib, b_angle] = molix_plot_phi(resfile,varargin)
+% plot ballooning representation at the specified iput by time2plot for
+% each ky modes
+
+
+set(0,'defaulttextInterpreter','latex')
+set(0, 'defaultAxesTickLabelInterpreter','latex');
+set(0, 'defaultLegendInterpreter','latex');
+set(0,'defaultAxesFontSize',16);
+set(0, 'DefaultLineLineWidth', 1.5);
+
+
+% read data and attributes
+coordkx = h5read(resfile,'/data/var2d/phi/coordkx');
+dimskx = size(coordkx);
+Nkx = (dimskx(1) - 1)/2;
+% Nkx = (length(coordkx)-1)/2
+coordz = h5read(resfile,'/data/var2d/phi/coordz');
+coordky = h5read(resfile,'/data/var2d/phi/coordky');
+Nky = length(coordky);
+Nz = length(coordz);
+coordtime =h5read(resfile,'/data/var2d/time');
+Nt = length(coordtime);
+pmax = double(h5readatt(resfile,'/data/input','pmaxi'));
+jmax = double(h5readatt(resfile,'/data/input','jmaxi'));
+epsilon = double(h5readatt(resfile,'/data/input','inv. aspect ratio'));
+safety_fac = double(h5readatt(resfile,'/data/input','safety factor'));
+nu = double(h5readatt(resfile,'/data/input','nu'));
+shear = double(h5readatt(resfile,'/data/input','magnetic shear'));
+
+RTi= double(h5readatt(resfile,'/data/input','RTi'));
+RTe= double(h5readatt(resfile,'/data/input','RTe'));
+Rn= double(h5readatt(resfile,'/data/input','Rni'));
+
+% read electrostatic pot.
+if(nargin == 1 )
+ % plot end of the simulation
+ if( Nt >1)
+ [~,idxte] = min(abs(coordtime -coordtime(end-1)));
+ else
+ idxte =0;
+ end
+elseif(nargin ==2)
+ time2plot = varargin{1};
+ [~,idxte] = min(abs(coordtime -time2plot));
+end
+
+iframe = idxte-1;
+
+dataset = ['/data/var2d/phi/',num2str(iframe,'%06d')];
+phi = h5read(resfile,dataset);
+
+% read eigenvalu
+% dataset = ['/data/var2d/eig/growth_rate'];
+% growth_rate = h5read(resfile,dataset);
+% dataset = ['/data/var2d/eig/freq'];
+% frequency = h5read(resfile,dataset);
+
+dataset = ['/data/var2d/time'];
+time2d= h5read(resfile,dataset);
+
+% Apply baollooning tranform
+for iky=1:Nky
+ dims = size(phi.real);
+ phib_real = zeros( dims(1)*Nz ,1);
+ phib_imag= phib_real;
+ b_angle = phib_real;
+
+ midpoint = floor((dims(1)*Nz )/2)+1;
+ for ip =1: dims(1)
+ start_ = (ip -1)*Nz +1;
+ end_ = ip*Nz;
+ phib_real(start_:end_) = phi.real(ip,iky,:);
+ phib_imag(start_:end_) = phi.imaginary(ip,iky, :);
+ end
+
+ % Define ballooning angle
+ idx = -Nkx:1:Nkx;
+ for ip =1: dims(1)
+ for iz=1:Nz
+ ii = Nz*(ip -1) + iz;
+ b_angle(ii) =coordz(iz) + 2*pi*idx(ip);
+ end
+ end
+
+ % normalize real and imaginary parts at chi =0
+ [~,idxLFS] = min(abs(b_angle -0));
+ phib = phib_real(:) + 1i * phib_imag(:);
+ % Normalize to the outermid plane
+ phib_norm = phib(:);% / phib( idxLFS) ;
+ phib_real_norm(:) = real( phib_norm);%phib_real(:)/phib_real(idxLFS);
+ phib_imag_norm(:) = imag( phib_norm);%phib_imag(:)/ phib_imag(idxLFS);
+%
+%
+% figure; hold all;
+% plot(b_angle/pi, phib_real_norm,'-b');
+% plot(b_angle/pi, phib_imag_norm ,'-r');
+% plot(b_angle/pi, sqrt(phib_real_norm .^2 + phib_imag_norm.^2),'--k');
+% legend('real','imag','norm')
+% xlabel('$\chi / \pi$')
+% ylabel('$\phi_B (\chi)$');
+% % title(['$(P,J) =(',num2str(pmax),', ', num2str(jmax),'$), $\nu =',num2str(nu),...
+% % '$, $\epsilon = ',num2str(epsilon),'$, $k_y = ', num2str(coordky( iky)),'$, $q =',num2str(safety_fac),'$, $s = ', num2str(shear),'$, $R_N = ', ...
+% % num2str(Rn),'$, $R_{T_i} = ', num2str(RTi),'$, $N_z =',num2str(Nz),'$']);
+% title(['molix, t=',num2str(coordtime(idxte))]);
+% %set(gca,'Yscale','log')
+ %
+
+% %Check symmetry of the mode at the outter mid plane
+% figure; hold all;
+% right = phib_real(midpoint+1:end);
+% left = fliplr(phib_real(1:midpoint-1)');
+% up_down_symm = right(1:end) - left(1:end-1)';
+% %plot(b_angle(midpoint+1:end)/pi,phib_real(midpoint+1:end),'-xb');
+% plot(b_angle(midpoint+1:end)/pi,up_down_symm ,'-xb');
+ %plot(abs(b_angle(1:midpoint-1))/pi,phib_real(1:midpoint-1),'-xb');
+ %
+ %
+ % figure; hold all
+ % plot(b_angle/pi, phib_imag.^2 + phib_real.^2 ,'xk');
+ % %set(gca,'Yscale','log')
+ % xlabel('$\chi / \pi$')
+ % ylabel('$\phi_B (\chi)$');
+ % title(['$(P,J) =(',num2str(pmax),', ', num2str(jmax),'$), $\nu =',num2str(nu),'$, $\epsilon = ',num2str(epsilon),'$, $q =',num2str(safety_fac),'$, $s = ', num2str(shear),'$, $k_y =',num2str(ky),'$']);
+end
+
+end
\ No newline at end of file
diff --git a/wk/plot_cosol_mat.m b/wk/plot_cosol_mat.m
index 6c5dcf67..e410d054 100644
--- a/wk/plot_cosol_mat.m
+++ b/wk/plot_cosol_mat.m
@@ -93,7 +93,7 @@ subplot(224)
%% Single eigenvalue analysis
% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_20_P_4_J_2_N_50_kpm_4.0/scanfiles_00005/self.0.h5';
-mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_20_P_6_J_3_N_50_kpm_4.0/scanfiles_00036self.0.h5';
+mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_20_P_6_J_3_N_50_kpm_4.0/scanfiles_00042/self.0.h5';
matidx = 01;
diff --git a/wk/plot_phi_ballooning.m b/wk/plot_phi_ballooning.m
index 24011c82..47bd05e8 100644
--- a/wk/plot_phi_ballooning.m
+++ b/wk/plot_phi_ballooning.m
@@ -1,5 +1,6 @@
-phi_real=(real(PHI(:,:,:,end)));
-phi_imag=(imag(PHI(:,:,:,end)));
+[~,it] = min(abs(Ts3D - time_2_plot));
+phi_real=(real(PHI(:,:,:,it)));
+phi_imag=(imag(PHI(:,:,:,it)));
% Apply baollooning tranform
for iky=2
dims = size(phi_real);
@@ -29,20 +30,22 @@ for iky=2
[~,idxLFS] = min(abs(b_angle -0));
phib = phib_real(:) + 1i * phib_imag(:);
% Normalize to the outermid plane
- phib_norm = phib / phib( idxLFS) ;
- phib_real_norm = real( phib_norm);%phib_real(:)/phib_real(idxLFS);
- phib_imag_norm = imag( phib_norm);%phib_imag(:)/ phib_imag(idxLFS);
+ phib_norm = phib(:);% / phib( idxLFS) ;
+ phib_real_norm(:) = real( phib_norm);%phib_real(:)/phib_real(idxLFS);
+ phib_imag_norm(:) = imag( phib_norm);%phib_imag(:)/ phib_imag(idxLFS);
figure; hold all;
plot(b_angle/pi, phib_real_norm,'-b');
plot(b_angle/pi, phib_imag_norm ,'-r');
- plot(b_angle/pi, phib_real_norm .^2 + phib_imag_norm.^2,'-k');
+ plot(b_angle/pi, sqrt(phib_real_norm .^2 + phib_imag_norm.^2),'-k');
+ legend('real','imag','norm')
xlabel('$\chi / \pi$')
ylabel('$\phi_B (\chi)$');
- title(['$(P,J) =(',num2str(PMAXI),', ', num2str(JMAXI),'$), $\nu =',num2str(NU),...
- '$, $\epsilon = ',num2str(eps),'$, $k_y = ', num2str(ky( iky)),'$, $q =',num2str(q0),'$, $s = ', num2str(shear),'$, $R_N = ', ...
- num2str(K_N),'$, $R_{T_i} = ', num2str(K_T),'$, $N_z =',num2str(Nz),'$']);
+ title(['HeLaZ(-) molix(o) benchmark, t=',num2str(Ts3D(it))]);
+% title(['HeLaZ,$(P,J) =(',num2str(PMAXI),', ', num2str(JMAXI),'$), $\nu =',num2str(NU),...
+% '$, $\epsilon = ',num2str(eps),'$, $k_y = ', num2str(ky( iky)),'$, $q =',num2str(q0),'$, $s = ', num2str(shear),'$, $R_N = ', ...
+% num2str(K_N),'$, $R_{T_i} = ', num2str(K_T),'$, $N_z =',num2str(Nz),'$']);
%set(gca,'Yscale','log')
%
diff --git a/wk/shearless_linear_fluxtube.m b/wk/shearless_linear_fluxtube.m
index 70311b3b..4fa96d08 100644
--- a/wk/shearless_linear_fluxtube.m
+++ b/wk/shearless_linear_fluxtube.m
@@ -21,11 +21,11 @@ q0 = 2.7; % safety factor
shear = 0.0; % magnetic shear
eps = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 5; % Maximal time unit
+TMAX = 10; % Maximal time unit
DT = 1e-3; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
-SPS3D = 1; % Sampling per time unit for 2D arrays
+SPS3D = 10; % Sampling per time unit for 2D arrays
SPS5D = 1/100; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
@@ -87,10 +87,10 @@ for i = 1:Nparam
JMAXE = JA(i); JMAXI = JA(i);
DT = DTA(i);
setup
- system(['rm fort*.90']);
+% system(['rm fort*.90']);
% Run linear simulation
if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0 > out.txt; cd ../../../wk'])
end
% Load and process results
%%
@@ -98,10 +98,6 @@ for i = 1:Nparam
load_results
end
-if 1
-%% Plot
-plot_phi_ballooning
-end
end
if 0
--
GitLab