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

fennecs/matlab/ForceBalancespline.m
a2847b9464fageom_testpardiso

Last Change
Raw File

ForceBalancespline.m
View Options

function [Forces, Density]=ForceBalancespline(M,it,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 2
contourscale=0.6;
norm=false;
log=false;
zlims=[-inf inf];
case 3
log=false;
contourscale=0.6;
zlims=[-inf inf];
case 4
contourscale=0.6;
zlims=[-inf inf];
case 5
zlims=[-inf inf];
case 6
otherwise
error("Invalid number of arguments")
end
contourcolor=[0 0 0];%[255 20 147]/255;
fluiduThet=M.fluidUTHET(:,:,it);
Density.fluiduThet=fluiduThet;
N=M.N(:,:,it);
[R,~]=meshgrid(M.rgrid,M.zgrid);
Rinv=1./R';
Rinv(isinf(Rinv))=0;
Density.N=N;
maxdens=max(N(:));
densitycontour=contour(M.zgrid,M.rgrid,mean(N,3),[contourscale contourscale]*maxdens);
Forces.Eforcer=-M.qe*N.*M.Er(:,:,it);
Forces.Eforcez=-M.qe*N.*M.Ez(:,:,it);
Bforcer=zeros(size(N,1),size(N,2),size(N,3));
Bforcez=zeros(size(N,1),size(N,2),size(N,3));
inertforce=zeros(size(N,1),size(N,2),size(N,3));
Pressforcer=zeros(size(N,1),size(N,2),size(N,3));
Pressforcez=zeros(size(N,1),size(N,2),size(N,3));
parfor j=1:size(N,3)
Bforcer(:,:,j)=-M.qe*N(:,:,j).*fluiduThet(:,:,j).*M.Bz';
Bforcez(:,:,j)=M.qe*N(:,:,j).*fluiduThet(:,:,j).*M.Br';
inertforce(:,:,j)=M.me*N(:,:,j).*fluiduThet(:,:,j).^2.*Rinv;
Pressforcer(:,:,j)=-( squeeze(M.Presstens.der(1,:,:,it(j),[1 0]))...
+ squeeze(M.Presstens(1,:,:,it(j)) - M.Presstens(4,:,:,it(j))).*Rinv...
+ squeeze(M.Presstens.der(3,:,:,it(j),[0 1])) );
Pressforcez(:,:,j)=-( squeeze(M.Presstens.der(3,:,:,it(j),[1 0]))...
+ squeeze(M.Presstens(3,:,:,it(j))).*Rinv...
+ squeeze(M.Presstens.der(6,:,:,it(j),[0 1])) );
end
Forces.Bforcer=Bforcer;
Forces.Bforcez=Bforcez;
Forces.inertforce=inertforce;
Forces.Pressforcer=Pressforcer;
Forces.Pressforcez=Pressforcez;
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(Eforcer~=0);
Br=Forces.Bforcer(Bforcer~=0);
inert=Forces.inertforce(inertforce~=0);
pressr=Forces.Pressforcer(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}]',false);
%
% %zlim=[Fmin Fmax];
% zlim=[-inf inf];
%
% ax2=subplot(4,1,2);
% plotvalue(ax2,mean(Forces.Bforcer,3),'Magnetic force density','F_{Br} [Nm^{-3}]',false,zlim)
%
% ax3=subplot(4,1,3);
% plotvalue(ax3,mean(Forces.inertforce,3),'Centrifugal force density','F_{i} [Nm^{-3}]',false,zlim)
%
% ax4=subplot(4,1,4);
% plotvalue(ax4,mean(Forces.Pressforcer,3),'Pressure force density','F_{pr} [Nm^{-3}]',false,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);
% reference=mean(Forces.Bforcer,3);
% plotvalue(ax1,mean(Forces.Eforcer,3)./reference,'Relative Electric force density', 'F_{Er}/F_{Br}',false);
%
% %zlim=[Fmin Fmax];
% zlim=[-inf inf];
%
% ax2=subplot(4,1,2);
% plotvalue(ax2,mean(Forces.Bforcer,3)./reference,'Relative Magnetic force density','F_{Br}/F_{Br}',false)
%
% ax3=subplot(4,1,3);
% plotvalue(ax3,mean(Forces.inertforce,3)./reference,'Relative Centrifugal force density','F_i/F_{Br}',false)
%
% ax4=subplot(4,1,4);
% plotvalue(ax4,mean(Forces.Pressforcer,3)./reference,'Relative Pressure force density','F_{pr}/F_{Br}',false)
%
% 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}]',false);
%
% zlim=[-inf inf];
%
% ax2=subplot(3,1,2);
% plotvalue(ax2,mean(Forces.Bforcez,3),'Magnetic force density','F_{Bz} [Nm^{-3}]',false,zlim)
%
% ax3=subplot(3,1,3);
% plotvalue(ax3,mean(Forces.Pressforcez,3),'Pressure force density','F_{pz} [Nm^{-3}]',false,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}]',false,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}]',false,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')
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(2:end),abs(matrix(2:end,:)));
else
h=surface(axis,M.zgrid,M.rgrid(2:end),matrix(2: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')
end
title(axis,titl)
grid(axis, 'on')
view(axis,2)
end
end
c4science · Help