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

fennecs/matlab/load_espic2d.m
50e38dd6641fmaster

Last Change
Raw File

load_espic2d.m
View Options

function Data = load_espic2d( file, Data, type, partstime)
% Old function used to load simulations data from espic2d HDF5 results
% This function does not handle b-splines data apart from the density and load the full file to memory
% This can be problematic for big files
if(nargin==1)
type='all';
else
filedata=dir(file);
if (isempty(filedata))
error("File: ""%s"" doesn't exist",file)
end
timestamp=filedata.date;
try
if (strcmp(Data.file,file) && ~max(timestamp > Data.timestamp) && (strcmp(Data.type,type)||strcmp(Data.type,'all')) )
return
end
catch
end
if(nargin<3)
type='all';
end
if(nargin<4)
partstime=[0 1];
end
end
Data.type=type;
Data.file=file;
% Physical constants
Data.vlight=299792458;
Data.qe=1.60217662E-19;
Data.me=9.109383E-31;
Data.eps_0=8.85418781762E-12;
Data.kb=1.38064852E-23;
%Data.job_time = h5readatt(file,'/data/input.00/','job_time');
%Data.extra_time = h5readatt(file,'/data/input.00/','extra_time');
Data.dt = h5readatt(file,'/data/input.00/','dt');
Data.tmax= h5readatt(file,'/data/input.00/','tmax');
Data.nrun = h5readatt(file,'/data/input.00/','nrun');
Data.nlres = strcmp(h5readatt(file,'/data/input.00/','nlres'),'y');
Data.nlsave = strcmp(h5readatt(file,'/data/input.00/','nlsave'),'y');
Data.nlclassical =strcmp(h5readatt(file,'/data/input.00/','nlclassical'),'y');
Data.nlPhis =strcmp(h5readatt(file,'/data/input.00/','nlPhis'),'y');
Data.nz = h5readatt(file,'/data/input.00/','nz');
Data.nnr = h5read(file,'/data/input.00/nnr');
Data.femorder = h5read(file,'/data/input.00/femorder');
Data.lz = h5read(file,'/data/input.00/lz');
Data.lr = h5read(file,'/data/input.00/radii');
Data.qsim = h5readatt(file,'/data/input.00/','qsim');
Data.msim = h5readatt(file,'/data/input.00/','msim');
Data.nplasma = h5readatt(file,'/data/input.00/','nplasma');
Data.potinn = h5readatt(file,'/data/input.00/','potinn');
Data.potout = h5readatt(file,'/data/input.00/','potout');
Data.B0 = h5readatt(file,'/data/input.00/','B0');
Data.Rcurv = h5readatt(file,'/data/input.00/','Rcurv');
Data.width = h5readatt(file,'/data/input.00/','width');
Data.n0 = h5readatt(file,'/data/input.00/','n0');
Data.temp = h5readatt(file,'/data/input.00/','temp');
try
Data.it0 = h5readatt(file,'/data/input.00/','it0d');
Data.it1 = h5readatt(file,'/data/input.00/','it2d');
Data.it2 = h5readatt(file,'/data/input.00/','itparts');
catch
Data.it0 = h5readatt(file,'/data/input.00/','it0');
Data.it1 = h5readatt(file,'/data/input.00/','it1');
Data.it1 = h5readatt(file,'/data/input.00/','it2');
end
Data.omepe=sqrt(abs(Data.n0)*Data.qe^2/(Data.me*Data.eps_0));
Data.omece=Data.qe*Data.B0/Data.me;
% Normalizations
Data.tnorm=1/Data.omepe;
Data.rnorm=Data.vlight*Data.tnorm;
Data.bnorm=Data.B0;
Data.enorm=Data.vlight*Data.bnorm;
Data.phinorm=Data.enorm*Data.rnorm;
Data.vnorm=Data.vlight;
% Grid data
Data.rgrid= h5read(file, '/data/var1d/rgrid')*Data.rnorm;
Data.zgrid= h5read(file, '/data/var1d/zgrid')*Data.rnorm;
Data.dz=(Data.zgrid(end)-Data.zgrid(1))/double(Data.nz);
Data.dr=(Data.lr(2:end)-Data.lr(1:end-1))./double(Data.nnr);
Br = h5read(file,'/data/fields/Br')*Data.bnorm;
Data.Br= reshape(Br,length(Data.zgrid),length(Data.rgrid));
Bz = h5read(file,'/data/fields/Bz')*Data.bnorm;
Data.Bz= reshape(Bz,length(Data.zgrid),length(Data.rgrid));
try
Atheta = h5read(file,'/data/fields/Athet')*Data.bnorm;
Data.Athet= reshape(Atheta,length(Data.zgrid),length(Data.rgrid));
catch
end
Data.B=sqrt(Data.Bz.^2+Data.Br.^2);
clear Br Bz
Data.femorder = h5read(file,'/data/input.00/femorder');
Data.ngauss = h5read(file,'/data/input.00/ngauss');
Data.plasmadim = h5read(file,'/data/input.00/plasmadim');
Data.radii = h5read(file,'/data/input.00/radii');
Data.epot = h5read(file,'/data/var0d/epot');
Data.ekin = h5read(file,'/data/var0d/ekin');
Data.etot = h5read(file,'/data/var0d/etot');
try
Data.etot0 = h5read(file,'/data/var0d/etot0');
Data.eerr = Data.etot-Data.etot0;
catch
Data.eerr = Data.etot-Data.etot(2);
end
if(strcmp(type,'all') || strcmp(type,'fields'))
pot = h5read(file,'/data/fields/pot')*Data.phinorm;
Data.pot= reshape(pot,length(Data.zgrid),length(Data.rgrid),size(pot,2));
clear pot
Data.pot=permute(Data.pot,[2,1,3]);
Er = h5read(file,'/data/fields/Er')*Data.enorm;
Data.Er= reshape(Er,length(Data.zgrid),length(Data.rgrid),size(Er,2));
clear Er
Data.Er=permute(Data.Er,[2,1,3]);
Ez = h5read(file,'/data/fields/Ez')*Data.enorm;
Data.Ez= reshape(Ez,length(Data.zgrid),length(Data.rgrid),size(Ez,2));
clear Ez
Data.Ez=permute(Data.Ez,[2,1,3]);
try
dens = h5read(file,'/data/fields/moments',[1 1 1],[1,Inf,Inf]);
kr=Data.femorder(2)+1;
Data.knotsr=augknt(Data.rgrid,kr);
kz=Data.femorder(1)+1;
Data.knotsz=augknt(Data.zgrid,kz);
zvol=fnder(spmak(Data.knotsz,ones(1,length(Data.zgrid))), -1 );
rvol=fnder(spmak(Data.knotsr,2*pi*Data.rgrid'), -1 );
ZVol=diff(fnval(zvol,Data.knotsz));
RVol=diff(fnval(rvol,Data.knotsr));
Data.V=ZVol(2:end-1)*RVol(2:end-1)';
Data.invV=1./Data.V';
Data.invV(isinf(Data.invV))=0;
dens = h5read(file,'/data/fields/moments',[1 1 1],[1,Inf,Inf]);
Data.dens=reshape(dens,Data.nz+Data.femorder(1),sum(Data.nnr)+Data.femorder(2),[]);
Data.dens=permute(Data.dens,[2,1,3]);
clear dens
fluidvr = h5read(file,'/data/fields/moments',[2 1 1],[1,Inf,Inf]);
Data.fluidvr=reshape(fluidvr,Data.nz+Data.femorder(1),sum(Data.nnr)+Data.femorder(2),[]);
clear fluidvr
Data.fluidvr=permute(Data.fluidvr,[2,1,3])*Data.vnorm./Data.dens;
Data.fluidvr(isinf(Data.fluidvr))=0;
fluidvthet = h5read(file,'/data/fields/moments',[3 1 1],[1,Inf,Inf]);
Data.fluidvthet=reshape(fluidvthet,Data.nz+Data.femorder(1),sum(Data.nnr)+Data.femorder(2),[]);
Data.fluidvthet=permute(Data.fluidvthet,[2,1,3])*Data.vnorm./Data.dens;
Data.fluidvthet(isinf(Data.fluidvthet))=0;
clear fluidvthet
fluidvz = h5read(file,'/data/fields/moments',[4 1 1],[1,Inf,Inf]);
Data.fluidvz=reshape(fluidvz,Data.nz+Data.femorder(1),sum(Data.nnr)+Data.femorder(2),[]);
Data.fluidvz=permute(Data.fluidvz,[2,1,3])*Data.vnorm./Data.dens;
Data.fluidvz(isinf(Data.fluivz))=0;
clear fluidvz
vrvr = h5read(file,'/data/fields/moments',[5 1 1],[1,Inf,Inf]);
Data.fluidvz=reshape(vrvr,Data.nz+Data.femorder(1),sum(Data.nnr)+Data.femorder(2),[]);
Data.fluidvz=permute(Data.fluidvz,[2,1,3]);
%Data.Presstens=zeroes(6,size(Data.zgrid,1),size(Data.rgrid,1), size(moments,3));
%Data.Presstens(1,:,:,:)=reshape(moments(5,:,:),Data.nz+Data.femorder(1),sum(Data.nnr)+Data.femorder(2),[])
%Data.Presstens=Data.me*Data.vlight^2*reshape(h5read(file,'/data/fields/presstens'),6,length(Data.zgrid)-1,length(Data.rgrid)-1,[]);
%Data.Presstens=permute(Data.Presstens,[1,3,2,4]);
clear moments;
Data.N=zeros(length(Data.rgrid),length(Data.zgrid),size(Data.dens,3));
for i=1:size(Data.dens,3)
Data.dens(:,:,i)=Data.dens(:,:,i).*Data.invV*abs(Data.qsim/Data.qe);
Data.N(:,:,i)=fnval(spmak({Data.knotsr,Data.knotsz},Data.dens(:,:,i)),{Data.rgrid,Data.zgrid});
end
%Data.moments=reshape(Data.moments,10,length(Data.zgrid),length(Data.rgrid),[]);
%Data.moments=permute(Data.moments,[1,3,2,4]);
catch
N = h5read(file,'/data/fields/partdensity');
Data.N=reshape(N,length(Data.zgrid)-1,length(Data.rgrid)-1,[]);
Data.N=permute(Data.N,[2,1,3]);
Data.Npartscell=Data.N;
Vcell=(Data.zgrid(2:end)-Data.zgrid(1:end-1))*((Data.rgrid(2:end).^2-Data.rgrid(1:end-1).^2)*pi)';
for i=1:size(Data.N,3)
Data.N(:,:,i)=abs(Data.qsim/Data.qe)*Data.N(:,:,i)./Vcell';
end
Data.N=padarray(Data.N,[1,1],0,'post');
Data.Npartscell=padarray(Data.Npartscell,[1,1],0,'post');
clear Vcell
clear N
Data.fluidUR=reshape(h5read(file,'/data/fields/fluidur'),length(Data.zgrid)-1,length(Data.rgrid)-1,[]);
Data.fluidUR=permute(Data.fluidUR,[2,1,3])*Data.vnorm;
Data.fluidUR=padarray(Data.fluidUR,[1,1],0,'post');
Data.fluidUZ=reshape(h5read(file,'/data/fields/fluiduz'),length(Data.zgrid)-1,length(Data.rgrid)-1,[]);
Data.fluidUZ=permute(Data.fluidUZ,[2,1,3])*Data.vnorm;
Data.fluidUZ=padarray(Data.fluidUZ,[1,1],0,'post');
Data.fluidUTHET=reshape(h5read(file,'/data/fields/fluiduthet'),length(Data.zgrid)-1,length(Data.rgrid)-1,[]);
Data.fluidUTHET=permute(Data.fluidUTHET,[2,1,3])*Data.vnorm;
Data.fluidUTHET=padarray(Data.fluidUTHET,[1,1],0,'post');
Data.Presstens=Data.me*Data.vlight^2*reshape(h5read(file,'/data/fields/presstens'),6,length(Data.zgrid)-1,length(Data.rgrid)-1,[]);
Data.Presstens=permute(Data.Presstens,[1,3,2,4]);
Data.Presstens=padarray(Data.Presstens,[0,1,1,0],0,'post');
for i=1:6
Data.Presstens(i,:,:,:)=squeeze(Data.Presstens(i,:,:,:)).*Data.N;
end
end
try
Data.t2d = h5read(file,'/data/fields/time');
catch
Data.t2d=Data.dt*(0:size(Data.N,3)-1)*double(Data.it1);
end
if(strcmp(Data.type,'parts'))
Data.type='all';
end
end
if(strcmp(type,'all') || strcmp(type,'parts'))
Data.R = h5read(file,'/data/part/R')*Data.rnorm;
partsbegin=floor(partstime(1)*size(Data.R,2))+1;
partsend=floor(partstime(2)*size(Data.R,2));
Data.R=Data.R(:,partsbegin:partsend);
Data.Z = h5read(file,'/data/part/Z',[1,partsbegin],[size(Data.R,1) partsend])*Data.rnorm;
try
Data.THET = h5read(file,'/data/part/THET',[1,partsbegin],[size(Data.R,1) partsend]);
catch
clear Data.THET
end
try
Data.Rindex=h5read(file,'/data/part/Rindex',[1,partsbegin],[size(Data.R,1) partsend]);
Data.Zindex=h5read(file,'/data/part/Zindex',[1,partsbegin],[size(Data.R,1) partsend]);
catch
clear Data.Rindex Data.Zindex
end
UR = h5read(file,'/data/part/UR',[1,partsbegin],[size(Data.R,1) partsend]);
UZ = h5read(file,'/data/part/UZ',[1,partsbegin],[size(Data.R,1) partsend]);
UTHET= h5read(file,'/data/part/UTHET',[1,partsbegin],[size(Data.R,1) partsend]);
if(Data.nlclassical)
Data.VR = UR*Data.vnorm;
clear UR
Data.VZ = UZ*Data.vnorm;
clear UZ
Data.VTHET = UTHET*Data.vnorm;
clear UTHET
else
GAMM = sqrt(1+UR.^2+UZ.^2+UTHET.^2);
Data.VR = UR./GAMM*Data.vnorm;
clear UR
Data.VZ = UZ./GAMM*Data.vnorm;
clear UZ
Data.VTHET = UTHET./GAMM*Data.vnorm;
clear UTHET
clear GAMM
end
%Data.VPERP = sqrt(Data.VR.*Data.VR+Data.VTHET.*Data.VTHET);
partepot = sign(Data.qsim)*h5read(file,'/data/part/pot',[1,partsbegin],[size(Data.R,1) partsend])*Data.phinorm*Data.qe;
Data.H=0.5*Data.me*(Data.VR.^2+Data.VTHET.^2+Data.VZ.^2)+partepot;
clear partepot
Data.P=Data.R.*(Data.VTHET*Data.me+sign(Data.qsim)*Data.qe*Athet(Data.R,Data.Z));
try
partindex = h5read(file,'/data/part/partindex',[1,partsbegin],[size(Data.R,1) partsend]);
[~,Indices]=sort(partindex);
for i=1:size(Indices,2)
Data.H(:,i)=Data.H(Indices(:,i),i);
Data.P(:,i)=Data.P(Indices(:,i),i);
Data.R(:,i)=Data.R(Indices(:,i),i);
Data.Z(:,i)=Data.Z(Indices(:,i),i);
Data.Rindex(:,i)=Data.Rindex(Indices(:,i),i);
Data.Zindex(:,i)=Data.Zindex(Indices(:,i),i);
Data.VR(:,i)=Data.VR(Indices(:,i),i);
Data.VZ(:,i)=Data.VZ(Indices(:,i),i);
Data.VTHET(:,i)=Data.VTHET(Indices(:,i),i);
Data.THET(:,i)=Data.THET(Indices(:,i),i);
end
catch
end
clear partindex;
try
Data.tpart = h5read(file,'/data/part/time');
catch
Data.tpart=Data.dt*(0:size(Data.R,2)-1)*double(Data.it2);
end
if(strcmp(Data.type,'fields'))
Data.type='all';
end
end
Data.t0d=Data.dt.*double(0:length(Data.epot)-1);
filedata=dir(file);
Data.timestamp=filedata.date;
disp('M reloaded')
function Atheta=Athet(R,Z)
Atheta=0.5*Data.B0*(R-Data.width/pi*(Data.Rcurv-1)/(Data.Rcurv+1)...
.*besseli(1,2*pi*R/Data.width).*cos(2*pi*Z/Data.width));
end
end
c4science · Help