Homec4science
Diffusion FENNECS lite (3Dmpimode) Browse

fennecslite/matlab/@fennecslitehdf5/displaytotcurrevol_geom.m
6b37315661933Dmpimode

Last Change
Raw File

displaytotcurrevol_geom.m
View Options

function displaytotcurrevol_geom(obj,timesteps,toptitle,scalet,dens,subdiv,nmean, rsplit)
% Computes and display the time evolution of the outgoing
% currents at time obj.t2d(timesteps)
%-scalet=true scales the time by the ellastic collision
% frequency
%- dens = true plot the time evolution of the maximum electron
% density in the simulation domain otherwise plot the total
% number of electrons in the domain
% also plot in a subplot the color coded boundary corresponding
% to each current
% - subdiv allows the subdivision of individual boundaries by
% defining relative length e.g. [1 1] divides the boundary in
% two boundaries of equal length
% - nmean defines the window time average for the current and
% density
% - rsplit defines the radial position, for the calculation of the
% time evolution of the density, at which a new cloud is considered.
% This helps tracking the density in two clouds at the same
% time.
if nargin<2
timesteps=1:length(obj.t2d);
end
if nargin<3
toptitle="";
end
if nargin<4
scalet=1;
end
if nargin <5
dens=true;
end
if nargin<6
subdiv=1;
end
if nargin<7
nmean=1;
end
if nargin<8
rsplit=0.07;
%rsplit=0.014;
end
if scalet==1
if obj.neutcol.present
vexb0=(obj.Ezxt(:,:,1).*obj.Br'-obj.Erxt(:,:,1).*obj.Bz')./(obj.B'.^2);
vexb0(obj.geomweight(:,:,1)<=0)=0;
potwell=obj.PotentialWell(0);
vexb0(isnan(potwell))=NaN;
vexb0=mean(abs(vexb0(:)),'omitnan');
E=0.5*obj.msim/obj.weight*vexb0^2/obj.qe;
taucol=1/(obj.neutcol.neutdens*vexb0*(obj.sigio(E)+obj.sigmela(E)+obj.sigmio(E)))
try
Sio_S=1e17*(obj.neutcol.neutdens*vexb0*obj.sigio(E))/(obj.maxwellsrce.frequency*obj.weight/(pi*(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2)*diff(obj.maxwellsrce.zlim)));
catch
end
tlabel='t/\tau_d [-]';
else
taucol=2*pi/obj.omece;
tlabel='t/\tau_ce [-]';
end
elseif scalet==2
if obj.neutcol.present && obj.psupply.active
taucol=1/(obj.neutcol.neutdens/obj.psupply.expdens*1000);
tlabel='t [ms]';
else
taucol=2*pi/obj.omece;
tlabel='t/\tau_ce [-]';
end
else
taucol=1e-9;
tlabel='t [ns]';
end
if dens
N=obj.N(:,:,timesteps);
geomw=obj.geomweight(:,:,1);
geomw(geomw<0)=0;
geomw(geomw>0)=1;
N=N.*geomw;
%[~,idl]=max(N(:,:,end),[],'all','linear');
%[ir,iz]=ind2sub(size(geomw),idl);
%nmax=squeeze(max(max(N,[],1),[],2));
tn=(obj.t2d(timesteps));
nrhalf=find(obj.rgrid>rsplit,1,'first');
if(nrhalf<length(obj.rgrid))
nmax=zeros(2,length(tn));
nmax(1,:)=squeeze(max(max(N(1:nrhalf,:,:),[],1),[],2));
nmax(2,:)=squeeze(max(max(N(nrhalf+1:end,:,:),[],1),[],2));
else
nmax(:)=squeeze(max(max(N(:,:,:),[],1),[],2));
end
nlabel='n_{e,max} [m^{-3}]';
ndlabel='n_{e,max}';
else
t0dst=find(obj.t0d>=obj.t2d(timesteps(1)),1,'first');
t0dlst=find(obj.t0d<=obj.t2d(timesteps(end)),1,'last');
tn=obj.t0d(t0dst:t0dlst);
nmax=obj.npart(t0dst:t0dlst)*obj.weight;
nlabel='Nb e^-';
ndlabel='Nb e^-';
end
% We get the current at the desired time steps for each
% boundary
[currents,pos]=obj.OutCurrents(timesteps,subdiv);
% The currents are normalized to the simulated pressure
P=1;
if obj.neutcol.present
P=obj.neutcol.neutdens*obj.kb*300/100;% pressure at room temperature in mbar
end
currents=currents/P;
if scalet==2
currents=currents*obj.psupply.expdens*obj.kb*300/100*1000;
end
% Start of the figure
f=figure('Name',sprintf('%s Charges',obj.name));
tiledlayout(2,1)
nexttile
% Plot the evolution of nb of particles or density
yyaxis right
if size(nmax,2)>1
nmax=nmax';
end
for i=1:size(nmax,2)
p=plot(tn/taucol,nmax(:,i),'b','linewidth',1.8,'Displayname',sprintf('%s, %d',ndlabel,i));
hold on
end
ylabel(nlabel)
axl=gca;
axl.YAxis(2).Color=p.Color;
axl.YAxis(1).Color='k';
ylim([0 inf])
if(obj.B(1,1)>obj.B(end,1))
lname='HFS';
rname='LFS';
else
lname='LFS';
rname='HFS';
end
% Plot of the different current contributions
yyaxis( 'left');
map=colormap(lines);
set(axl,'linestyleorder',{'-',':','--','*','+'},...
'ColorOrder',map(2:7,:), 'NextPlot','replacechildren')
p(1)=plot(axl,obj.t2d(timesteps)/taucol,movmean(currents(1,:),nmean),'Displayname',lname,'linewidth',1.8);
hold on
p(2)=plot(axl,obj.t2d(timesteps)/taucol,movmean(currents(2,:),nmean),'Displayname',rname,'linewidth',1.8);
% Plot the currents
for i=3:size(currents,1)
p(i)=plot(axl,obj.t2d(timesteps)/taucol,movmean(currents(i,:),nmean),'Displayname',sprintf('border %i',i-2),'linewidth',1.8);
end
plot(axl,obj.t2d(timesteps)/taucol,movmean(sum(currents(:,:),1,'omitnan'),nmean),'k-','Displayname','total','linewidth',1.8);
xlabel(tlabel)
if obj.neutcol.present && scalet~=2
ylabel('I/p_n [A/mbar]')
else
ylabel('I [mA]')
end
grid on
set(gca,'fontsize',12)
ax.YAxis(1).Color='black';
%legend('Orientation','horizontal','location','north','numcolumns',3)
if ~isempty(toptitle)
title(toptitle)
end
% Plot of the geometry and of the different boundaries
% considered for the calculation of the currents
ax2=nexttile;
%% Magnetic field lines
Blines=obj.rAthet;
levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),15);
Blines(obj.geomweight(:,:,1)<0)=NaN;
[~,h1]=contour(ax2,obj.zgrid*1000,obj.rgrid*1000,Blines,real(levels),'-.','color','k','linewidth',1.1,'Displayname','Magnetic field lines');
hold on
%% Boundaries
geomw=obj.geomweight(:,:,1);
geomw(geomw<1e-6)=0;
geomw(geomw>1e-6)=NaN;
[c1,hContour]=contourf(ax2,obj.zgrid*1000,obj.rgrid*1000,geomw, [0 0]);
hold on
dirw=obj.dirichletweight(:,:,1);
geomw=obj.geomweight(:,:,1);
dirw(geomw<=1e-6&dirw>=1e-6)=0;
dirw(dirw>0)=NaN;
[c1,hContour2]=contourf(ax2,obj.zgrid*1000,obj.rgrid*1000,dirw, [0 0]);
drawnow;
grid on;
for i=1:length(pos)
plot(ax2,pos{i}(1,:)*1000,pos{i}(2,:)*1000,'linestyle',p(i+2).LineStyle,...
'color',p(i+2).Color,'marker',p(i+2).Marker,...
'displayname',sprintf('border %i',i),'linewidth',1.8)
hold on
end
title('Domain')
plot(ax2,ones(size(obj.rgrid))*obj.zgrid(1)*1000,obj.rgrid*1000,'linestyle',p(1).LineStyle,...
'color',p(1).Color,'marker',p(1).Marker,...
'displayname',lname,'linewidth',1.8)
plot(ax2,ones(size(obj.rgrid))*obj.zgrid(end)*1000,obj.rgrid*1000,'linestyle',p(2).LineStyle,...
'color',p(2).Color,'marker',p(2).Marker,...
'displayname',rname,'linewidth',1.8)
xlabel('z [mm]')
ylabel('r [mm]')
grid on
set(gca,'fontsize',12)
% We add the coloring for the solid parts of the geometry
% Metallic parts
hFills=hContour.FacePrims;
[hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor'
try
hFills(1).ColorData = uint8([150;150;150;255]);
for idx = 2 : numel(hFills)
hFills(idx).ColorData(4) = 0; % default=255
end
catch
end
% Dielectrics
hFills=hContour2.FacePrims;
[hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor'
try
hFills(1).ColorData = uint8([255;165;0;255]);
for idx = 2 : numel(hFills)
hFills(idx).ColorData(4) = 0; % default=255
end
catch
end
%legend('Orientation','horizontal','location','north','numcolumns',4)
% Output to console the mean total current at the end of the
% time considered
fprintf('mean total current: %f [A/mbar]\n',mean(sum(currents(:,max(1,size(currents,2)-50):end),1,'omitnan')));
%if nargin <3
% sgtitle(sprintf('\\phi_b-\\phi_a=%.2g kV, B=%f T',(obj.potout-obj.potinn)/1e3,mean(obj.B(:))))
%elseif ~isempty(toptitle)
% sgtitle(toptitle)
%end
xlim([obj.zgrid(1) obj.zgrid(end)]*1e3)
ylim([obj.rgrid(1) obj.rgrid(end)]*1e3)
% Save to file
if length(subdiv)>1
obj.savegraph(f,sprintf('%s/%s_totIEvol%i%i_div%i',obj.folder,obj.name,scalet,dens,nmean),[16 14]);
else
obj.savegraph(f,sprintf('%s/%s_totIEvol%i%i_%i',obj.folder,obj.name,scalet,dens,nmean),[16 14]);
end
end
c4science · Help