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

fennecs/matlab/FBspline.m
ddfe457d9de0geom_testmumps

Last Change
Raw File

FBspline.m
View Options

function [Forces, Continuity]=FBspline(M,it,Forces, Continuity, Density,norm,log,contourscale, zlims)
%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
if it==':'
it=floor(0.95*size(M.t2d)):size(M.t2d);
end
switch nargin
case 5
contourscale=0.6;
norm=false;
log=false;
zlims=[-inf inf];
case 6
log=false;
contourscale=0.6;
zlims=[-inf inf];
case 7
contourscale=0.6;
zlims=[-inf inf];
case 8
zlims=[-inf inf];
case 9
otherwise
error("Invalid number of arguments")
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.inertforce(:)),max(Forces.Pressforcer(:))]);
Fmaxz=max([max(Forces.Eforcez(:)),max(Forces.Bforcez(:)),max(Forces.Pressforcez(:))]);
if log
Er=Forces.Eforcer(Forces.Eforcer~=0);
Br=Forces.Bforcer(Forces.Bforcer~=0);
inert=Forces.inertforce(Forces.inertforce~=0);
pressr=Forces.Pressforcer(Forces.Pressforcer~=0);
Fminr=min([min(abs(Er(:))),min(abs(Br(:))),min(abs(inert(:))),min(abs(pressr(:)))]);
else
Fminr=min([min(Forces.Eforcer(:)),min(Forces.Bforcer(:)),min(Forces.inertforce(:)),min(Forces.Pressforcer(:))]);
end
rgridmax=sum(M.nnr(1:2));
%% Radial forces
f=figure();
ax1=subplot(4,1,1);
plotvalue(ax1,mean(Forces.Eforcer,3),'Electric force density', 'F_{Er} [Nm^{-3}]',true);
%zlim=[Fmin Fmax];
zlim=[-inf inf];
ax2=subplot(4,1,2);
plotvalue(ax2,mean(Forces.Bforcer,3),'Magnetic force density','F_{Br} [Nm^{-3}]',true,zlim)
ax3=subplot(4,1,3);
plotvalue(ax3,mean(Forces.inertforce,3),'Centrifugal force density','F_{i} [Nm^{-3}]',true,zlim)
ax4=subplot(4,1,4);
plotvalue(ax4,mean(Forces.Pressforcer,3),'Pressure force density','F_{pr} [Nm^{-3}]',true,zlim)
linkaxes([ax1 ax2 ax3 ax4],'xy')
sgtitle(sprintf('Radial 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_balancer',name),'-dpdf','-fillpage')
%% Radial relative forces
f=figure();
ax1=subplot(4,1,1);
color_limits=[-1 1];
reference=Forces.Bforcer;
plotvalue(ax1,mean(Forces.Eforcer./reference,3),'Relative Electric force density', 'F_{Er}/F_{Br} [Nm^{-3}]',true,color_limits);
%zlim=[Fmin Fmax];
zlim=[-inf inf];
ax2=subplot(4,1,2);
plotvalue(ax2,mean(Forces.Bforcer./reference,3),'Relative Magnetic force density','F_{Br}/F_{Br} [Nm^{-3}]',true,color_limits)
ax3=subplot(4,1,3);
plotvalue(ax3,mean(Forces.inertforce./reference,3),'Relative Centrifugal force density','F_i/F_{Br} [Nm^{-3}]',true,color_limits)
ax4=subplot(4,1,4);
plotvalue(ax4,mean(Forces.Pressforcer./reference,3),'Relative Pressure force density','F_{pr}/F_{Br} [Nm^{-3}]',true,color_limits)
linkaxes([ax1 ax2 ax3 ax4],'xy')
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)))
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[20 29];
f.PaperSize=papsize;
print(f,sprintf('%s_relatforce_balance',name),'-dpdf','-fillpage')
%% Axial forces
f=figure();
ax1=subplot(3,1,1);
plotvalue(ax1,mean(Forces.Eforcez,3),'Electric force density', 'F_{Ez} [Nm^{-3}]',true);
zlim=[-inf inf];
ax2=subplot(3,1,2);
plotvalue(ax2,mean(Forces.Bforcez,3),'Magnetic force density','F_{Bz} [Nm^{-3}]',true,zlim)
ax3=subplot(3,1,3);
plotvalue(ax3,mean(Forces.Pressforcez,3),'Pressure force density','F_{pz} [Nm^{-3}]',true,zlim)
linkaxes([ax1 ax2 ax3],'xy')
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)))
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[20 29];
f.PaperSize=papsize;
print(f,sprintf('%sforce_balancez',name),'-dpdf','-fillpage')
%% Forces sums
f2=figure();
ax4=subplot(3,1,1);
totforcer=mean(Forces.Eforcer+Forces.Bforcer+Forces.inertforce+Forces.Pressforcer,3);
Forcemaxr=cat(3,mean(Forces.Eforcer,3),mean(Forces.Bforcer,3),mean(Forces.inertforce,3),mean(Forces.Pressforcer,3));
Forcemaxr=max(Forcemaxr,[],3);
plotvalue(ax4,totforcer./Forcemaxr,'Total radial force density','F [Nm^{-3}]',true,zlims)
ax5=subplot(3,1,2);
totforcez=mean(Forces.Eforcez+Forces.Bforcez+Forces.Pressforcez,3);
Forcemaxz=cat(3,mean(Forces.Eforcez,3),mean(Forces.Bforcez,3),mean(Forces.Pressforcez,3));
Forcemaxz=max(Forcemaxz,[],3);
plotvalue(ax5,totforcez./Forcemaxz,'Total axial force density','F [Nm^{-3}]',true,zlims)
ax6=subplot(3,1,3);
surface(ax6,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(ax6,densitycontour(1,2:end),densitycontour(2,2:end),zlevel,'-.','linewidth',2,'color',contourcolor)
xlim(ax6,[M.zgrid(1) M.zgrid(end)])
ylim(ax6,[M.rgrid(1) M.rgrid(rgridmax)])
colormap(ax6,'parula')
xlabel(ax6,'z [m]')
ylabel(ax6,'r [m]')
title(ax6,'Density')
c = colorbar(ax6);
c.Label.String= 'n[m^{-3}]';
view(ax6,2)
sgtitle(sprintf('Radial 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')
%% Continuity
f3=figure();
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(Continuity.ndivu,3),'n\nabla(u)','n\nabla(u)[m^{-3}s^{-1}]',true,zlim)
ax9=subplot(4,1,3);
plotvalue(ax9,mean(Continuity.ugradn,3),'u\nabla(n)','u\nabla(n)[m^{-3}s^{-1}]',true,zlim)
ax10=subplot(4,1,4);
plotvalue(ax10,mean(Continuity.ndivu+Continuity.ugradn+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')
function plotvalue(axis,matrix,titl, clab, rm1strow, zlims)
if nargin < 5
rm1strow=false;
end
if nargin>5
zmin=zlims(1);
zmax=zlims(2);
else
zmin=-inf;
zmax=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]')
xlim(axis,[M.zgrid(1) M.zgrid(end)])
ylim(axis,[M.rgrid(1) M.rgrid(rgridmax)])
colormap(axis,'jet')
c = colorbar(axis);
c.Label.String=clab;
caxis(axis,[zmin zmax])
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