Newer
Older
function [L,LX,LY] = LLXLY_IDS_corsica(L_LIU, a)
%% Load data drom IDS
% Get geomtrical information from the existing LIUQE
L = struct();
L.P = L_LIU.P;
fieldlist = {'rl';'zl'; 'rw';'zw'; 'Twa'; 'rv';'zv'; 'Tvs'; 'rf'; 'zf';'rm'; 'zm'; 'am'};
for ii= 1:numel(fieldlist)
tmpname = fieldlist{ii};
L.G.(tmpname) = L_LIU.G.(tmpname);
end
% Set new grid
L.G.rx = a.equilibrium.time_slice{1}.profiles_2d{1}.grid.dim1;
L.G.zx = a.equilibrium.time_slice{1}.profiles_2d{1}.grid.dim2;
L.P.pq = sqrt((a.equilibrium.time_slice{1}.profiles_1d.psi - a.equilibrium.time_slice{1}.global_quantities.psi_axis)/(a.equilibrium.time_slice{1}.global_quantities.psi_boundary - a.equilibrium.time_slice{1}.global_quantities.psi_axis));
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
% Recompute the geometrical quantities for the CORSICA grid
L.G = liug(L.G,L.P);
L = liuc(L.P, L.G);
[rrx, zzx] = meshgrid(L.G.rx, L.G.zx);
%% Get LX structure
% {Ia, Ip, rBt} available
% Ff, Bm, Ft recomputed
LX.t = a.pf_active.coil{1}.current.time;
LX.nt = numel(LX.t);
LX.Ip = zeros(LX.nt,1);
for ii=1:LX.nt
LX.Ip(ii) = a.equilibrium.time_slice{ii}.global_quantities.ip;
end
% Ia is the current per number of turns
LX.Ia = zeros(L.G.na, LX.nt);
for ii= 1:L.G.na
LX.Ia(ii,:) = a.pf_active.coil{ii}.current.data/sum(L.G.Twa(:,ii));
end
LX.rBt = zeros(LX.nt,1);
LX.rBt(:) = a.equilibrium.vacuum_toroidal_field.r0*a.equilibrium.vacuum_toroidal_field.b0;
%% Synthetic data that needs to be recomputed
LY.Ix = zeros( L.nzx, L.nrx ,LX.nt);
for tt=1:LX.nt
LY.Ix(:,:,tt) = a.equilibrium.time_slice{tt}.profiles_2d{1}.j_tor./rrx*L.drx*L.dzx*1e-2; % Wrong definition in IDS Corsica
% Ix_Corsica(:,:,tt) = a.equilibrium.time_slice{tt}.profiles_2d{1}.j_tor; % Wrong definition of CORSICA
% Ix_Corsica(:,:,tt) = -Ix_Corsica(:,:,tt)*abs(LX.Ip(tt)/sum(sum(Ix_Corsica(:,:,tt))));
end
LX.Bm = zeros(L.G.nm, LX.nt);
for tt= 1:LX.nt
LX.Bm(:,tt) = L.G.Bmx*reshape(LY.Ix(:,:,tt), L.nrx*L.nzx,1 ) + L.G.Bma*LX.Ia(:,tt);
end
LX.Ff = zeros(L.G.nf, LX.nt);
for tt= 1:LX.nt
LX.Ff(:,tt) = L.G.Mfx*reshape(LY.Ix(:,:,tt), L.nrx*L.nzx,1 ) + L.G.Mfa*LX.Ia(:,tt);
end
LX.Is = zeros(0,LX.nt);
LX.shot = 5000;
LY.Fx = zeros(L.nzx,L.nrx ,LX.nt);
for tt = 1:LX.nt
LY.Fx(:,:,tt) = a.equilibrium.time_slice{tt}.profiles_2d{1}.psi;
end
% Normalized flux of the R,Z grid
LY.Opy = zeros(L.nzy, L.nry, LX.nt);
LY.FNy = zeros(L.nzy, L.nry, LX.nt);
for tt=1:LX.nt
[LY.rA(tt),LY.zA(tt),LY.FA(tt),dr2FA,dz2FA,drzFA,rX,zX,~,~,~,~,rl,zl,...
LY.rB(tt),LY.zB(tt),LY.FB(tt),LY.lB(tt),LY.lX(tt),LY.Opy(:,:,tt)] = meqpdom(LY.Fx(:,:,tt),LX.Ip(tt),L.P.isaddl,L);
LY.FNy(:,:,tt) = (LY.Fx(2:end-1,2:end-1,tt) - LY.FA(tt))/(LY.FB(tt) - LY.FA(tt));
end
%% Recompute the plasma current as a check. There is a bit of a difference
psiN = (a.equilibrium.time_slice{tt}.profiles_1d.psi - a.equilibrium.time_slice{tt}.global_quantities.psi_axis)/(a.equilibrium.time_slice{tt}.global_quantities.psi_boundary - a.equilibrium.time_slice{tt}.global_quantities.psi_axis);
for tt=1:LX.nt
pprime = a.equilibrium.time_slice{tt}.profiles_1d.dpressure_dpsi;
p = spline(psiN,pprime);
pprimey = ppval(p, LY.FNy(:,:,tt));
pprimey = pprimey.*LY.Opy(:,:,tt);
ttprime = a.equilibrium.time_slice{tt}.profiles_1d.f_df_dpsi;
T2d2 = cumtrapz(a.equilibrium.time_slice{tt}.profiles_1d.psi,ttprime);
T2d2 = T2d2 - T2d2(end);
p = spline(psiN,ttprime);
ttprimey = ppval(p, LY.FNy(:,:,tt));
ttprimey = ttprimey.*LY.Opy(:,:,tt);
p = spline(psiN,T2d2);
Ty = ppval(p, LY.FNy(:,:,tt));
Ty = Ty.*LY.Opy(:,:,tt);
L.rymap = repmat(L.ry',numel(L.zy),1); % Matriz map ry
LX.Ft(tt) = sum(sum( Ty/LX.rBt(tt)*L.drx*L.dzx./L.rymap ));
Iy = 2*pi*(L.rymap.*pprimey+1./L.rymap.*ttprimey/mu0)*L.drx*L.dzx; % Recomputed plasma current
LY.Ip(tt) = sum(sum(Iy));
LY.FA(tt) = a.equilibrium.time_slice{tt}.global_quantities.psi_axis;
LY.FB(tt) = a.equilibrium.time_slice{tt}.global_quantities.psi_boundary;
LY.Wk(tt) = a.equilibrium.time_slice{tt}.global_quantities.w_mhd;
LY.PpQ(:,tt) = pprime;
LY.TTpQ(:,tt) = ttprime;
LY.psiN(:,tt) = psiN;
LY.PQ(:,tt) = a.equilibrium.time_slice{tt}.profiles_1d.pressure;
end
LY.t = LX.t;
% Put all vector in raw form
flist = fieldnames(LX);
for field = flist
LX.(field{1}) = v2rw(LX.(field{1}));
end
% Put all vector in raw form
flist = fieldnames(LY);
for field = flist
LY.(field{1}) = v2rw(LY.(field{1}));
end
end
% Put all vectors in raw form
function vout = v2rw(vin)
if numel(size(vin))==2 && sum(ismember(1,size(vin)))==1 && ~ischar(vin)
vout = reshape(vin, numel(vin),1);
else
vout = vin;
end
end