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

fennecs/matlab/ForceBalancespline.m
dd5453d97d74master

Last Change
Raw File

ForceBalancespline.m
View Options

function [Forces, Density]=ForceBalancespline(M,it,fdens)
%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 strcmp(it,':')
it=floor(0.95*size(M.t2d)):size(M.t2d);
end
if nargin<3
fdens=true;
end
fluiduThet=M.fluidUTHET(:,:,it);
Density.fluiduThet=fluiduThet;
me=M.me;
qe=M.qe;
N=M.N(:,:,it);
[R,~]=meshgrid(M.rgrid,M.zgrid);
Rinv=1./R';
Rinv(isinf(Rinv))=0;
Density.N=N;
Eforcer=-M.qe*M.Er(:,:,it);
Eforcez=-M.qe*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));
Bforcethet=zeros(size(N,1),size(N,2),size(N,3));
inertforcer=zeros(size(N,1),size(N,2),size(N,3));
inertforcez=zeros(size(N,1),size(N,2),size(N,3));
inertforcethet=zeros(size(N,1),size(N,2),size(N,3));
Pressforcer=zeros(size(N,1),size(N,2),size(N,3));
Pressforcethet=zeros(size(N,1),size(N,2),size(N,3));
Pressforcez=zeros(size(N,1),size(N,2),size(N,3));
durdt=zeros(size(N,1),size(N,2),size(N,3)-2);
duthetdt=zeros(size(N,1),size(N,2),size(N,3)-2);
duzdt=zeros(size(N,1),size(N,2),size(N,3)-2);
time=M.t2d(it);
parfor j=1:size(N,3)
Bforcer(:,:,j)=-qe.*fluiduThet(:,:,j).*M.Bz';
Bforcethet(:,:,j)=-qe.*(M.fluidUZ(:,:,it(j)).*M.Br'-M.fluidUR(:,:,it(j)).*M.Bz');
Bforcez(:,:,j)=qe.*fluiduThet(:,:,j).*M.Br';
inertforcer(:,:,j)=-me.*(-fluiduThet(:,:,j).^2.*Rinv+M.fluidUR(:,:,it(j)).*M.fluidUR.der(:,:,it(j),[1 0])+M.fluidUZ(:,:,it(j)).*M.fluidUR.der(:,:,it(j),[0 1]));
inertforcethet(:,:,j)=-me.*(M.fluidUR(:,:,it(j)).*fluiduThet(:,:,j).*Rinv+M.fluidUR(:,:,it(j)).*M.fluidUTHET.der(:,:,it(j),[1 0])+M.fluidUZ(:,:,it(j)).*M.fluidUTHET.der(:,:,it(j),[0 1]));
inertforcez(:,:,j)=-me.*(M.fluidUR(:,:,it(j)).*M.fluidUZ.der(:,:,it(j),[1 0])+M.fluidUZ(:,:,it(j)).*M.fluidUZ.der(:,:,it(j),[0 1]));
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])) )...
./N(:,:,j);
Pressforcethet(:,:,j)=-( squeeze(M.Presstens.der(2,:,:,it(j),[1 0]))...
+ 2*squeeze(M.Presstens(2,:,:,it(j))).*Rinv...
+ squeeze(M.Presstens.der(5,:,:,it(j),[0 1])) )...
./N(:,:,j);
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])) )...
./N(:,:,j);
end
parfor j=1:size(durdt,3)
durdt(:,:,j)=me*(M.fluidUR(:,:,it(j+2))-M.fluidUR(:,:,it(j)))/(time(j+2)-time(j));
duthetdt(:,:,j)=me*(M.fluidUTHET(:,:,it(j+2))-M.fluidUTHET(:,:,it(j)))/(time(j+2)-time(j));
duzdt(:,:,j)=me*(M.fluidUZ(:,:,it(j+2))-M.fluidUZ(:,:,it(j)))/(time(j+2)-time(j));
end
Forces.it=it;
if(~fdens)
Forces.Eforcer=Eforcer;
Forces.Eforcez=Eforcez;
Forces.Bforcer=Bforcer;
Forces.Bforcethet=Bforcethet;
Forces.Bforcez=Bforcez;
Forces.inertforcer=inertforcer;
Forces.inertforcethet=inertforcethet;
Forces.inertforcez=inertforcez;
Forces.Pressforcer=Pressforcer;
Forces.Pressforcethet=Pressforcethet;
Forces.Pressforcez=Pressforcez;
Forces.durdt=durdt;
Forces.duthetdt=duthetdt;
Forces.duzdt=duzdt;
else
Forces.Eforcer=Eforcer.*N;
Forces.Eforcez=Eforcez.*N;
Forces.Bforcer=Bforcer.*N;
Forces.Bforcethet=Bforcethet.*N;
Forces.Bforcez=Bforcez.*N;
Forces.inertforcer=inertforcer.*N;
Forces.inertforcethet=inertforcethet.*N;
Forces.inertforcez=inertforcez.*N;
Forces.Pressforcer=Pressforcer.*N;
Forces.Pressforcethet=Pressforcethet.*N;
Forces.Pressforcez=Pressforcez.*N;
Forces.durdt=durdt.*N(:,:,2:end-1);
Forces.duthetdt=duthetdt.*N(:,:,2:end-1);
Forces.duzdt=duzdt.*N(:,:,2:end-1);
end
end
c4science · Help