diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m index d548c76..377c80f 100644 --- a/matlab/plot/plot_metric.m +++ b/matlab/plot/plot_metric.m @@ -1,30 +1,41 @@ -function [ fig ] = plot_metric( data ) +function [ fig ] = plot_metric( data, options ) names = {'Jacobian','gradxB','gradyB','gradzB','gradz_coeff',... 'gxx','gxy','gyy','gyz','gzz','hatB','hatR','hatZ'}; geo_arrays = zeros(2,data.Nz,numel(names)); for i_ = 1:numel(names) namae = names{i_}; geo_arrays(:,:,i_) = h5read(data.outfilenames{end},['/data/metric/',namae])'; end -fig = figure; -subplot(311) - for i = 1:5 - plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on; - end - xlim([min(data.z),max(data.z)]); legend('show'); title('GYACOMO geometry'); +NPLOT = options.SHOW_FLUXSURF + options.SHOW_METRICS; +if NPLOT > 0 + fig = figure; + if options.SHOW_METRICS + subplot(311) + for i = 1:5 + plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on; + end + xlim([min(data.z),max(data.z)]); legend('show'); title('GYACOMO geometry'); -subplot(312) - for i = 6:10 - plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on; - end - xlim([min(data.z),max(data.z)]); legend('show'); + subplot(312) + for i = 6:10 + plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on; + end + xlim([min(data.z),max(data.z)]); legend('show'); -subplot(313) - for i = 11:13 - plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on; + subplot(313) + for i = 11:13 + plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on; + end + xlim([min(data.z),max(data.z)]); legend('show'); end - xlim([min(data.z),max(data.z)]); legend('show'); + if options.SHOW_FLUXSURF + R = [geo_arrays(1,:,12),geo_arrays(1,1,12)]; + Z = [geo_arrays(1,:,13),geo_arrays(1,1,13)]; + plot(R,Z); + axis equal + end +end end diff --git a/matlab/setup.m b/matlab/setup.m index 2661841..3b8be6c 100644 --- a/matlab/setup.m +++ b/matlab/setup.m @@ -1,187 +1,190 @@ %% _______________________________________________________________________ SIMDIR = ['../results/',SIMID,'/']; % Grid parameters GRID.pmaxe = PMAXE; % Electron Hermite moments GRID.jmaxe = JMAXE; % Electron Laguerre moments GRID.pmaxi = PMAXI; % Ion Hermite moments GRID.jmaxi = JMAXI; % Ion Laguerre moments GRID.Nx = NX; % x grid resolution GRID.Lx = LX; % x length GRID.Nexc = NEXC; % to extend Lx when s>0 GRID.Ny = NY; % y '' GRID.Ly = LY; % y '' GRID.Nz = NZ; % z resolution GRID.Npol = NPOL; % z resolution if SG; GRID.SG = '.true.'; else; GRID.SG = '.false.';end; % Geometry GEOM.geom = ['''',GEOMETRY,'''']; GEOM.q0 = Q0; % q factor GEOM.shear = SHEAR; % shear GEOM.eps = EPS; % inverse aspect ratio +GEOM.kappa = KAPPA; % elongation +GEOM.delta = DELTA; % triangularity +GEOM.zeta = ZETA; % squareness % Model parameters MODEL.CLOS = CLOS; MODEL.NL_CLOS = NL_CLOS; MODEL.LINEARITY = ['''',LINEARITY,'''']; MODEL.KIN_E = KIN_E; if KIN_E; MODEL.KIN_E = '.true.'; else; MODEL.KIN_E = '.false.';end; MODEL.beta = BETA; MODEL.mu_x = MU_X; MODEL.mu_y = MU_Y; MODEL.N_HD = N_HD; MODEL.mu_z = MU_Z; MODEL.mu_p = MU_P; MODEL.mu_j = MU_J; MODEL.nu = NU; % hyper diffusive coefficient nu for HW % temperature ratio T_a/T_e MODEL.tau_e = TAU; MODEL.tau_i = TAU; % mass ratio sqrt(m_a/m_i) MODEL.sigma_e = SIGMA_E; MODEL.sigma_i = 1.0; % charge q_a/e MODEL.q_e =-1.0; MODEL.q_i = 1.0; if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end; % gradients L_perp/L_x MODEL.K_Ni = K_Ni; MODEL.K_Ne = K_Ne; MODEL.K_Ti = K_Ti; MODEL.K_Te = K_Te; MODEL.GradB = GRADB; % Magnetic gradient MODEL.CurvB = CURVB; % Magnetic curvature MODEL.lambdaD = LAMBDAD; % Collision parameters COLL.collision_model = ['''',CO,'''']; if (GKCO); COLL.gyrokin_CO = '.true.'; else; COLL.gyrokin_CO = '.false.';end; if (ABCO); COLL.interspecies = '.true.'; else; COLL.interspecies = '.false.';end; COLL.mat_file = '''null'''; switch CO case 'SG' COLL.mat_file = '''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'''; % COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5'''; % COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5'''; case 'LR' COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'''; case 'LD' % COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5'''; % COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5'''; % COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5'''; COLL.mat_file = '''../../../iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5'''; end COLL.coll_kcut = COLL_KCUT; % Time integration and intialization parameters TIME_INTEGRATION.numerical_scheme = '''RK4'''; INITIAL.INIT_OPT = ['''',INIT_OPT,'''']; INITIAL.ACT_ON_MODES = ['''',ACT_ON_MODES,'''']; INITIAL.init_background = BCKGD0; INITIAL.init_noiselvl = NOISE0; INITIAL.iseed = 42; % Naming and creating input file CONAME = CO; if GKCO CONAME = [CONAME,'GK']; else CONAME = [CONAME,'DK']; end if ~ABCO CONAME = [CONAME,'aa']; end if (CLOS == 0); CLOSNAME = 'Trunc.'; elseif(CLOS == 1); CLOSNAME = 'Clos. 1'; elseif(CLOS == 2); CLOSNAME = 'Clos. 2'; end % Hermite-Laguerre degrees naming if (PMAXE == PMAXI) && (JMAXE == JMAXI) HLdeg_ = ['_',num2str(PMAXE+1),'x',num2str(JMAXE+1)]; else HLdeg_ = ['_Pe_',num2str(PMAXE+1),'_Je_',num2str(JMAXE+1),... '_Pi_',num2str(PMAXI+1),'_Ji_',num2str(JMAXI+1)]; end % temp. dens. drives drives_ = []; if abs(K_Ni) > 0; drives_ = [drives_,'_kN_',num2str(K_Ni)]; end; if abs(K_Ti) > 0; drives_ = [drives_,'_kT_',num2str(K_Ti)]; end; % collision coll_ = ['_nu_%1.1e_',CONAME]; coll_ = sprintf(coll_,NU); % nonlinear lin_ = []; if ~LINEARITY; lin_ = '_lin'; end adiabe_ = []; if ~KIN_E; adiabe_ = '_adiabe'; end % resolution and boxsize res_ = [num2str(GRID.Nx),'x',num2str(GRID.Ny)]; if (LX ~= LY) geo_ = ['_Lx_',num2str(LX),'_Ly_',num2str(LY)]; else geo_ = ['_L_',num2str(LX)]; end if (NZ > 1) %3D case res_ = [res_,'x',num2str(NZ)]; if abs(Q0) > 0 geo_ = [geo_,'_q0_',num2str(Q0)]; end if abs(EPS) > 0 geo_ = [geo_,'_e_',num2str(EPS)]; end if abs(SHEAR) > 0 geo_ = [geo_,'_s_',num2str(SHEAR)]; end end switch GEOMETRY case 'circular' geo_ = [geo_,'_circ_']; end % put everything together in the param character chain u_ = '_'; % underscore variable PARAMS = [res_,HLdeg_,geo_,drives_,coll_,lin_,adiabe_]; BASIC.RESDIR = [SIMDIR,PARAMS,'/']; BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',SIMDIR(4:end),PARAMS,'/']; BASIC.PARAMS = PARAMS; BASIC.SIMID = SIMID; BASIC.nrun = 1e8; BASIC.dt = DT; BASIC.tmax = TMAX; %time normalized to 1/omega_pe BASIC.maxruntime = str2num(CLUSTER.TIME(1:2))*3600 ... + str2num(CLUSTER.TIME(4:5))*60 ... + str2num(CLUSTER.TIME(7:8)); % Outputs parameters OUTPUTS.nsave_0d = floor(1.0/SPS0D/DT); OUTPUTS.nsave_1d = -1; OUTPUTS.nsave_2d = floor(1.0/SPS2D/DT); OUTPUTS.nsave_3d = floor(1.0/SPS3D/DT); OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT); if W_DOUBLE; OUTPUTS.write_doubleprecision = '.true.'; else; OUTPUTS.write_doubleprecision = '.false.';end; if W_GAMMA; OUTPUTS.write_gamma = '.true.'; else; OUTPUTS.write_gamma = '.false.';end; if W_HF; OUTPUTS.write_hf = '.true.'; else; OUTPUTS.write_hf = '.false.';end; if W_PHI; OUTPUTS.write_phi = '.true.'; else; OUTPUTS.write_phi = '.false.';end; if W_NA00; OUTPUTS.write_Na00 = '.true.'; else; OUTPUTS.write_Na00 = '.false.';end; if W_NAPJ; OUTPUTS.write_Napj = '.true.'; else; OUTPUTS.write_Napj = '.false.';end; if W_SAPJ; OUTPUTS.write_Sapj = '.true.'; else; OUTPUTS.write_Sapj = '.false.';end; if W_DENS; OUTPUTS.write_dens = '.true.'; else; OUTPUTS.write_dens = '.false.';end; if W_TEMP; OUTPUTS.write_temp = '.true.'; else; OUTPUTS.write_temp = '.false.';end; OUTPUTS.job2load = JOB2LOAD; %% Create directories if ~exist(SIMDIR, 'dir') mkdir(SIMDIR) end if ~exist(BASIC.RESDIR, 'dir') mkdir(BASIC.RESDIR) end if ~exist(BASIC.MISCDIR, 'dir') mkdir(BASIC.MISCDIR) end %% Compile and WRITE input file INPUT = write_fort90(OUTPUTS,GRID,GEOM,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC); nproc = 1; MAKE = 'cd ..; make; cd wk'; % system(MAKE); %% disp(['Set up ',SIMID]); disp([res_,geo_,HLdeg_]); if JOB2LOAD>=0 disp(['- restarting from JOBNUM = ',num2str(JOB2LOAD)]); else disp(['- starting from T = 0']); end diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m index d42c780..29304ef 100644 --- a/matlab/write_fort90.m +++ b/matlab/write_fort90.m @@ -1,105 +1,108 @@ function [INPUT] = write_fort90(OUTPUTS,GRID,GEOM,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC) % Write the input script "fort.90" with desired parameters INPUT = ['fort_',sprintf('%2.2d',OUTPUTS.job2load+1),'.90']; fid = fopen(INPUT,'wt'); fprintf(fid,'&BASIC\n'); fprintf(fid,[' nrun = ', num2str(BASIC.nrun),'\n']); fprintf(fid,[' dt = ', num2str(BASIC.dt),'\n']); fprintf(fid,[' tmax = ', num2str(BASIC.tmax),'\n']); fprintf(fid,[' maxruntime = ', num2str(BASIC.maxruntime),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&GRID\n'); fprintf(fid,[' pmaxe = ', num2str(GRID.pmaxe),'\n']); fprintf(fid,[' jmaxe = ', num2str(GRID.jmaxe),'\n']); fprintf(fid,[' pmaxi = ', num2str(GRID.pmaxi),'\n']); fprintf(fid,[' jmaxi = ', num2str(GRID.jmaxi),'\n']); fprintf(fid,[' Nx = ', num2str(GRID.Nx),'\n']); fprintf(fid,[' Lx = ', num2str(GRID.Lx),'\n']); fprintf(fid,[' Nexc = ', num2str(GRID.Nexc),'\n']); fprintf(fid,[' Ny = ', num2str(GRID.Ny),'\n']); fprintf(fid,[' Ly = ', num2str(GRID.Ly),'\n']); fprintf(fid,[' Nz = ', num2str(GRID.Nz),'\n']); fprintf(fid,[' Npol = ', num2str(GRID.Npol),'\n']); fprintf(fid,[' SG = ', GRID.SG,'\n']); fprintf(fid,'/\n'); fprintf(fid,'&GEOMETRY\n'); fprintf(fid,[' geom = ', GEOM.geom,'\n']); fprintf(fid,[' q0 = ', num2str(GEOM.q0),'\n']); fprintf(fid,[' shear = ', num2str(GEOM.shear),'\n']); fprintf(fid,[' eps = ', num2str(GEOM.eps),'\n']); +fprintf(fid,[' kappa = ', num2str(GEOM.kappa),'\n']); +fprintf(fid,[' delta = ', num2str(GEOM.delta),'\n']); +fprintf(fid,[' zeta = ', num2str(GEOM.zeta),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&OUTPUT_PAR\n'); fprintf(fid,[' nsave_0d = ', num2str(OUTPUTS.nsave_0d),'\n']); fprintf(fid,[' nsave_1d = ', num2str(OUTPUTS.nsave_1d),'\n']); fprintf(fid,[' nsave_2d = ', num2str(OUTPUTS.nsave_2d),'\n']); fprintf(fid,[' nsave_3d = ', num2str(OUTPUTS.nsave_3d),'\n']); fprintf(fid,[' nsave_5d = ', num2str(OUTPUTS.nsave_5d),'\n']); fprintf(fid,[' write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']); fprintf(fid,[' write_gamma = ', OUTPUTS.write_gamma,'\n']); fprintf(fid,[' write_hf = ', OUTPUTS.write_hf,'\n']); fprintf(fid,[' write_phi = ', OUTPUTS.write_phi,'\n']); fprintf(fid,[' write_Na00 = ', OUTPUTS.write_Na00,'\n']); fprintf(fid,[' write_Napj = ', OUTPUTS.write_Napj,'\n']); fprintf(fid,[' write_Sapj = ', OUTPUTS.write_Sapj,'\n']); fprintf(fid,[' write_dens = ', OUTPUTS.write_dens,'\n']); fprintf(fid,[' write_temp = ', OUTPUTS.write_temp,'\n']); fprintf(fid,[' job2load = ', num2str(OUTPUTS.job2load),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&MODEL_PAR\n'); fprintf(fid,' ! Collisionality\n'); fprintf(fid,[' CLOS = ', num2str(MODEL.CLOS),'\n']); fprintf(fid,[' NL_CLOS = ', num2str(MODEL.NL_CLOS),'\n']); fprintf(fid,[' LINEARITY = ', MODEL.LINEARITY,'\n']); fprintf(fid,[' KIN_E = ', MODEL.KIN_E,'\n']); fprintf(fid,[' mu_x = ', num2str(MODEL.mu_x),'\n']); fprintf(fid,[' mu_y = ', num2str(MODEL.mu_y),'\n']); fprintf(fid,[' N_HD = ', num2str(MODEL.N_HD),'\n']); fprintf(fid,[' mu_z = ', num2str(MODEL.mu_z),'\n']); fprintf(fid,[' mu_p = ', num2str(MODEL.mu_p),'\n']); fprintf(fid,[' mu_j = ', num2str(MODEL.mu_j),'\n']); fprintf(fid,[' nu = ', num2str(MODEL.nu),'\n']); fprintf(fid,[' tau_e = ', num2str(MODEL.tau_e),'\n']); fprintf(fid,[' tau_i = ', num2str(MODEL.tau_i),'\n']); fprintf(fid,[' sigma_e = ', num2str(MODEL.sigma_e),'\n']); fprintf(fid,[' sigma_i = ', num2str(MODEL.sigma_i),'\n']); fprintf(fid,[' q_e = ', num2str(MODEL.q_e),'\n']); fprintf(fid,[' q_i = ', num2str(MODEL.q_i),'\n']); fprintf(fid,[' K_Ne = ', num2str(MODEL.K_Ne),'\n']); fprintf(fid,[' K_Ni = ', num2str(MODEL.K_Ni),'\n']); fprintf(fid,[' K_Te = ', num2str(MODEL.K_Te),'\n']); fprintf(fid,[' K_Ti = ', num2str(MODEL.K_Ti),'\n']); fprintf(fid,[' GradB = ', num2str(MODEL.GradB),'\n']); fprintf(fid,[' CurvB = ', num2str(MODEL.CurvB),'\n']); fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']); fprintf(fid,[' beta = ', num2str(MODEL.beta),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&COLLISION_PAR\n'); fprintf(fid,[' collision_model = ', COLL.collision_model,'\n']); fprintf(fid,[' gyrokin_CO = ', COLL.gyrokin_CO,'\n']); fprintf(fid,[' interspecies = ', COLL.interspecies,'\n']); fprintf(fid,[' mat_file = ', COLL.mat_file,'\n']); fprintf(fid,[' collision_kcut = ', num2str(COLL.coll_kcut),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&INITIAL_CON\n'); fprintf(fid,[' INIT_OPT = ', INITIAL.INIT_OPT,'\n']); fprintf(fid,[' ACT_ON_MODES = ', INITIAL.ACT_ON_MODES,'\n']); fprintf(fid,[' init_background = ', num2str(INITIAL.init_background),'\n']); fprintf(fid,[' init_noiselvl = ', num2str(INITIAL.init_noiselvl),'\n']); fprintf(fid,[' iseed = ', num2str(INITIAL.iseed),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&TIME_INTEGRATION_PAR\n'); fprintf(fid,[' numerical_scheme = ', TIME_INTEGRATION.numerical_scheme,'\n']); fprintf(fid,'/'); fclose(fid); system(['cp fort*.90 ',BASIC.RESDIR,'/.']); end diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90 index 4921117..6f478b7 100644 --- a/src/geometry_mod.F90 +++ b/src/geometry_mod.F90 @@ -1,432 +1,452 @@ module geometry ! computes geometrical quantities ! Adapted from B.J.Frei MOLIX code (2021) use prec_const use model use grid use array use fields use basic use calculus, ONLY: simpson_rule_z use miller, ONLY: set_miller_parameters, get_miller implicit none PRIVATE ! Geometry input parameters CHARACTER(len=16), & PUBLIC, PROTECTED :: geom REAL(dp), PUBLIC, PROTECTED :: q0 = 1.4_dp ! safety factor REAL(dp), PUBLIC, PROTECTED :: shear = 0._dp ! magnetic field shear REAL(dp), PUBLIC, PROTECTED :: eps = 0.18_dp ! inverse aspect ratio REAL(dp), PUBLIC, PROTECTED :: alpha_MHD = 0 ! shafranov shift effect alpha = -q2 R dbeta/dr ! parameters for Miller geometry REAL(dp), PUBLIC, PROTECTED :: kappa = 1._dp ! elongation REAL(dp), PUBLIC, PROTECTED :: s_kappa = 0._dp ! r normalized derivative skappa = r/kappa dkappa/dr REAL(dp), PUBLIC, PROTECTED :: delta = 0._dp ! triangularity REAL(dp), PUBLIC, PROTECTED :: s_delta = 0._dp ! '' sdelta = r/sqrt(1-delta2) ddelta/dr REAL(dp), PUBLIC, PROTECTED :: zeta = 0._dp ! squareness REAL(dp), PUBLIC, PROTECTED :: s_zeta = 0._dp ! '' szeta = r dzeta/dr ! GENE unused additional parameters for miller_mod REAL(dp), PUBLIC, PROTECTED :: edge_opt = 0 ! meant to redistribute the points in z REAL(dp), PUBLIC, PROTECTED :: major_R = 1 ! major radius REAL(dp), PUBLIC, PROTECTED :: major_Z = 0 ! vertical elevation REAL(dp), PUBLIC, PROTECTED :: dpdx_pm_geom = 0 ! amplitude mag. eq. pressure grad. REAL(dp), PUBLIC, PROTECTED :: C_y = 0 ! defines y coordinate : Cy (q theta - phi) REAL(dp), PUBLIC, PROTECTED :: C_xy = 0 ! defines x coordinate : B = Cxy Vx x Vy ! Geometrical auxiliary variables LOGICAL, PUBLIC, PROTECTED :: SHEARED = .false. ! flag for shear magn. geom or not ! Curvature REAL(dp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z, odd/even p ! Jacobian REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Jacobian ! dimensions: z, odd/even p COMPLEX(dp), PUBLIC, PROTECTED :: iInt_Jacobian ! Inverse integrated Jacobian ! Metric 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 ! Relative magnetic field strength REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB, hatB_NL ! Relative strength of major radius REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ ! Some geometrical coefficients REAL(dp), PUBLIC, DIMENSION(:,:) , ALLOCATABLE :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ] ! Array to map the index of mode (kx,ky,-pi) to (kx+2pi*s*ky,ky,pi) for sheared periodic boundary condition INTEGER, PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ikx_zBC_L, ikx_zBC_R ! Functions PUBLIC :: geometry_readinputs, geometry_outputinputs,& eval_magnetic_geometry, set_ikx_zBC_map CONTAINS SUBROUTINE geometry_readinputs ! Read the input parameters IMPLICIT NONE NAMELIST /GEOMETRY/ geom, q0, shear, eps,& kappa, s_kappa,delta, s_delta, zeta, s_zeta ! For miller READ(lu_in,geometry) IF(shear .NE. 0._dp) SHEARED = .true. END SUBROUTINE geometry_readinputs subroutine eval_magnetic_geometry ! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient implicit none REAL(dp) :: kx,ky COMPLEX(dp), DIMENSION(izs:ize) :: integrant INTEGER :: fid + real(dp) :: G1,G2,G3,Cx,Cy ! Allocate arrays CALL geometry_allocate_mem ! 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 SELECT CASE(geom) CASE('s-alpha') IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry' call eval_salphaB_geometry CASE('Z-pinch') IF( my_id .eq. 0 ) WRITE(*,*) 'Z-pinch geometry' call eval_zpinch_geometry SHEARED = .FALSE. CASE('miller') IF( my_id .eq. 0 ) WRITE(*,*) 'Miller geometry' 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,& Ckxky,gradz_coeff) CASE DEFAULT STOP 'geometry not recognized!!' END SELECT ENDIF ! ! 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 eo = 0,1 DO iz = izgs,izge DO iky = ikys, ikye 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 ENDDO ENDDO ENDDO + ! Curvature operator (Frei et al. 2022 eq 2.15) + DO iz = izgs,izge + 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) + + DO iky = ikys, ikye + ky = kyarray(iky) + DO ikx= ikxs, ikxe + kx = kxarray(ikx) + Ckxky(iky, ikx, iz,eo) = (Cx*kx + Cy*ky) + ENDDO + ENDDO + ! coefficient in the front of parallel derivative + gradz_coeff(iz,eo) = 1._dp / jacobian(iz,eo) / hatB(iz,eo) + ! Factor in front of the nonlinear term + hatB_NL(iz,eo) = Jacobian(iz,eo)& + *(gxx(iz,eo)*gyy(iz,eo) - gxy(iz,eo)**2)/hatB(iz,eo) + ENDDO ENDDO - ! Factor in front of the nonlinear term - hatB_NL(izgs:izge,0:1) = Jacobian(izgs:izge,0:1)& - *(gxx(izgs:izge,0:1)*gyy(izgs:izge,0:1) - gxy(izgs:izge,0:1)**2)/hatB(izgs:izge,0:1) + ! set the mapping for parallel boundary conditions CALL set_ikx_zBC_map two_third_kpmax = 2._dp/3._dp * MAXVAL(kparray) ! ! Compute the inverse z integrated Jacobian (useful for flux averaging) integrant = Jacobian(izs:ize,0) ! Convert into complex array CALL simpson_rule_z(integrant,iInt_Jacobian) iInt_Jacobian = 1._dp/iInt_Jacobian ! reverse it END SUBROUTINE eval_magnetic_geometry ! !-------------------------------------------------------------------------------- ! SUBROUTINE eval_salphaB_geometry ! evaluate s-alpha geometry model implicit none REAL(dp) :: z, kx, ky, Gx, Gy alpha_MHD = 0._dp parity: DO eo = 0,1 zloop: DO iz = izgs,izge z = zarray(iz,eo) ! metric gxx(iz,eo) = 1._dp gxy(iz,eo) = shear*z - alpha_MHD*SIN(z) 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 dxdR(iz,eo)= COS(z) dxdZ(iz,eo)= SIN(z) ! Relative strengh of radius hatR(iz,eo) = 1._dp + eps*COS(z) hatZ(iz,eo) = 1._dp + eps*SIN(z) ! toroidal coordinates Rc (iz,eo) = hatR(iz,eo) 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) ! 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 ! 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) ENDDO zloop ENDDO parity END SUBROUTINE eval_salphaB_geometry ! !-------------------------------------------------------------------------------- ! SUBROUTINE eval_zpinch_geometry implicit none REAL(dp) :: z, kx, ky, alpha_MHD alpha_MHD = 0._dp parity: DO eo = 0,1 zloop: DO iz = izgs,izge z = zarray(iz,eo) ! metric gxx(iz,eo) = 1._dp gxy(iz,eo) = 0._dp gxz(iz,eo) = 0._dp gyy(iz,eo) = 1._dp gyz(iz,eo) = 0._dp gzz(iz,eo) = 1._dp dxdR(iz,eo)= COS(z) dxdZ(iz,eo)= SIN(z) ! Relative strengh of radius hatR(iz,eo) = 1._dp hatZ(iz,eo) = 1._dp ! toroidal coordinates Rc (iz,eo) = hatR(iz,eo) phic(iz,eo) = z Zc (iz,eo) = hatZ(iz,eo) ! Jacobian Jacobian(iz,eo) = 1._dp ! Relative strengh of modulus of B hatB (iz,eo) = 1._dp hatB_NL(iz,eo) = 1._dp ! Derivative of the magnetic field strenght gradxB(iz,eo) = 0._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 ! Curvature operator DO iky = ikys, ikye ky = kyarray(iky) DO ikx= ikxs, ikxe kx = kxarray(ikx) Ckxky(iky, ikx, iz,eo) = -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) ENDDO zloop ENDDO parity END SUBROUTINE eval_zpinch_geometry ! !-------------------------------------------------------------------------------- ! subroutine eval_1D_geometry ! evaluate 1D perp geometry model implicit none REAL(dp) :: z, kx, ky parity: DO eo = 0,1 zloop: DO iz = izs,ize z = zarray(iz,eo) ! metric gxx(iz,eo) = 1._dp gxy(iz,eo) = 0._dp gyy(iz,eo) = 1._dp ! Relative strengh of radius hatR(iz,eo) = 1._dp ! Jacobian Jacobian(iz,eo) = 1._dp ! Relative strengh of modulus of B hatB(iz,eo) = 1._dp ! Curvature operator DO iky = ikys, ikye ky = kyarray(iky) DO ikx= ikxs, ikxe kx = kxarray(ikx) Ckxky(ikx, iky, iz,eo) = -kx ! .. 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 ENDDO zloop ENDDO parity END SUBROUTINE eval_1D_geometry ! !-------------------------------------------------------------------------------- ! SUBROUTINE set_ikx_zBC_map IMPLICIT NONE REAL :: shift, kx_shift ! For periodic CHI BC or 0 dirichlet LOGICAL :: PERIODIC_CHI_BC = .TRUE. ALLOCATE(ikx_zBC_R(ikys:ikye,ikxs:ikxe)) ALLOCATE(ikx_zBC_L(ikys:ikye,ikxs:ikxe)) ! No periodic connection for extension of the domain IF(Nexc .GT. 1) PERIODIC_CHI_BC = .TRUE. !! No shear case (simple id mapping) !3 | 1 2 3 4 5 6 | ky = 3 dky !2 ky | 1 2 3 4 5 6 | ky = 2 dky !1 A | 1 2 3 4 5 6 | ky = 1 dky !0 | -> kx | 1____2____3____4____5____6 | ky = 0 dky !(e.g.) kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=free) DO iky = ikys,ikye DO ikx = ikxs,ikxe ikx_zBC_L(iky,ikx) = ikx ikx_zBC_R(iky,ikx) = ikx ENDDO ENDDO ! Modify connection map only at border of z IF(SHEARED) THEN ! connection map BC of the RIGHT boundary (z=pi*Npol-dz) (even NZ) !3 | 4 x x x 2 3 | ky = 3 dky !2 ky | 3 4 x x 1 2 | ky = 2 dky !1 A | 2 3 4 x 6 1 | ky = 1 dky !0 | -> kx | 1____2____3____4____5____6 | ky = 0 dky !kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=2pi*shear*npol*dky) ! connection map BC of the RIGHT boundary (z=pi*Npol-dz) (ODD NZ) !3 | x x x 2 3 | ky = 3 dky !2 ky | 3 x x 1 2 | ky = 2 dky !1 A | 2 3 x 5 1 | ky = 1 dky !0 | -> kx | 1____2____3____4____5 | ky = 0 dky !kx = 0 0.1 0.2 -0.2 -0.1 (dkx=2pi*shear*npol*dky) IF(contains_zmax) THEN ! Check if the process is at the end of the FT DO iky = ikys,ikye shift = 2._dp*PI*shear*kyarray(iky)*Npol DO ikx = ikxs,ikxe kx_shift = kxarray(ikx) + shift ! We use EPSILON() to treat perfect equality case IF( ((kx_shift-EPSILON(kx_shift)) .GT. kx_max) .AND. (.NOT. PERIODIC_CHI_BC) )THEN ! outside of the frequ domain ikx_zBC_R(iky,ikx) = -99 ELSE ikx_zBC_R(iky,ikx) = ikx+(iky-1)*Nexc IF( ikx_zBC_R(iky,ikx) .GT. Nkx ) & ikx_zBC_R(iky,ikx) = ikx_zBC_R(iky,ikx) - Nkx ENDIF ENDDO ENDDO ENDIF ! connection map BC of the LEFT boundary (z=-pi*Npol) !3 | x 5 6 1 x x | ky = 3 dky !2 ky | 5 6 1 2 x x | ky = 2 dky !1 A | 6 1 2 3 x 5 | ky = 1 dky !0 | -> kx | 1____2____3____4____5____6 | ky = 0 dky !(e.g.) kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=2pi*shear*npol*dky) IF(contains_zmin) THEN ! Check if the process is at the start of the FT DO iky = ikys,ikye shift = 2._dp*PI*shear*kyarray(iky)*Npol DO ikx = ikxs,ikxe kx_shift = kxarray(ikx) - shift ! We use EPSILON() to treat perfect equality case IF( ((kx_shift+EPSILON(kx_shift)) .LT. kx_min) .AND. (.NOT. PERIODIC_CHI_BC) ) THEN ! outside of the frequ domain ikx_zBC_L(iky,ikx) = -99 ELSE ikx_zBC_L(iky,ikx) = ikx-(iky-1)*Nexc IF( ikx_zBC_L(iky,ikx) .LT. 1 ) & ikx_zBC_L(iky,ikx) = ikx_zBC_L(iky,ikx) + Nkx ENDIF ENDDO ENDDO ENDIF ELSE ENDIF END SUBROUTINE set_ikx_zBC_map ! !-------------------------------------------------------------------------------- ! SUBROUTINE geometry_allocate_mem ! Curvature and geometry CALL allocate_array( Ckxky, ikys,ikye, ikxs,ikxe,izgs,izge,0,1) CALL allocate_array( Jacobian,izgs,izge, 0,1) CALL allocate_array( gxx,izgs,izge, 0,1) CALL allocate_array( gxy,izgs,izge, 0,1) CALL allocate_array( gxz,izgs,izge, 0,1) 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( hatB,izgs,izge, 0,1) CALL allocate_array( hatB_NL,izgs,izge, 0,1) CALL allocate_array( hatR,izgs,izge, 0,1) CALL allocate_array( hatZ,izgs,izge, 0,1) CALL allocate_array( Rc,izgs,izge, 0,1) CALL allocate_array( phic,izgs,izge, 0,1) CALL allocate_array( Zc,izgs,izge, 0,1) CALL allocate_array( dxdR,izgs,izge, 0,1) CALL allocate_array( dxdZ,izgs,izge, 0,1) call allocate_array(gradz_coeff,izgs,izge, 0,1) CALL allocate_array( kparray, ikys,ikye, ikxs,ikxe,izgs,izge,0,1) END SUBROUTINE geometry_allocate_mem SUBROUTINE geometry_outputinputs(fidres, str) ! Write the input parameters to the results_xx.h5 file USE futils, ONLY: attach USE prec_const IMPLICIT NONE INTEGER, INTENT(in) :: fidres CHARACTER(len=256), INTENT(in) :: str CALL attach(fidres, TRIM(str),"geometry", geom) CALL attach(fidres, TRIM(str), "q0", q0) CALL attach(fidres, TRIM(str), "shear", shear) CALL attach(fidres, TRIM(str), "eps", eps) END SUBROUTINE geometry_outputinputs end module geometry diff --git a/src/miller_mod.F90 b/src/miller_mod.F90 index b33e3a5..eb3a9e5 100644 --- a/src/miller_mod.F90 +++ b/src/miller_mod.F90 @@ -1,667 +1,667 @@ !! This source has been adapted from GENE https://genecode.org/ !! !>Implementation of the local equilibrium model of [R.L. Miller et al., PoP 5, 973 (1998) !>and [J. Candy, PPCF 51, 105009 (2009)] MODULE miller USE prec_const USE basic ! use coordinates,only: gcoor, get_dzprimedz USE grid ! use discretization USE lagrange_interpolation ! use par_in, only: beta, sign_Ip_CW, sign_Bt_CW, Npol USE model implicit none public:: get_miller, set_miller_parameters public:: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta public:: thetaShift public:: mMode, nMode public:: thetak, thetad public:: aSurf, Delta2, Delta3, theta2, theta3, Raxis, Zaxis public:: Deltam, Deltan, s_Deltam, s_Deltan, thetam, thetan public:: cN_m, sN_m, cNdr_m, sNdr_m private real(dp) :: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta INTEGER :: mMode, nMode real(dp) :: thetaShift real(dp) :: thetak, thetad real(dp) :: aSurf, Delta2, Delta3, theta2, theta3, Raxis, Zaxis real(dp) :: Deltam, Deltan, s_Deltam, s_Deltan, thetam, thetan INTEGER, PARAMETER :: IND_M=32 real(dp), DIMENSION(0:IND_M-1) :: cN_m, sN_m, cNdr_m, sNdr_m CONTAINS !>Set defaults for miller parameters subroutine set_miller_parameters(kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_) real(dp), INTENT(IN) :: kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_ rho = -1.0 kappa = kappa_ s_kappa = s_kappa_ delta = delta_ s_delta = s_delta_ zeta = zeta_ s_zeta = s_zeta_ drR = 0.0 drZ = 0.0 thetak = 0.0 thetad = 0.0 aSurf = 0.54 Delta2 = 1.0 Delta3 = 1.0 theta2 = 0.0 theta3 = 0.0 Raxis = 1.0 Zaxis = 0.0 mMode = 2 nMode = 3 Deltam = 1.0 Deltan = 1.0 s_Deltam = 0.0 s_Deltan = 0.0 thetam = 0.0 thetan = 0.0 cN_m = 0.0 sN_m = 0.0 cNdr_m = 0.0 sNdr_m = 0.0 end subroutine set_miller_parameters !>Get Miller metric, magnetic field, jacobian etc. subroutine get_miller(trpeps,major_R,major_Z,q0,shat,amhd,edge_opt,& C_y,C_xy,dpdx_pm_geom,gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,dBdx_,dBdy_,& Bfield_,jacobian_,dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_,Ckxky_,gradz_coeff_) !!!!!!!!!!!!!!!! GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real(dp), INTENT(INOUT) :: trpeps ! eps in gyacomo (inverse aspect ratio) real(dp), INTENT(INOUT) :: major_R ! major radius real(dp), INTENT(INOUT) :: major_Z ! major Z real(dp), INTENT(INOUT) :: q0 ! safetyfactor real(dp), INTENT(INOUT) :: shat ! safetyfactor real(dp), INTENT(INOUT) :: amhd ! alpha mhd real(dp), INTENT(INOUT) :: edge_opt ! alpha mhd real(dp), INTENT(INOUT) :: dpdx_pm_geom ! amplitude mag. eq. pressure grad. real(dp), INTENT(INOUT) :: C_y, C_xy real(dp), dimension(izgs:izge,0:1), INTENT(INOUT) :: & gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,& dBdx_,dBdy_,Bfield_,jacobian_,& dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_, & gradz_coeff_ real(dp), dimension(ikys:ikye,ikxs:ikxe,izgs:izge,0:1), INTENT(INOUT) :: Ckxky_ ! No parameter in gyacomo yet integer :: ikxgs =1 ! the left ghost gpdisc%pi1gl real(dp) :: sign_Ip_CW=1 ! current sign (only normal current) real(dp) :: sign_Bt_CW=1 ! current sign (only normal current) !!!!!!!!!!!!!! END GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Auxiliary variables for curvature computation real(dp) :: G1,G2,G3,Cx,Cy,ky,kx integer:: np, np_s, Npol_ext, Npol_s real(dp), dimension(500*(Npol+2)):: R,Z,R_rho,Z_rho,R_theta,Z_theta,R_theta_theta,Z_theta_theta,dlp,Rc,cosu,sinu,Bphi real(dp), dimension(500*(Npol+2)):: drRcirc, drRelong, drRelongTilt, drRtri, drRtriTilt, drZcirc, drZelong, drZelongTilt real(dp), dimension(500*(Npol+2)):: drZtri, drZtriTilt, dtdtRcirc, dtdtRelong, dtdtRelongTilt, dtdtRtri, dtdtRtriTilt real(dp), dimension(500*(Npol+2)):: dtdtZcirc, dtdtZelong, dtdtZelongTilt, dtdtZtri, dtdtZtriTilt, dtRcirc, dtRelong real(dp), dimension(500*(Npol+2)):: dtRelongTilt, dtRtri, dtRtriTilt, dtZcirc, dtZelong, dtZelongTilt, dtZtri, dtZtriTilt real(dp), dimension(500*(Npol+2)):: Rcirc, Relong, RelongTilt, Rtri, RtriTilt, Zcirc, Zelong, ZelongTilt, Ztri, ZtriTilt real(dp), dimension(500*(Npol+2)):: drrShape, drrAng, drxAng, dryAng, dtdtrShape, dtdtrAng, dtdtxAng real(dp), dimension(500*(Npol+2)):: dtdtyAng, dtrShape, dtrAng, dtxAng, dtyAng, rShape, rAng, xAng, yAng real(dp), dimension(500*(Npol+2)):: theta, thAdj, J_r, B, Bp, D0, D1, D2, D3, nu, chi, psi1, nu1 real(dp), dimension(500*(Npol+2)):: tmp_reverse, theta_reverse, tmp_arr real(dp), dimension(500*(Npol+1)):: theta_s, thAdj_s, chi_s, theta_s_reverse real(dp), dimension(500*(Npol+1)):: R_s, Z_s, R_theta_s, Z_theta_s, Rc_s, cosu_s, sinu_s, Bphi_s, B_s, Bp_s, dlp_s real(dp), dimension(500*(Npol+1)):: dtRcirc_s, dtRelong_s, dtRelongTilt_s, dtRtri_s, dtRtriTilt_s, dtZcirc_s real(dp), dimension(500*(Npol+1)):: dtZelong_s, dtZelongTilt_s, dtZtri_s, dtZtriTilt_s, Rcirc_s, Relong_s, RelongTilt_s real(dp), dimension(500*(Npol+1)):: Rtri_s, RtriTilt_s, Zcirc_s, Zelong_s, ZelongTilt_s, Ztri_s, ZtriTilt_s, dtrShape_s real(dp), dimension(500*(Npol+1)):: dtrAng_s, dtxAng_s, dtyAng_s, rShape_s, rAng_s, xAng_s, yAng_s real(dp), dimension(500*(Npol+1)):: psi1_s, nu1_s, dchidx_s, dB_drho_s, dB_dl_s, dnu_drho_s, dnu_dl_s, grad_nu_s real(dp), dimension(500*(Npol+1)):: gxx, gxy, gxz, gyy, gyz, gzz, dtheta_dchi_s, dBp_dchi_s, jacobian, dBdx, dBdz real(dp), dimension(500*(Npol+1)):: g_xx, g_xy, g_xz, g_yy, g_yz, g_zz, tmp_arr_s, dxdR_s, dxdZ_s, K_x, K_y !tmp_arr2 real(dp), dimension(1:Nz):: gxx_out,gxy_out,gxz_out,gyy_out,gyz_out,gzz_out,Bfield_out,jacobian_out, dBdx_out, dBdz_out, chi_out real(dp), dimension(1:Nz):: R_out, Z_out, dxdR_out, dxdZ_out real(dp):: d_inv, drPsi, dxPsi, dq_dx, dq_dpsi, R0, Z0, B0, F, D0_full, D1_full, D2_full, D3_full !real(dp) :: Lnorm, Psi0 ! currently module-wide defined anyway real(dp):: pprime, ffprime, D0_mid, D1_mid, D2_mid, D3_mid, dx_drho, pi, mu_0, dzprimedz real(dp):: rho_a, psiN, drpsiN, CN2, CN3, Rcenter, Zcenter, drRcenter, drZcenter logical:: bMaxShift integer:: i, k, iBmax Npol_ext = Npol+2 Npol_s = Npol+1 np = 500*Npol_ext np_s = 500*Npol_s rho = trpeps*major_R if (rho.le.0.0) stop 'flux surface radius not defined' trpeps = rho/major_R q0 = sign_Ip_CW * sign_Bt_CW * abs(q0) R0=major_R B0=1.0*sign_Bt_CW F=R0*B0 Z0=major_Z pi = acos(-1.0) mu_0=4.0*pi theta=linspace(-pi*Npol_ext,pi*Npol_ext-2._dp*pi*Npol_ext/np,np) d_inv=asin(delta) thetaShift = 0.0 iBmax = 1 !flux surface parametrization thAdj = theta + thetaShift if (zeta/=0.0 .or. s_zeta/=0.0) 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 !dl/dtheta dlp=(R_theta**2+Z_theta**2)**0.5 !curvature radius Rc=dlp**3*(R_theta*Z_theta_theta-Z_theta*R_theta_theta)**(-1) ! some useful quantities (see papers for definition of u) cosu=Z_theta/dlp sinu=-R_theta/dlp !Jacobian J_r = (dPsi/dr) J_psi = (dPsi/dr) / [(nabla fz x nabla psi)* nabla theta] ! = R * (dR/drho dZ/dtheta - dR/dtheta dZ/drho) = R dlp / |nabla r| J_r=R*(R_rho*Z_theta-R_theta*Z_rho) !From definition of q = 1/(2 pi) int (B nabla fz) / (B nabla theta) dtheta: !dPsi/dr = sign_Bt sign_Ip / (2 pi q) int F / R^2 J_r dtheta ! = F / (2 pi |q|) int J_r/R^2 dtheta tmp_arr=J_r/R**2 drPsi=sign_Ip_CW*F/(2.*pi*Npol_ext*q0)*sum(tmp_arr)*2*pi*Npol_ext/np !dlp_int(tmp_arr,1.0) !Poloidal field (Bp = Bvec * nabla l) Bp=sign_Ip_CW * drPsi / J_r * dlp !toroidal field Bphi=F/R !total modulus of Bfield B=sqrt(Bphi**2+Bp**2) bMaxShift = .false. ! if (thetaShift==0.0.and.trim(magn_geometry).ne.'miller_general') then if (thetaShift==0.0) then do i = 2,np-1 if (B(iBmax) x = r dx_drho=1. !drPsi/Psi0*Lnorm*q0 if (my_id==0) write(*,"(A,ES12.4)") 'Using radial coordinate with dx/dr = ',dx_drho dxPsi = drPsi/dx_drho C_y = dxPsi*sign_Ip_CW C_xy = abs(B0*dxPsi/C_y) if (my_id==0) then write(*,"(A,ES12.4,A,ES12.4,A,ES12.4)") & "Setting C_xy = ",C_xy,' C_y = ', C_y," C_x' = ", 1./dxPsi write(*,'(A,ES12.4)') "B_unit/Bref conversion factor = ", q0/rho*drPsi write(*,'(A,ES12.4)') "dPsi/dr = ", drPsi if (thetaShift.ne.0.0) write(*,'(A,ES12.4)') "thetaShift = ", thetaShift endif !--------shear is expected to be defined as rho/q*dq/drho--------! dq_dx=shat*q0/rho/dx_drho dq_dpsi=dq_dx/dxPsi pprime=-amhd/q0**2/R0/(2*mu_0)*B0**2/drPsi !neg. dpdx normalized to magnetic pressure for pressure term dpdx_pm_geom=amhd/q0**2/R0/dx_drho !first coefficient of psi in varrho expansion psi1 = R*Bp*sign_Ip_CW !integrals for ffprime evaluation do i=1,np tmp_arr=(2./Rc-2.*cosu/R)/(R*psi1**2) D0(i)=-F*dlp_int_ind(tmp_arr,dlp,i) tmp_arr=B**2*R/psi1**3 D1(i)=-dlp_int_ind(tmp_arr,dlp,i)/F tmp_arr=mu_0*R/psi1**3 D2(i)=-dlp_int_ind(tmp_arr,dlp,i)*F tmp_arr=1./(R*psi1) D3(i)=-dlp_int_ind(tmp_arr,dlp,i)*F enddo tmp_arr=(2./Rc-2.*cosu/R)/(R*psi1**2) D0_full=-F*dlp_int(tmp_arr,dlp) tmp_arr=B**2*R/psi1**3 D1_full=-dlp_int(tmp_arr,dlp)/F tmp_arr=mu_0*R/psi1**3 D2_full=-dlp_int(tmp_arr,dlp)*F tmp_arr=1./(R*psi1) D3_full=-dlp_int(tmp_arr,dlp)*F D0_mid=D0(np/2+1) D1_mid=D1(np/2+1) D2_mid=D2(np/2+1) D3_mid=D3(np/2+1) ffprime=-(sign_Ip_CW*dq_dpsi*2.*pi*Npol_ext+D0_full+D2_full*pprime)/D1_full if (my_id==0) then write(*,'(A,ES12.4)') "ffprime = ", ffprime endif D0=D0-D0_mid D1=D1-D1_mid D2=D2-D2_mid nu=D3-D3_mid nu1=psi1*(D0+D1*ffprime+D2*pprime) !straight field line angle defined on equidistant theta grid !alpha = fz + nu = - (q chi - fz) => chi = -nu / q chi=-nu/q0 !correct small scaling error (<0.5%, due to finite integration resolution) chi=chi*(maxval(theta)-minval(theta))/(maxval(chi)-minval(chi)) !new grid equidistant in straight field line angle chi_s = linspace(-pi*Npol_s,pi*Npol_s-2*pi*Npol_s/np_s,np_s) if (sign_Ip_CW.lt.0.0) then !make chi increasing function to not confuse lag3interp tmp_reverse = chi(np:1:-1) theta_reverse = theta(np:1:-1) call lag3interp(theta_reverse,tmp_reverse,np,theta_s,chi_s,np_s) theta_s_reverse = theta_s(np_s:1:-1) else !lag3interp(y_in,x_in,n_in,y_out,x_out,n_out) call lag3interp(theta,chi,np,theta_s,chi_s,np_s) endif dtheta_dchi_s=deriv_fd(theta_s,chi_s,np_s) !arrays equidistant in straight field line angle thAdj_s = theta_s + thetaShift if (zeta/=0.0 .or. s_zeta/=0.0) then R_s = R0 + rho*Cos(thAdj_s + d_inv*Sin(thAdj_s)) Z_s = Z0 + kappa*rho*Sin(thAdj_s + zeta*Sin(2*thAdj_s)) R_theta_s = -(dtheta_dchi_s*rho*(1 + d_inv*Cos(thAdj_s))*Sin(thAdj_s + d_inv*Sin(thAdj_s))) Z_theta_s = dtheta_dchi_s*kappa*rho*(1 + 2*zeta*Cos(2*thAdj_s))*Cos(thAdj_s + zeta*Sin(2*thAdj_s)) else Rcirc_s = rho*Cos(thAdj_s - thetad + thetak) Zcirc_s = rho*Sin(thAdj_s - thetad + thetak) Relong_s = Rcirc_s Zelong_s = Zcirc_s + (-1 + kappa)*rho*Sin(thAdj_s - thetad + thetak) RelongTilt_s = Relong_s*Cos(thetad - thetak) - Zelong_s*Sin(thetad - thetak) ZelongTilt_s = Zelong_s*Cos(thetad - thetak) + Relong_s*Sin(thetad - thetak) Rtri_s = RelongTilt_s - rho*Cos(thAdj_s) + rho*Cos(thAdj_s + delta*Sin(thAdj_s)) Ztri_s = ZelongTilt_s RtriTilt_s = Rtri_s*Cos(thetad) + Ztri_s*Sin(thetad) ZtriTilt_s = Ztri_s*Cos(thetad) - Rtri_s*Sin(thetad) R_s = R0 + RtriTilt_s Z_s = Z0 + ZtriTilt_s dtRcirc_s = -(rho*Sin(thAdj_s - thetad + thetak)) dtZcirc_s = rho*Cos(thAdj_s - thetad + thetak) dtRelong_s = dtRcirc_s dtZelong_s = dtZcirc_s + (-1 + kappa)*rho*Cos(thAdj_s - thetad + thetak) dtRelongTilt_s = dtRelong_s*Cos(thetad - thetak) - dtZelong_s*Sin(thetad - thetak) dtZelongTilt_s = dtZelong_s*Cos(thetad - thetak) + dtRelong_s*Sin(thetad - thetak) dtRtri_s = dtRelongTilt_s + rho*Sin(thAdj_s) & - rho*(1 + delta*Cos(thAdj_s))*Sin(thAdj_s + delta*Sin(thAdj_s)) dtZtri_s = dtZelongTilt_s dtRtriTilt_s = dtRtri_s*Cos(thetad) + dtZtri_s*Sin(thetad) dtZtriTilt_s = dtZtri_s*Cos(thetad) - dtRtri_s*Sin(thetad) R_theta_s = dtheta_dchi_s*dtRtriTilt_s Z_theta_s = dtheta_dchi_s*dtZtriTilt_s endif if (sign_Ip_CW.lt.0.0) then call lag3interp(nu1,theta,np,tmp_arr_s,theta_s_reverse,np_s) nu1_s = tmp_arr_s(np_s:1:-1) call lag3interp(Bp,theta,np,tmp_arr_s,theta_s_reverse,np_s) Bp_s = tmp_arr_s(np_s:1:-1) call lag3interp(dlp,theta,np,tmp_arr_s,theta_s_reverse,np_s) dlp_s = tmp_arr_s(np_s:1:-1) call lag3interp(Rc,theta,np,tmp_arr_s,theta_s_reverse,np_s) Rc_s = tmp_arr_s(np_s:1:-1) else call lag3interp(nu1,theta,np,nu1_s,theta_s,np_s) call lag3interp(Bp,theta,np,Bp_s,theta_s,np_s) call lag3interp(dlp,theta,np,dlp_s,theta_s,np_s) call lag3interp(Rc,theta,np,Rc_s,theta_s,np_s) endif psi1_s = R_s*Bp_s*sign_Ip_CW dBp_dchi_s=deriv_fd(Bp_s,chi_s,np_s) Bphi_s=F/R_s B_s=sqrt(Bphi_s**2+Bp_s**2) cosu_s=Z_theta_s/dlp_s/dtheta_dchi_s sinu_s=-R_theta_s/dlp_s/dtheta_dchi_s !radial derivative of straight field line angle dchidx_s=-(nu1_s/psi1_s*dxPsi+chi_s*dq_dx)/q0 !Bfield derivatives in Mercier-Luc coordinates (varrho,l,fz) dB_drho_s=-1./B_s*(F**2/R_s**3*cosu_s+Bp_s**2/Rc_s+mu_0*psi1_s*pprime) dB_dl_s=1./B_s*(Bp_s*dBp_dchi_s/dtheta_dchi_s/dlp_s+F**2/R_s**3*sinu_s) dnu_drho_s=nu1_s dnu_dl_s=-F/(R_s*psi1_s) grad_nu_s=sqrt(dnu_drho_s**2+dnu_dl_s**2) !contravariant metric coefficients (varrho,l,fz)->(x,y,z) gxx=(psi1_s/dxPsi)**2 gxy=-psi1_s/dxPsi*C_y*sign_Ip_CW*nu1_s gxz=-psi1_s/dxPsi*(nu1_s+psi1_s*dq_dpsi*chi_s)/q0 gyy=C_y**2*(grad_nu_s**2+1/R_s**2) gyz=sign_Ip_CW*C_y/q0*(grad_nu_s**2+dq_dpsi*nu1_s*psi1_s*chi_s) gzz=1./q0**2*(grad_nu_s**2+2.*dq_dpsi*nu1_s*psi1_s*chi_s+(dq_dpsi*psi1_s*chi_s)**2) jacobian=1./sqrt(gxx*gyy*gzz + 2.*gxy*gyz*gxz - gxz**2*gyy - gyz**2*gxx - gzz*gxy**2) !covariant metric coefficients g_xx=jacobian**2*(gyy*gzz-gyz**2) g_xy=jacobian**2*(gxz*gyz-gxy*gzz) g_xz=jacobian**2*(gxy*gyz-gxz*gyy) g_yy=jacobian**2*(gxx*gzz-gxz**2) g_yz=jacobian**2*(gxz*gxy-gxx*gyz) g_zz=jacobian**2*(gxx*gyy-gxy**2) !Bfield derivatives !dBdx = e_x * nabla B = J (nabla y x nabla z) * nabla B dBdx=jacobian*C_y/(q0*R_s)*(F/(R_s*psi1_s)*dB_drho_s+(nu1_s+dq_dpsi*chi_s*psi1_s)*dB_dl_s) dBdz=1./B_s*(Bp_s*dBp_dchi_s-F**2/R_s**3*R_theta_s) !curvature terms (these are just local and will be recalculated in geometry.F90) K_x = (0.-g_yz/g_zz*dBdz) K_y = (dBdx-g_xz/g_zz*dBdz) !(R,Z) derivatives for visualization dxdR_s = dx_drho/drPsi*psi1_s*cosu_s dxdZ_s = dx_drho/drPsi*psi1_s*sinu_s if (edge_opt==0.0) then !gene z-grid chi_out=linspace(-pi*Npol,pi*Npol-2*pi*Npol/Nz,Nz) else !new parallel coordinate chi_out==zprime !see also tracer_aux.F90 if (Npol>1) STOP "ERROR: Npol>1 has not been implemented for edge_opt=\=0.0" do k=izs,ize chi_out(k)=sinh((-pi+k*2.*pi/Nz)*log(edge_opt*pi+sqrt(edge_opt**2*pi**2+1))/pi)/edge_opt enddo !transform metrics according to chain rule do k=1,np_s !>dz'/dz conversion for edge_opt as function of z if (edge_opt.gt.0) then dzprimedz = edge_opt*pi/log(edge_opt*pi+sqrt((edge_opt*pi)**2+1))/& sqrt((edge_opt*chi_s(k))**2+1) else dzprimedz = 1.0 endif gxz(k)=gxz(k)*dzprimedz gyz(k)=gyz(k)*dzprimedz gzz(k)=gzz(k)*dzprimedz**2 jacobian(k)=jacobian(k)/dzprimedz dBdz(k)=dBdz(k)/dzprimedz enddo endif !edge_opt !interpolate down to GENE z-grid call lag3interp(gxx,chi_s,np_s,gxx_out,chi_out,Nz) call lag3interp(gxy,chi_s,np_s,gxy_out,chi_out,Nz) call lag3interp(gxz,chi_s,np_s,gxz_out,chi_out,Nz) call lag3interp(gyy,chi_s,np_s,gyy_out,chi_out,Nz) call lag3interp(gyz,chi_s,np_s,gyz_out,chi_out,Nz) call lag3interp(gzz,chi_s,np_s,gzz_out,chi_out,Nz) call lag3interp(B_s,chi_s,np_s,Bfield_out,chi_out,Nz) call lag3interp(jacobian,chi_s,np_s,jacobian_out,chi_out,Nz) call lag3interp(dBdx,chi_s,np_s,dBdx_out,chi_out,Nz) call lag3interp(dBdz,chi_s,np_s,dBdz_out,chi_out,Nz) call lag3interp(R_s,chi_s,np_s,R_out,chi_out,Nz) call lag3interp(Z_s,chi_s,np_s,Z_out,chi_out,Nz) call lag3interp(dxdR_s,chi_s,np_s,dxdR_out,chi_out,Nz) call lag3interp(dxdZ_s,chi_s,np_s,dxdZ_out,chi_out,Nz) ! Fill the geom arrays with the results do eo=0,1 gxx_(izs:ize,eo) =gxx_out(izs:ize) gyy_(izs:ize,eo) =gyy_out(izs:ize) gxz_(izs:ize,eo) =gxz_out(izs:ize) gyz_(izs:ize,eo) =gyz_out(izs:ize) dBdx_(izs:ize,eo) =dBdx_out(izs:ize) dBdy_(izs:ize,eo) =0. gxy_(izs:ize,eo) =gxy_out(izs:ize) gzz_(izs:ize,eo) =gzz_out(izs:ize) Bfield_(izs:ize,eo) =Bfield_out(izs:ize) jacobian_(izs:ize,eo) =jacobian_out(izs:ize) dBdz_(izs:ize,eo) =dBdz_out(izs:ize) R_hat_(izs:ize,eo) =R_out(izs:ize) Z_hat_(izs:ize,eo) =Z_out(izs:ize) dxdR_(izs:ize,eo) = dxdR_out(izs:ize) dxdZ_(izs:ize,eo) = dxdZ_out(izs:ize) !! UPDATE GHOSTS VALUES (since the miller function in GENE does not) CALL update_ghosts_z(gxx_(:,eo)) CALL update_ghosts_z(gyy_(:,eo)) CALL update_ghosts_z(gxz_(:,eo)) CALL update_ghosts_z(gxy_(:,eo)) CALL update_ghosts_z(gzz_(:,eo)) CALL update_ghosts_z(Bfield_(:,eo)) CALL update_ghosts_z(dBdx_(:,eo)) CALL update_ghosts_z(dBdy_(:,eo)) CALL update_ghosts_z(dBdz_(:,eo)) CALL update_ghosts_z(jacobian_(:,eo)) CALL update_ghosts_z(R_hat_(:,eo)) CALL update_ghosts_z(Z_hat_(:,eo)) CALL update_ghosts_z(dxdR_(:,eo)) CALL update_ghosts_z(dxdZ_(:,eo)) ! Curvature operator (Frei et al. 2022 eq 2.15) DO iz = izgs,izge - G1 = gxx_(iz,eo)*gyy_(iz,eo)-gxy_(iz,eo)*gxy_(iz,eo) - G2 = gxx_(iz,eo)*gyz_(iz,eo)-gxy_(iz,eo)*gxz_(iz,eo) - G3 = gxy_(iz,eo)*gyz_(iz,eo)-gyy_(iz,eo)*gxz_(iz,eo) + 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*dBdy_(iz,eo) + G2*dBdz_(iz,eo))/Bfield_(iz,eo) Cy = (G3*dBdz_(iz,eo) - G1*dBdx_(iz,eo))/Bfield_(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)/Bfield_(iz,eo) + Ckxky_(iky, ikx, iz,eo) = (Cx*kx + Cy*ky) ENDDO ENDDO + ! coefficient in the front of parallel derivative + gradz_coeff_(iz,eo) = 1._dp / jacobian_(iz,eo) / Bfield_(iz,eo) ENDDO - ! coefficient in the front of parallel derivative - gradz_coeff_(iz,eo) = 1._dp / Jacobian_(iz,eo) / Bfield_(iz,eo) - enddo + ENDDO contains SUBROUTINE update_ghosts_z(fz_) IMPLICIT NONE ! INTEGER, INTENT(IN) :: nztot_ REAL(dp), DIMENSION(izgs:izge), INTENT(INOUT) :: fz_ REAL(dp), DIMENSION(-2:2) :: buff INTEGER :: status(MPI_STATUS_SIZE), count IF(Nz .GT. 1) THEN IF (num_procs_z .GT. 1) THEN CALL MPI_BARRIER(MPI_COMM_WORLD,ierr) count = 1 ! one point to exchange !!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!! CALL mpi_sendrecv(fz_(ize), count, MPI_DOUBLE, nbr_U, 0, & ! Send to Up the last buff(-1), count, MPI_DOUBLE, nbr_D, 0, & ! Receive from Down the first-1 comm0, status, ierr) CALL mpi_sendrecv(fz_(ize-1), count, MPI_DOUBLE, nbr_U, 0, & ! Send to Up the last buff(-2), count, MPI_DOUBLE, nbr_D, 0, & ! Receive from Down the first-2 comm0, status, ierr) !!!!!!!!!!! Send ghost to down neighbour !!!!!!!!!!!!!!!!!!!!!! CALL mpi_sendrecv(fz_(izs), count, MPI_DOUBLE, nbr_D, 0, & ! Send to Down the first buff(+1), count, MPI_DOUBLE, nbr_U, 0, & ! Recieve from Up the last+1 comm0, status, ierr) CALL mpi_sendrecv(fz_(izs+1), count, MPI_DOUBLE, nbr_D, 0, & ! Send to Down the first buff(+2), count, MPI_DOUBLE, nbr_U, 0, & ! Recieve from Up the last+2 comm0, status, ierr) ELSE buff(-1) = fz_(ize ) buff(-2) = fz_(ize-1) buff(+1) = fz_(izs ) buff(+2) = fz_(izs+1) ENDIF fz_(ize+1) = buff(+1) fz_(ize+2) = buff(+2) fz_(izs-1) = buff(-1) fz_(izs-2) = buff(-2) ENDIF END SUBROUTINE update_ghosts_z !> Generate an equidistant array from min to max with n points function linspace(min,max,n) result(out) real(dp):: min, max integer:: n real(dp), dimension(n):: out do i=1,n out(i)=min+(i-1)*(max-min)/(n-1) enddo end function linspace !> Weighted average real(dp) function average(var,weight) real(dp), dimension(np):: var, weight average=sum(var*weight)/sum(weight) end function average !> full theta integral with weight function dlp real(dp) function dlp_int(var,dlp) real(dp), dimension(np):: var, dlp dlp_int=sum(var*dlp)*2*pi*Npol_ext/np end function dlp_int !> theta integral with weight function dlp, up to index 'ind' real(dp) function dlp_int_ind(var,dlp,ind) real(dp), dimension(np):: var, dlp integer:: ind dlp_int_ind=0. if (ind.gt.1) then dlp_int_ind=dlp_int_ind+var(1)*dlp(1)*pi*Npol_ext/np dlp_int_ind=dlp_int_ind+(sum(var(2:ind-1)*dlp(2:ind-1)))*2*pi*Npol_ext/np dlp_int_ind=dlp_int_ind+var(ind)*dlp(ind)*pi*Npol_ext/np endif end function dlp_int_ind !> 1st derivative with 2nd order finite differences function deriv_fd(y,x,n) result(out) integer, intent(in) :: n real(dp), dimension(n):: x,y,out,dx !call lag3deriv(y,x,n,out,x,n) out=0. do i=2,n-1 out(i)=out(i)-y(i-1)/2 out(i)=out(i)+y(i+1)/2 enddo out(1)=y(2)-y(1) out(n)=y(n)-y(n-1) dx=x(2)-x(1) out=out/dx end function deriv_fd end subroutine get_miller END MODULE miller diff --git a/testcases/miller_example/bug/fort_01.90 b/testcases/miller_example/fort_01.90 similarity index 98% rename from testcases/miller_example/bug/fort_01.90 rename to testcases/miller_example/fort_01.90 index 7deb16c..9e3c980 100644 --- a/testcases/miller_example/bug/fort_01.90 +++ b/testcases/miller_example/fort_01.90 @@ -1,86 +1,86 @@ &BASIC nrun = 100000000 dt = 0.01 - tmax = 400 + tmax = 100 maxruntime = 356400 / &GRID pmaxe = 2 jmaxe = 1 pmaxi = 2 jmaxi = 1 Nx = 128 Lx = 100 Ny = 64 Ly = 120 Nz = 16 Npol = 1 SG = .false. / &GEOMETRY geom = 'miller' q0 = 1.4 shear = 0.0 eps = 0.18 / &OUTPUT_PAR nsave_0d = 50 nsave_1d = -1 nsave_2d = Inf nsave_3d = 100 nsave_5d = 500 write_doubleprecision = .false. write_gamma = .true. write_hf = .true. write_phi = .true. write_Na00 = .true. write_Napj = .true. write_Sapj = .false. write_dens = .true. write_temp = .true. job2load = 0 / &MODEL_PAR ! Collisionality CLOS = 0 NL_CLOS = -1 LINEARITY = 'nonlinear' KIN_E = .false. mu_x = 0.1 mu_y = 0.1 N_HD = 4 mu_z = 0.2 mu_p = 0 mu_j = 0 nu = 0.2 tau_e = 1 tau_i = 1 sigma_e = 0.023338 sigma_i = 1 q_e = -1 q_i = 1 K_Ni = 2.22 K_Ti = 6.92 K_Ne = 0 K_Te = 0 GradB = 1 CurvB = 1 lambdaD = 0 beta = 0 / &COLLISION_PAR collision_model = 'DG' gyrokin_CO = .false. interspecies = .true. mat_file = 'null' / &INITIAL_CON INIT_OPT = 'ppj' ACT_ON_MODES = 'none' init_background = 0 init_noiselvl = 1e-3 iseed = 42 / &TIME_INTEGRATION_PAR numerical_scheme = 'RK4' / diff --git a/testcases/miller_example/bug/fort_00.90 b/testcases/miller_example/fort_02.90 similarity index 96% rename from testcases/miller_example/bug/fort_00.90 rename to testcases/miller_example/fort_02.90 index 728a345..e0af42a 100644 --- a/testcases/miller_example/bug/fort_00.90 +++ b/testcases/miller_example/fort_02.90 @@ -1,86 +1,86 @@ &BASIC nrun = 100000000 - dt = 0.01 + dt = 0.0075 tmax = 100 maxruntime = 356400 / &GRID pmaxe = 2 jmaxe = 1 pmaxi = 2 jmaxi = 1 Nx = 128 Lx = 100 Ny = 64 Ly = 120 Nz = 16 Npol = 1 SG = .false. / &GEOMETRY geom = 'miller' q0 = 1.4 shear = 0.0 eps = 0.18 / &OUTPUT_PAR nsave_0d = 50 nsave_1d = -1 nsave_2d = Inf nsave_3d = 100 nsave_5d = 500 write_doubleprecision = .false. write_gamma = .true. write_hf = .true. write_phi = .true. write_Na00 = .true. write_Napj = .true. write_Sapj = .false. write_dens = .true. write_temp = .true. - job2load = -1 + job2load = 1 / &MODEL_PAR ! Collisionality CLOS = 0 NL_CLOS = -1 LINEARITY = 'nonlinear' KIN_E = .false. mu_x = 0.1 mu_y = 0.1 N_HD = 4 - mu_z = 0.2 + mu_z = 1.0 mu_p = 0 mu_j = 0 nu = 0.2 tau_e = 1 tau_i = 1 sigma_e = 0.023338 sigma_i = 1 q_e = -1 q_i = 1 K_Ni = 2.22 K_Ti = 6.92 K_Ne = 0 K_Te = 0 GradB = 1 CurvB = 1 lambdaD = 0 beta = 0 / &COLLISION_PAR collision_model = 'DG' gyrokin_CO = .false. interspecies = .true. mat_file = 'null' / &INITIAL_CON INIT_OPT = 'ppj' ACT_ON_MODES = 'none' init_background = 0 init_noiselvl = 1e-3 iseed = 42 / &TIME_INTEGRATION_PAR numerical_scheme = 'RK4' / diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m index 65ba29b..54353a9 100644 --- a/wk/analysis_gyacomo.m +++ b/wk/analysis_gyacomo.m @@ -1,228 +1,230 @@ %% UNCOMMENT FOR TUTORIAL % gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory % resdir = '.'; %Name of the directory where the results are located % JOBNUMMIN = 00; JOBNUMMAX = 10; %% addpath(genpath([gyacomodir,'matlab'])) % ... add addpath(genpath([gyacomodir,'matlab/plot'])) % ... add addpath(genpath([gyacomodir,'matlab/compute'])) % ... add addpath(genpath([gyacomodir,'matlab/load'])) % ... add %% Load the results LOCALDIR = [gyacomodir,resdir,'/']; MISCDIR = ['/misc/gyacomo_outputs/',resdir,'/']; %For long term storage system(['mkdir -p ',MISCDIR]); system(['mkdir -p ',LOCALDIR]); CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD); system(CMD); % Load outputs from jobnummin up to jobnummax data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing data.localdir = LOCALDIR; data.FIGDIR = LOCALDIR; %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 1 %% Space time diagramm (fig 11 Ivanov 2020) % data.scale = 1;%/(data.Nx*data.Ny)^2; -i_ = 1; +i_ = 3; disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))]) disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))]) options.TAVG_0 = data.TJOB_SE(i_);%0.4*data.Ts3D(end); options.TAVG_1 = data.TJOB_SE(i_+1);%0.9*data.Ts3D(end); % Averaging times duration options.NCUT = 4; % Number of cuts for averaging and error estimation options.NMVA = 1; % Moving average for time traces % options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) options.INTERP = 0; options.RESOLUTION = 256; fig = plot_radial_transport_and_spacetime(data,options); save_figure(data,fig,'.png') end if 0 %% statistical transport averaging for i_ = 1:2:21 % i_ = 3; disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))]) disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))]) options.T = [data.TJOB_SE(i_) data.TJOB_SE(i_+1)]; options.NPLOTS = 0; fig = statistical_transport_averaging(data,options); end end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options options.INTERP = 1; options.POLARPLOT = 0; -options.NAME = '\phi'; -% options.NAME = '\omega_z'; -% options.NAME = 'N_e^{00}'; +% options.NAME = '\phi'; +options.NAME = '\omega_z'; +% options.NAME = 'N_i^{00}'; % options.NAME = 'v_y'; % options.NAME = 'n_i^{NZ}'; % options.NAME = '\Gamma_x'; % options.NAME = 'n_i'; -options.PLAN = 'xy'; +options.PLAN = 'kxky'; % options.NAME = 'f_i'; % options.PLAN = 'sx'; options.COMP = 'avg'; % options.TIME = data.Ts5D(end-30:end); % options.TIME = data.Ts3D; options.TIME = [000:0.1:7000]; data.EPS = 0.1; data.a = data.EPS * 2000; options.RESOLUTION = 256; create_film(data,options,'.gif') end if 1 %% 2D snapshots % Options options.INTERP = 1; options.POLARPLOT = 0; options.AXISEQUAL = 1; options.NAME = '\phi'; % options.NAME = '\psi'; % options.NAME = 'n_i'; % options.NAME = 'N_i^{00}'; % options.NAME = 'T_i'; % options.NAME = '\Gamma_x'; % options.NAME = 'k^2n_e'; options.PLAN = 'xy'; % options.NAME 'f_i'; % options.PLAN = 'sx'; options.COMP = 'avg'; options.TIME = [1000 1800 2500 3000 4000]; data.a = data.EPS * 2e3; fig = photomaton(data,options); % save_figure(data,fig) end if 0 %% 3D plot on the geometry options.INTERP = 0; options.NAME = '\phi'; options.PLANES = [1]; options.TIME = [30]; options.PLT_MTOPO = 1; options.PLT_FTUBE = 0; data.EPS = 0.4; data.rho_o_R = 3e-3; % Sound larmor radius over Machine size ratio fig = show_geometry(data,options); save_figure(data,fig,'.png') end if 0 %% Kinetic distribution function sqrt(xy) (GENE vsp) options.SPAR = linspace(-3,3,32)+(6/127/2); options.XPERP = linspace( 0,6,32); % options.SPAR = gene_data.vp'; % options.XPERP = gene_data.mu'; options.iz = 'avg'; options.T = [250 600]; options.PLT_FCT = 'pcolor'; options.ONED = 0; options.non_adiab = 0; options.SPECIE = 'i'; options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene fig = plot_fa(data,options); % save_figure(data,fig,'.png') end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 0; options.ST = 1; options.PLOT_TYPE = 'space-time'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [1000]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); % save_figure(data,fig,'.png'); end if 0 %% Time averaged spectrum options.TIME = [2000 3000]; options.NORM =1; options.NAME = '\phi'; % options.NAME = 'N_i^{00}'; % options.NAME ='\Gamma_x'; options.PLAN = 'kxky'; options.COMPZ = 'avg'; options.OK = 0; options.COMPXY = 'avg'; % avg/sum/max/zero/ 2D plot otherwise options.COMPT = 'avg'; options.PLOT = 'semilogy'; fig = spectrum_1D(data,options); % save_figure(data,fig,'.png') end if 0 %% 1D real plot options.TIME = [50 100 200]; options.NORM = 0; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME ='\Gamma_x'; % options.NAME ='s_y'; options.COMPX = 'avg'; options.COMPY = 'avg'; options.COMPZ = 1; options.COMPT = 1; options.MOVMT = 1; fig = real_plot_1D(data,options); % save_figure(data,fig,'.png') end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = [0.1 0.2 0.3 0.4]; options.TIME = [00:1200]; options.NMA = 1; options.NMODES = 5; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end if 0 %% ZF caracteristics (space time diagrams) TAVG_0 = 1200; TAVG_1 = 1500; % Averaging times duration % chose your field to plot in spacetime diag (uzf,szf,Gx) fig = ZF_spacetime(data,TAVG_0,TAVG_1); save_figure(data,fig,'.png') end if 0 %% Metric infos -fig = plot_metric(data); +options.SHOW_FLUXSURF = 1; +options.SHOW_METRICS = 0; +fig = plot_metric(data,options); end if 0 %% linear growth rate for 3D fluxtube trange = [0 100]; nplots = 1; lg = compute_fluxtube_growth_rate(data,trange,nplots); end if 0 %% linear growth rate for 3D Zpinch trange = [5 15]; options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 1; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig,'.png') end diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m index 189cc0e..779dd11 100644 --- a/wk/header_3D_results.m +++ b/wk/header_3D_results.m @@ -1,64 +1,64 @@ % Directory of the code "mypathtoHeLaZ/HeLaZ/" gyacomodir = '/home/ahoffman/gyacomo/'; % Directory of the simulation (from results) % if 1% Local results % resdir ='volcokas/64x32x16x5x3_kin_e_npol_1'; %% Dimits % resdir ='shearless_cyclone/128x64x16x5x3_Dim_90'; % resdir ='shearless_cyclone/128x64x16x9x5_Dim_scan/128x64x16x9x5_Dim_60'; % resdir ='shearless_cyclone/128x64x16x5x3_Dim_scan/128x64x16x5x3_Dim_70'; % resdir ='shearless_cyclone/64x32x16x5x3_Dim_scan/64x32x16x5x3_Dim_70'; %% AVS % resdir = 'volcokas/64x32x16x5x3_kin_e_npol_1'; % resdir = 'volcokas/64x32x16x5x3_kin_e_npol_1'; % resdir = 'shearless_cyclone/64x32x80x5x3_CBC_Npol_5_kine'; % resdir = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; % resdir = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine'; % resdir = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; %% shearless CBC % resdir ='shearless_cyclone/64x32x16x5x3_CBC_080'; % resdir ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100'; % resdir ='shearless_cyclone/64x32x16x5x3_CBC_120'; % resdir ='shearless_cyclone/64x32x16x9x5_CBC_080'; % resdir ='shearless_cyclone/64x32x16x9x5_CBC_100'; % resdir ='shearless_cyclone/64x32x16x9x5_CBC_120'; % resdir = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK'; %% CBC % resdir = 'CBC/64x32x16x5x3'; % resdir = 'CBC/64x128x16x5x3'; % resdir = 'CBC/128x64x16x5x3'; % resdir = 'CBC/96x96x16x3x2_Nexc_6'; % resdir = 'CBC/128x96x16x3x2'; % resdir = 'CBC/192x96x16x3x2'; % resdir = 'CBC/192x96x24x13x7'; % resdir = 'CBC/kT_11_128x64x16x5x3'; % resdir = 'CBC/kT_9_256x128x16x3x2'; % resdir = 'CBC/kT_4.5_128x64x16x13x3'; % resdir = 'CBC/kT_4.5_192x96x24x13x7'; % resdir = 'CBC/kT_4.5_128x64x16x13x7'; % resdir = 'CBC/kT_4.5_128x96x24x15x5'; % resdir = 'CBC/kT_5.3_192x96x24x13x7'; % resdir = 'CBC/kT_13_large_box_128x64x16x5x3'; % resdir = 'CBC/kT_11_96x64x16x5x3_ky_0.02'; % resdir = 'CBC/kT_scan_128x64x16x5x3'; % resdir = 'CBC/kT_scan_192x96x16x3x2'; % resdir = 'CBC/kT_13_96x96x16x3x2_Nexc_6'; % resdir = 'dbg/nexc_dbg'; % resdir = 'CBC/NM_F4_kT_4.5_192x64x24x6x4'; % resdir = 'CBC_Ke_EM/192x96x24x5x3'; % resdir = 'CBC_Ke_EM/96x48x16x5x3'; % resdir = 'CBC_Ke_EM/minimal_res'; %% KBM % resdir = 'NL_KBM/192x64x24x5x3'; %% Linear CBC % resdir = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe'; -resdir = 'testcases/miller_example_bug'; +resdir = 'testcases/miller_example'; JOBNUMMIN = 00; JOBNUMMAX = 10; run analysis_gyacomo diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m index cffb9cc..ab514b2 100644 --- a/wk/lin_ITG.m +++ b/wk/lin_ITG.m @@ -1,184 +1,189 @@ %% QUICK RUN SCRIPT % This script create a directory in /results and run a simulation directly % from matlab framework. It is meant to run only small problems in linear % for benchmark and debugging purpose since it makes matlab "busy" % SIMID = 'lin_ITG'; % Name of the simulation RUN = 1; % To run or just to load addpath(genpath('../matlab')) % ... add default_plots_options HELAZDIR = '/home/ahoffman/gyacomo/'; +% EXECNAME = 'gyacomo_1.0'; +% EXECNAME = 'gyacomo_dbg'; EXECNAME = 'gyacomo'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS -NU = 0.0; % Collision frequency +NU = 0.05; % Collision frequency TAU = 1.0; % e/i temperature ratio K_Ne = 2.22; % ele Density ''' K_Te = 6.96; % ele Temperature ''' K_Ni = 2.22; % ion Density gradient drive K_Ti = 6.96; % ion Temperature ''' SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) % SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons BETA = 0.0; % electron plasma beta %% GRID PARAMETERS P = 4; J = P/2; PMAXE = P; % Hermite basis size of electrons JMAXE = J; % Laguerre " PMAXI = P; % " ions JMAXI = J; % " -NX = 11; % real space x-gridpoints -NY = 32; % '' y-gridpoints +NX = 20; % real space x-gridpoints +NY = 2; % '' y-gridpoints LX = 2*pi/0.8; % Size of the squared frequency domain -LY = 2*pi/0.05; % Size of the squared frequency domain -NZ = 16; % number of perpendicular planes (parallel grid) +LY = 2*pi/0.3; % Size of the squared frequency domain +NZ = 32; % number of perpendicular planes (parallel grid) NPOL = 1; SG = 0; % Staggered z grids option %% GEOMETRY -GEOMETRY= 's-alpha'; -% GEOMETRY= 'miller'; +% GEOMETRY= 's-alpha'; +GEOMETRY= 'miller'; Q0 = 1.4; % safety factor SHEAR = 0.8; % magnetic shear +KAPPA = 0.0; % elongation +DELTA = 0.0; % triangularity +ZETA = 0.0; % squareness NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) EPS = 0.18; % inverse aspect ratio %% TIME PARMETERS -TMAX = 50; % Maximal time unit -DT = 1e-2; % Time step +TMAX = 1; % Maximal time unit +DT = 3e-3; % Time step SPS0D = 1; % Sampling per time unit for 2D arrays SPS2D = 0; % Sampling per time unit for 2D arrays SPS3D = 5; % Sampling per time unit for 2D arrays SPS5D = 1/5; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints JOB2LOAD= -1; %% OPTIONS LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) % Collision operator % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) CO = 'DG'; GKCO = 0; % gyrokinetic operator ABCO = 1; % interspecies collisions INIT_ZF = 0; ZF_AMP = 0.0; CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS) KERN = 0; % Kernel model (0 : GK) INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom %% OUTPUTS W_DOUBLE = 1; W_GAMMA = 1; W_HF = 1; W_PHI = 1; W_NA00 = 1; W_DENS = 1; W_TEMP = 1; W_NAPJ = 1; W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % unused HD_CO = 0.0; % Hyper diffusivity cutoff ratio MU = 0.0; % Hyperdiffusivity coefficient INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; MU_X = MU; % MU_Y = MU; % N_HD = 4; MU_Z = 2.0; % MU_P = 0.0; % MU_J = 0.0; % LAMBDAD = 0.0; NOISE0 = 0.0e-5; % Init noise amplitude BCKGD0 = 1.0; % Init background GRADB = 1.0; CURVB = 1.0; %%------------------------------------------------------------------------- %% RUN setup % system(['rm fort*.90']); % Run linear simulation if RUN % system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk']) system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk']) end %% Load results %% filename = [SIMID,'/',PARAMS,'/']; LOCALDIR = [HELAZDIR,'results/',filename,'/']; % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 00; data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing %% Short analysis if 1 %% linear growth rate (adapted for 2D zpinch and fluxtube) trange = [0.5 1]*data.Ts3D(end); nplots = 3; lg = compute_fluxtube_growth_rate(data,trange,nplots); [gmax, kmax] = max(lg.g_ky(:,end)); [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky); msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg); msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg); end if 1 %% Ballooning plot options.time_2_plot = [120]; -options.kymodes = [0.9]; +options.kymodes = [0.3]; options.normalized = 1; % options.field = 'phi'; fig = plot_ballooning(data,options); end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 1; options.ST = 1; options.PLOT_TYPE = 'space-time'; % options.PLOT_TYPE = 'Tavg-1D'; % options.PLOT_TYPE = 'Tavg-2D'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [0 50]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); save_figure(data,fig) end if 0 %% linear growth rate for 3D Zpinch (kz fourier transform) trange = [0.5 1]*data.Ts3D(end); options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 0; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig) end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = 1; options.TIME = [0:1000]; options.NMA = 1; options.NMODES = 1; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end if 1 %% RH TEST ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end); [~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D)); plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3)); figure plot(data.Ts3D(it0:it1), plt(data.PHI)); xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$') title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky))) end