Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (master) Browse

fennecs/matlab/FBspline_article.m
dcbebf0c577cmaster

Last Change
Raw File

FBspline_article.m
View Options

function [Forces, Continuity]=FBspline_article(M,Forces, Continuity, Density,norm,log,parper,contourscale, zlims, rlims, clims)
%ForceBalance Show the radial force balance
% Plot the three radial forces for the given time-step it or averaged
% over the range of time steps defined in it
it=Forces.it;
switch nargin
case 4
parper=false;
contourscale=0.2;
norm=false;
log=false;
zlims=[-inf inf];
case 5
parper=false;
log=false;
contourscale=0.2;
zlims=[-inf inf];
case 6
parper=false;
contourscale=0.2;
zlims=[-inf inf];
case 7
contourscale=0.2;
zlims=[-inf inf];
case 8
zlims=[-inf inf];
case 9
rlims=[];
case 10
end
if nargin<10
rlims=[];
end
if nargin<11
climsr=[-230000 190000];
climsz=[-15000 13000];
climsthet=[-1100 1800];
else
climsr=clims(1,:);
climsthet=clims(2,:);
climsz=clims(3,:);
end
folder='Article_Forces';
[status, msg, msgID] = mkdir(folder);
contourcolor=[0 0 0];%[255 20 147]/255;
N=Density.N;
n=mean(N,3);
maxdens=max(N(:));
densitycontour=contourc(M.zgrid,M.rgrid,mean(N,3),[contourscale contourscale]*maxdens);
geometriccontour=contourc(M.zgrid,M.rgrid,M.geomweight(:,:,1),[0 0]);
Fmaxr=max([max(Forces.Eforcer(:)),max(Forces.Bforcer(:)),max(Forces.inertforcer(:)),max(Forces.Pressforcer(:))]);
Fmaxz=max([max(Forces.Eforcez(:)),max(Forces.Bforcez(:)),max(Forces.Pressforcez(:)),max(Forces.inertforcez(:))]);
FEr=mean(Forces.Eforcer,3);
FBr=mean(Forces.Bforcer,3);
FPr=mean(Forces.Pressforcer,3);
FIr=mean(Forces.inertforcer,3);
Fdr=mean(Forces.Dragforcer,3);
FBthet=mean(Forces.Bforcethet,3);
%FPthet=zeros(size(FBthet));%mean(Forces.Pressforcethet,3);
FPthet=mean(Forces.Pressforcethet,3);
FIthet=mean(Forces.inertforcethet,3);
Fdthet=mean(Forces.Dragforcethet,3);
FEz=mean(Forces.Eforcez,3);
FBz=mean(Forces.Bforcez,3);
FPz=mean(Forces.Pressforcez,3);
FIz=mean(Forces.inertforcez,3);
Fdz=mean(Forces.Dragforcez,3);
durdt=mean(Forces.durdt,3);
duthetdt=mean(Forces.duthetdt,3);
duzdt=mean(Forces.duzdt,3);
totforcer=(FEr+FBr+FPr+FIr+Fdr-durdt);
totforcethet=(FBthet+FPthet+FIthet+Fdthet-duthetdt);
totforcez=(FEz+FBz+FPz+FIz+Fdz-duzdt);
if log
Er=Forces.Eforcer(Forces.Eforcer~=0);
Br=Forces.Bforcer(Forces.Bforcer~=0);
inertr=Forces.inertforcer(Forces.inertforcer~=0);
pressr=Forces.Pressforcer(Forces.Pressforcer~=0);
Fminr=min([min(abs(Er(:))),min(abs(Br(:))),min(abs(inertr(:))),min(abs(pressr(:)))]);
else
Fminr=min([min(Forces.Eforcer(:)),min(Forces.Bforcer(:)),min(Forces.inertforcer(:)),min(Forces.Pressforcer(:))]);
end
rgridmax=sum(M.nnr(1:2));
if parper
costhet=(M.Br./M.B)';
sinthet=(M.Bz./M.B)';
FErd=FEr.*sinthet-FEz.*costhet;
FBrd=FBr.*sinthet-FBz.*costhet;
FPrd=FPr.*sinthet-FPz.*costhet;
FIrd=FIr.*sinthet-FIz.*costhet;
Fdrd=Fdr.*sinthet-Fdz.*costhet;
durd=durdt.*sinthet-duzdt.*costhet;
FBthetd=FBthet;
FPthetd=FPthet;
FIthetd=FIthet;
Fdthetd=Fdthet;
duthetd=duthetdt;
FEzd=FEr.*costhet+FEz.*sinthet;
FBzd=FBr.*costhet+FBz.*sinthet;
FPzd=FPr.*costhet+FPz.*sinthet;
FIzd=FIr.*costhet+FIz.*sinthet;
Fdzd=Fdr.*costhet+Fdz.*sinthet;
duzd=durdt.*costhet+duzdt.*sinthet;
totforcerd=totforcer.*sinthet-totforcez.*costhet;
totforcethetd=totforcethet;
totforcezd=totforcer.*costhet+totforcez.*sinthet;
else
FErd=FEr;
FBrd=FBr;
FPrd=FPr;
FIrd=FIr;
Fdrd=Fdr;
durd=durdt;
FBthetd=FBthet;
FPthetd=FPthet;
FIthetd=FIthet;
Fdthetd=Fdthet;
duthetd=duthetdt;
FEzd=FEz;
FBzd=FBz;
FPzd=FPz;
FIzd=FIz;
Fdzd=Fdz;
duzd=duzdt;
totforcerd=totforcer;
totforcethetd=totforcethet;
totforcezd=totforcez;
end
papsize=[12 17];
%% Radial forces
f=figure('Name','radial');
ax0=subplot(3,2,1);
plotvalue(ax0,durd,'mn\partial_tu_{r}', 'F_{a,r} [N/m^3]',true,climsr,zlims,rlims);
ax1=subplot(3,2,2);
plotvalue(ax1,FErd,'qnE_r', 'F_{E,r} [N/m^3]',true,climsr,zlims,rlims);
ax2=subplot(3,2,3);
plotvalue(ax2,FBrd,'qn(u_\thetaB_z)','F_{B,r} [N/m^3]',true,climsr,zlims,rlims);
ax3=subplot(3,2,4);
plotvalue(ax3,FIrd,'-mn(\bf{u\bullet\nabla})u_r','F_{i,r} [N/m^3]',true,climsr,zlims,rlims);
ax4=subplot(3,2,5);
plotvalue(ax4,FPrd,'-(\nablaP)_r','F_{p,r} [N/m^3]',true,climsr,zlims,rlims);
ax5=subplot(3,2,6);
plotvalue(ax5,Fdrd,'-mnn_n<\sigma_d v>u_r','F_{d,r} [N/m^3]',true,climsr,zlims,rlims);
linkaxes([ax0 ax1 ax2 ax3 ax4 ax5],'xy')
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
M.savegraph(f,sprintf('%s/%sforce_balancer_%i_art',folder,name,floor(mean(it))),papsize)
%% Azimuthal forces
f=figure('Name','azimuthal');
ax0=subplot(3,2,1);
plotvalue(ax0,duzd,'mn\partial_tu_{\theta}','F_{a,\theta} [N/m^3]',true,climsthet,zlims,rlims);
ax1=subplot(3,2, 2);
plotvalue(ax1,ones(size(FEzd)),'qnE_\theta', 'F_{E,\theta} [N/m^3]',true,climsthet,zlims,rlims);
ax2=subplot(3,2,3);
plotvalue(ax2,FBthetd,'qn(u_zB_r-u_rB_z)','F_{B,\theta} [N/m^3]',true,climsthet,zlims,rlims);
ax3=subplot(3,2,4);
plotvalue(ax3,FIthetd,'-mn(\bf{u\bullet\nabla})u_\theta','F_{i,\theta} [N/m^3]',true,climsthet,zlims,rlims);
ax4=subplot(3,2,5);
plotvalue(ax4,FPthetd,'-(\nablaP)_\theta','F_{p,\theta} [N/m^3]',true,climsthet,zlims,rlims);
ax5=subplot(3,2,6);
plotvalue(ax5,Fdthetd,'-mnn_n<\sigma_d v>u_\theta','F_{d,\theta} [N/m^3]',true,climsthet,zlims,rlims);
linkaxes([ax0 ax1 ax2 ax3 ax4 ax5],'xy')
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
f.PaperSize=papsize;
M.savegraph(f,sprintf('%s/%sforce_balancethet_%i_art',folder,name,floor(mean(it))),papsize)
%% Axial forces
f=figure('Name','axial');
%zlim=[Fmin Fmax];
%zlim=[-15 15];
ax0=subplot(3,2,1);
plotvalue(ax0,duzd,'mn\partial_tu_{z}','F_{a,z} [N/m^3]',true,climsz,zlims,rlims);
ax1=subplot(3,2,2);
plotvalue(ax1,FEzd,'qnE_z', 'F_{E,z} [N/m^3]',true,climsz,zlims,rlims);
ax2=subplot(3,2,3);
plotvalue(ax2,FBzd,'-qn(u_\thetaB_r)','F_{B,z} [N/m^3]',true,climsz,zlims,rlims);
ax3=subplot(3,2,4);
plotvalue(ax3,FIzd,'-mn(\bf{u\bullet\nabla})u_z','F_{i,z} [N/m^3]',true,climsz,zlims,rlims);
ax4=subplot(3,2,5);
plotvalue(ax4,FPzd,'-(\nablaP)_z','F_{p,z} [N/m^3]',true,climsz,zlims,rlims);
ax5=subplot(3,2,6);
plotvalue(ax5,Fdzd,'-mnn_n<\sigma_d v>u_z','F_{d,z} [N/m^3]',true,climsz,zlims,rlims);
linkaxes([ax0 ax1 ax2 ax3 ax4 ax5],'xy')
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
f.PaperSize=papsize;
M.savegraph(f,sprintf('%s/%sforce_balancez_%i_art',folder,name,floor(mean(it))),papsize)
function plotvalue(ax,matrix,titl, clab, rm1strow, clims, zlims, rlims)
if nargin < 5
rm1strow=false;
end
if nargin>5
cmin=clims(1);
cmax=clims(2);
else
cmin=-inf;
cmax=inf;
end
matrix(M.geomweight(:,:,1)<=0)=NaN;
if rm1strow
if log
h=contourf(ax,M.zgrid*1e3,M.rgrid(1:end)*1e3,abs(matrix(1:end,:)),30,'linecolor','none');
else
h=contourf(ax,M.zgrid*1e3,M.rgrid(1:end)*1e3,matrix(1:end,:),30,'linecolor','none');
end
else
if log
h=contourf(ax,M.zgrid*1e3,M.rgrid*1e3,abs(matrix),30,'linecolor','none');
else
h=contourf(ax,M.zgrid*1e3,M.rgrid*1e3,matrix,30,'linecolor','none');
end
end
%set(h,'edgecolor','none');
hold on;
if log
zpos=interp2(M.zgrid*1e3, M.rgrid*1e3, abs(matrix), densitycontour(1,2:end)*1e3, densitycontour(2,2:end)*1e3);
else
zpos=interp2(M.zgrid*1e3, M.rgrid*1e3, matrix, densitycontour(1,2:end)*1e3, densitycontour(2,2:end)*1e3);
end
%plot3(ax,densitycontour(1,2:end)*1e3,densitycontour(2,2:end)*1e3,zpos,'-.','linewidth',2,'color',contourcolor)
contour(M.zgrid*1e3,M.rgrid*1e3,mean(N,3),[contourscale contourscale]*maxdens,'-.','linewidth',2,'color',contourcolor);
contour(M.zgrid*1e3,M.rgrid*1e3,M.geomweight(:,:,1),[0 0],'--','linewidth',2,'color','red');
%plot(ax,geometriccontour(1,2:end)*1e3,geometriccontour(2,2:end)*1e3,'--','linewidth',2,'color','red')
xlabel(ax,'z [mm]')
ylabel(ax,'r [mm]')
if (nargin>6 && ~isempty(zlims))
xlim(ax,zlims*1e3)
else
xlim(ax,[M.zgrid(1) M.zgrid(end)]*1e3)
end
if (nargin>7 && ~isempty(rlims))
ylim(ax,rlims*1e3)
else
if M.conformgeom
ylim(ax,[M.rgrid(1) M.rgrid(rgridmax)]*1e3)
else
ylim(ax,[M.rgrid(1) M.rgrid(end)]*1e3)
end
end
colormap(ax,'jet')
c = colorbar(ax);
c.Label.String='[N/m^{3}]';
caxis(ax,[cmin cmax])
if norm
caxis(ax,[Fminr Fmaxr]);
end
if log
set(ax,'colorscale','log')
caxis(ax,'auto')
end
title(ax,clab(1:end-8))
grid(ax, 'on')
view(ax,2)
% pos=ax.Position;
% pos(3)=0.98*pos(3);
% pos(4)=0.95*pos(4);
% ax.Position=pos;
end
end
c4science · Help