Diffusion GYACOMO (OBSOLETE: LATEST VERSION ON gitlab.epfl.ch) (master) Browse

CBC_kT_PJ_scan.m
5bb85d4b05f5master

CBC_kT_PJ_scan.m
View Options

			addpath(genpath('../matlab')) % ... add
			default_plots_options
			HELAZDIR = '/home/ahoffman/HeLaZ/';
			EXECNAME = 'helaz3';
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			KT_a = [3:0.5:5];
			P_a_6 = [6 6 6 6 6 6 6];
			J_a_6 = [0 1 2 3 4 5 6];
			P_a_10 = 10*ones(1,6);
			J_a_10 = 5:10;
			% P_a_10b = 10*ones(1,7);
			% J_a_10b = 10:17;
			P_a_20 = 20*ones(1,11);
			J_a_20 = 10:20;

			P_a = [P_a_20]; J_a = [J_a_20];
			% KT_a = 5.0; P_a = 20; J_a = 20;

			g_max= zeros(numel(P_a),numel(KT_a));
			g_avg= g_max*0;
			g_std= g_max*0;
			k_max= g_max*0;

			CO = 'DG'; GKCO = 0;
			NU = 0.0;
			DT = 7e-3;
			TMAX = 40;
			ky_ = 0.15;
			SIMID = 'linear_CBC_kT_threshold'; % Name of the simulation
			RUN = 0;


			for i = 1:numel(P_a)
			P = P_a(i); J = J_a(i);
			j=1;
			%Set Up parameters
			for K_Ti = KT_a
			CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
			TAU = 1.0; % e/i temperature ratio
			K_Ni = 2.22; K_Ne = K_Ni;
			K_Te = K_Ti; % Temperature '''
			SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
			KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
			BETA = 0e-1; % electron plasma beta
			PMAXE = P; JMAXE = J;
			PMAXI = P; JMAXI = J;
			NX = 8; % real space x-gridpoints
			NY = 6; % '' y-gridpoints
			LX = 2*pi/0.15; % Size of the squared frequency domain
			LY = 2*pi/ky_;
			NZ = 24; % number of perpendicular planes (parallel grid)
			NPOL = 1; SG = 0; NEXC = 1;
			GEOMETRY= 's-alpha';
			Q0 = 1.4; % safety factor
			SHEAR = 0.8; % magnetic shear (Not implemented yet)
			EPS = 0.18; % inverse aspect ratio
			SPS0D = 1; SPS2D = 0; SPS3D = 5;SPS5D= 1/5; SPSCP = 0;
			JOB2LOAD= -1;
			LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
			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= 'phi'; % Start simulation with a noisy mom00/phi/allmom
			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;
			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 = 1.0; MU_P = 0.0; %
			MU_J = 0.0; LAMBDAD = 0.0;
			NOISE0 = 1.0e-4; % Init noise amplitude
			BCKGD0 = 0.0; % Init background
			GRADB = 1.0;CURVB = 1.0;

			%%-------------------------------------------------------------------------
			% RUN
			setup
			if RUN
			system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 2 2 1 0; cd ../../../wk'])
			end

			% Load results
			filename = [SIMID,'/',PARAMS,'/'];
			LOCALDIR = [HELAZDIR,'results/',filename,'/'];
			data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing

			%linear growth rate (adapted for 2D zpinch and fluxtube)
			trange = [0.5 1]*data.Ts3D(end);
			nplots = 0;
			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);

			g_max(i,j) = gmax;
			k_max(i,j) = kmax;

			[g_avg(i,j), ik_] = max(lg.avg_g);
			g_std(i,j) = max(lg.std_g(ik_));

			j = j + 1;
			if 0
			%% Verify gamma time trace
			figure
			for ik_ = 1:numel(lg.ky)
			plot(lg.trange(2:end),lg.g_ky(ik_,:)','DisplayName',['$k_y=',num2str(lg.ky(ik_)),'$']); hold on;
			end
			xlabel('$t$'); ylabel('$\gamma$');
			title(data.param_title); legend('show');
			drawnow
			end
			end
			end

			if 1
			%% PLOTS
			ERR_WEIGHT = 1/3; %weight of the error to compute marginal stability
			%% Superposed 1D plots
			sz_ = size(g_max);
			figure
			for i = 1:sz_(1)
			y_ = g_avg(i,:);
			e_ = g_std(i,:);

			y_ = y_.(y_-e_ERR_WEIGHT>0);
			e_ = e_ .* (y_>0);
			errorbar(KT_a,y_,e_,...
			'LineWidth',1.2,...
			'DisplayName',['(',num2str(P_a(i)),',',num2str(J_a(i)),')']);
			% plot(KT_a,y_,...
			% 'LineWidth',1.2,...
			% 'DisplayName',['(',num2str(P_a(i)),',',num2str(J_a(i)),')']);
			hold on;

			end
			title('Linear CBC $K_T$ threshold');
			legend('show'); xlabel('$K_T$'); ylabel('$\max_{k_y}(\gamma_k)$');
			drawnow
			%% Color map
			[NP__, KT__] = meshgrid(P_a+2*J_a, KT_a);
			% GG_ = g_avg;
			GG_ = g_avg .* (g_avg-g_std*ERR_WEIGHT > 0);
			figure;
			% pclr = pcolor(KT__,NP__,g_max'); set(pclr,'EdgeColor','none');
			pclr = imagesc(KT_a,1:numel(P_a),GG_);
			LABELS = [];
			for i_ = 1:numel(P_a)
			LABELS = [LABELS; '(',sprintf('%2.0f',P_a(i_)),',',sprintf('%2.0f',J_a(i_)),')'];
			end
			yticks(1:numel(P_a));
			yticklabels(LABELS);
			xlabel('$\kappa_T$'); ylabel('$(P,J)$');
			title('Linear ITG threshold in CBC');
			colormap(bluewhitered);
			%%
			%%
			end

gyacomo/wk/CBC_kT_PJ_scan.m5bb85d4b05f5master

CBC_kT_PJ_scan.mView Options

gyacomo/wk/CBC_kT_PJ_scan.m
5bb85d4b05f5master

CBC_kT_PJ_scan.m
View Options