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

fennecs/matlab/FBspline.m
118e68d1a0d5breakdown

Last Change
Raw File

FBspline.m
View Options

function [Forces, Continuity]=FBspline(M,Forces, Continuity, Density,norm,log,parper,contourscale, zlims, rlims)
%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.6;
norm=false;
log=false;
zlims=[-inf inf];
case 5
parper=false;
log=false;
contourscale=0.6;
zlims=[-inf inf];
case 6
parper=false;
contourscale=0.6;
zlims=[-inf inf];
case 7
contourscale=0.6;
zlims=[-inf inf];
case 8
zlims=[-inf inf];
case 9
rlims=[];
case 10
end
if nargin<10
rlims=[];
end
contourcolor=[0 0 0];%[255 20 147]/255;
N=Density.N;
maxdens=max(N(:));
densitycontour=contour(M.zgrid,M.rgrid,mean(N,3),[contourscale contourscale]*maxdens);
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=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
%% Radial forces
f=figure('Name','radial');
clims=[-inf inf];
ax0=subplot(6,1,1);
plotvalue(ax0,durd,'Mass flux time derivative', 'mn\partialu_{r}/\partialt [Nm^{-3}]',true,clims,zlims,rlims);
ax1=subplot(6,1,2);
plotvalue(ax1,FErd,'Electric force density', 'F_{E,r} [Nm^{-3}]',true,clims,zlims,rlims);
ax2=subplot(6,1,3);
plotvalue(ax2,FBrd,'Magnetic force density','F_{B,r} [Nm^{-3}]',true,clims,zlims,rlims);
ax3=subplot(6,1,4);
plotvalue(ax3,FIrd,'Inertial force density','F_{i,r} [Nm^{-3}]',true,clims,zlims,rlims);
ax4=subplot(6,1,5);
plotvalue(ax4,FPrd,'Pressure force density','F_{p,r} [Nm^{-3}]',true,clims,zlims,rlims);
ax5=subplot(6,1,6);
plotvalue(ax5,Fdrd,'Drag force density','F_{d,r} [Nm^{-3}]',true,clims,zlims,rlims);
linkaxes([ax0 ax1 ax2 ax3 ax4 ax5],'xy')
if parper
sgtitle(f,sprintf('Perp forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
else
sgtitle(f,sprintf('Radial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
end
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[20 29];
f.PaperSize=papsize;
print(f,sprintf('%sforce_balancer',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sforce_balancer',name))
%% Radial relative forces
f=figure('Name','radial relative');
ax1=subplot(5,1,1);
color_limits=[-1 1]*10;
referencer=FBrd;
plotvalue(ax1,FErd./referencer,'Relative Electric force density', 'F_{Er}/F_{Br}',true,color_limits,zlims,rlims);
%zlim=[Fmin Fmax];
zlim=[-inf inf];
ax2=subplot(5,1,2);
plotvalue(ax2,FBrd./referencer,'Relative Magnetic force density','F_{Br}/F_{Br}',true,color_limits,zlims,rlims);
ax3=subplot(5,1,3);
plotvalue(ax3,FIrd./referencer,'Relative Inertial force density','F_i/F_{Br}',true,color_limits,zlims,rlims);
ax4=subplot(5,1,4);
plotvalue(ax4,FPrd./referencer,'Relative Pressure force density','F_{pr}/F_{Br}',true,color_limits,zlims,rlims);
ax5=subplot(5,1,5);
plotvalue(ax5,Fdrd./referencer,'Relative Drag force density','F_{dr}/F_{Br}',true,color_limits,zlims,rlims);
linkaxes([ax1 ax2 ax3 ax4 ax5],'xy')
if parper
sgtitle(sprintf('Relative perp forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
else
sgtitle(sprintf('Relative radial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
end
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[20 29];
f.PaperSize=papsize;
print(f,sprintf('%s_relatforce_balancer',name),'-dpdf','-fillpage')
savefig(f,sprintf('%s_relatforce_balancer',name))
%% Azimuthal forces
f=figure('Name','azimuthal');
ax1=subplot(5,1,1);
plotvalue(ax1,duthetdt,'Mass flux time derivartive','mn\partialu_\theta/\partialt [Nm^{-3}]',true,clims,zlims,rlims);
ax2=subplot(5,1,2);
plotvalue(ax2,FBthetd,'Magnetic force density','F_{B,\theta} [Nm^{-3}]',true,clims,zlims,rlims);
ax3=subplot(5,1,3);
plotvalue(ax3,FIthetd,'Inertial force density','F_{i,\theta} [Nm^{-3}]',true,clims,zlims,rlims);
ax4=subplot(5,1,4);
plotvalue(ax4,FPthetd,'Pressure force density','F_{p,\theta} [Nm^{-3}]',true,clims,zlims,rlims);
ax5=subplot(5,1,5);
plotvalue(ax5,Fdthetd,'Drag force density','F_{d,\theta} [Nm^{-3}]',true,clims,zlims,rlims);
linkaxes([ax1 ax2 ax3 ax4 ax5],'xy')
sgtitle(sprintf('Azimuthal forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[20 29];
f.PaperSize=papsize;
print(f,sprintf('%sforce_balancethet',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sforce_balancethet',name))
%% Axial forces
f=figure('Name','axial');
%zlim=[Fmin Fmax];
clims=[-inf inf];
%zlim=[-15 15];
ax0=subplot(6,1,1);
plotvalue(ax0,duzd,'Mass flux time derivartive','mn\partialu_z/\partialt [Nm^{-3}]',true,clims,zlims,rlims);
ax1=subplot(6,1,2);
plotvalue(ax1,FEzd,'Electric force density', 'F_{E,z} [Nm^{-3}]',true,clims,zlims,rlims);
ax2=subplot(6,1,3);
plotvalue(ax2,FBzd,'Magnetic force density','F_{B,z} [Nm^{-3}]',true,clims,zlims,rlims);
ax3=subplot(6,1,4);
plotvalue(ax3,FIzd,'Inertial force density','F_{i,z} [Nm^{-3}]',true,clims,zlims,rlims);
ax4=subplot(6,1,5);
plotvalue(ax4,FPzd,'Pressure force density','F_{p,z} [Nm^{-3}]',true,clims,zlims,rlims);
ax5=subplot(6,1,6);
plotvalue(ax5,Fdzd,'Drag force density','F_{d,z} [Nm^{-3}]',true,clims,zlims,rlims);
linkaxes([ax0 ax1 ax2 ax3 ax4 ax5],'xy')
if parper
sgtitle(sprintf('Parallel forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
else
sgtitle(sprintf('Axial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
end
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[20 29];
f.PaperSize=papsize;
print(f,sprintf('%sforce_balancez',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sforce_balancez',name))
%% Axial relative forces
f=figure('Name','axial relative');
ax1=subplot(5,1,1);
color_limits=[-1 1]*10;
referencez=FPzd;
plotvalue(ax1,FEzd./referencez,'Relative Electric force density', 'F_{Ez}/F_{Ez}',true,color_limits,zlims,rlims);
%zlim=[Fmin Fmax];
zlim=[-inf inf];
ax2=subplot(5,1,2);
plotvalue(ax2,FBzd./referencez,'Relative Magnetic force density','F_{Bz}/F_{Ez}',true,color_limits,zlims,rlims);
ax3=subplot(5,1,3);
plotvalue(ax3,FIzd./referencez,'Relative Inertial force density','F_{Iz}/F_{Ez}',true,color_limits,zlims,rlims);
ax4=subplot(5,1,4);
plotvalue(ax4,FPzd./referencez,'Relative Pressure force density','F_{Pz}/F_{Ez}',true,color_limits,zlims,rlims);
ax5=subplot(5,1,5);
plotvalue(ax5,Fdzd./referencez,'Relative Drag force density','F_{dz}/F_{Ez}',true,color_limits,zlims,rlims);
linkaxes([ax1 ax2 ax3 ax4 ax5],'xy')
if parper
sgtitle(sprintf('Relative Parallel forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
else
sgtitle(sprintf('Relative axial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
end
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[20 29];
f.PaperSize=papsize;
print(f,sprintf('%s_relatforce_balancez',name),'-dpdf','-fillpage')
savefig(f,sprintf('%s_relatforce_balancez',name))
%% Forces sums
f2=figure('Name','total');
ax4=subplot(4,1,1);
maxforcerd=1;%cat(3,FErd,FBrd,FPrd,FIrd);
%maxforcerd=max(maxforcerd,[],3);
if parper
plotvalue(ax4,totforcerd./maxforcerd,'Total Perp force density','F [Nm^{-3}]',true,clims,zlims,rlims);
else
plotvalue(ax4,totforcerd./maxforcerd,'Total radial force density','F [Nm^{-3}]',true,clims,zlims,rlims);
end
ax5=subplot(4,1,2);
maxforcethetd=1;%cat(3,FEzd,FBzd,FPzd,FIzd);
%maxforcezd=max(maxforcezd,[],3);
%zlims=[-15 15];
plotvalue(ax5,totforcethetd./maxforcethetd,'Total azimuthal force density','F [Nm^{-3}]',true,clims,zlims,rlims);
ax6=subplot(4,1,3);
maxforcezd=1;%cat(3,FEzd,FBzd,FPzd,FIzd);
%maxforcezd=max(maxforcezd,[],3);
%zlims=[-15 15];
if parper
plotvalue(ax6,totforcezd./maxforcezd,'Total par force density','F [Nm^{-3}]',true,clims,zlims,rlims);
else
plotvalue(ax6,totforcezd./maxforcezd,'Total axial force density','F [Nm^{-3}]',true,clims,zlims,rlims);
end
ax7=subplot(4,1,4);
surface(ax7,M.zgrid,M.rgrid,mean(N,3),'edgecolor','none');
hold on;
zlevel=interp2(M.zgrid, M.rgrid, mean(N,3), densitycontour(1,2:end), densitycontour(2,2:end));
plot3(ax7,densitycontour(1,2:end),densitycontour(2,2:end),zlevel,'-.','linewidth',2,'color',contourcolor)
colormap(ax7,'parula')
xlabel(ax7,'z [m]')
ylabel(ax7,'r [m]')
title(ax7,'Density')
c = colorbar(ax7);
c.Label.String= 'n[m^{-3}]';
view(ax7,2)
linkaxes([ax4 ax5 ax6 ax7])
if (nargin>8 && ~isempty(zlims))
xlim(ax7,zlims)
else
xlim(ax7,[M.zgrid(1) M.zgrid(end)])
end
if (nargin>9 && ~isempty(rlims))
ylim(ax7,rlims)
else
if M.conformgeom
ylim(ax7,[M.rgrid(1) M.rgrid(rgridmax)])
else
ylim(ax7,[M.rgrid(1) M.rgrid(end)])
end
end
sgtitle(sprintf('Forces sum t=[%1.2g-%1.2g]s n_e=%1.2g m^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
f2.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f2.PaperUnits='centimeters';
papsize=[20 29];
f2.PaperSize=papsize;
print(f2,sprintf('%sforce_dens',name),'-dpdf','-fillpage')
savefig(f2,sprintf('%sforce_dens',name))
%% Continuity
if( ~isempty(Continuity))
f3=figure('Name','Continuity');
it=Continuity.it;
t=M.t2d(it);
ax7=subplot(4,1,1);
plotvalue(ax7,mean(Continuity.dndt,3),'\partialn/\partialt','\partialn/\partialt [m^{-3}s^{-1}]',true)
zlim=[max([min(min(mean(Continuity.ndivu(2:end,:,:),3))),min(min(mean(Continuity.ugradn(2:end,:,:),3)))]) min([max(max(mean(Continuity.ndivu(2:end,:,:),3))),max(max(mean(Continuity.ugradn(2:end,:,:),3)))])];
ax8=subplot(4,1,2);
plotvalue(ax8,mean(sum(Continuity.ndivu,4),3),'n\nabla(u)','n\nabla(u)[m^{-3}s^{-1}]',true,zlim)
ax9=subplot(4,1,3);
plotvalue(ax9,mean(sum(Continuity.ugradn,4),3),'u\nabla(n)','u\nabla(n)[m^{-3}s^{-1}]',true,zlim)
ax10=subplot(4,1,4);
plotvalue(ax10,mean(sum(Continuity.ndivu,4)+sum(Continuity.ugradn,4)+Continuity.dndt,3),'Total','total[m^{-3}s^{-1}]',true,zlim)
sgtitle(sprintf('Continuity t=[%1.2g-%1.2g]s n_e=%1.2g m^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
f3.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f3.PaperUnits='centimeters';
papsize=[20 19];
f3.PaperSize=papsize;
print(f3,sprintf('%scontinuity',name),'-dpdf','-fillpage')
savefig(f,sprintf('%scontinuity',name))
end
function plotvalue(axis,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
if rm1strow
if log
h=surface(axis,M.zgrid,M.rgrid(12:end),abs(matrix(12:end,:)));
else
h=surface(axis,M.zgrid,M.rgrid(12:end),matrix(12:end,:));
end
else
if log
h=surface(axis,M.zgrid,M.rgrid,abs(matrix));
else
h=surface(axis,M.zgrid,M.rgrid,matrix);
end
end
set(h,'edgecolor','none');
hold on;
if log
zpos=interp2(M.zgrid, M.rgrid, abs(matrix), densitycontour(1,2:end), densitycontour(2,2:end));
else
zpos=interp2(M.zgrid, M.rgrid, matrix, densitycontour(1,2:end), densitycontour(2,2:end));
end
plot3(axis,densitycontour(1,2:end),densitycontour(2,2:end),zpos,'-.','linewidth',2,'color',contourcolor)
xlabel(axis,'z [m]')
ylabel(axis,'r [m]')
if (nargin>6 && ~isempty(zlims))
xlim(axis,zlims)
else
xlim(axis,[M.zgrid(1) M.zgrid(end)])
end
if (nargin>7 && ~isempty(rlims))
ylim(axis,rlims)
else
if M.conformgeom
ylim(axis,[M.rgrid(1) M.rgrid(rgridmax)])
else
ylim(axis,[M.rgrid(1) M.rgrid(end)])
end
end
colormap(axis,'jet')
c = colorbar(axis);
c.Label.String=clab;
caxis(axis,[cmin cmax])
if norm
caxis(axis,[Fminr Fmaxr]);
end
if log
set(axis,'colorscale','log')
caxis(axis,'auto')
end
title(axis,titl)
grid(axis, 'on')
view(axis,2)
end
end
c4science · Help