Skip to content
Snippets Groups Projects
Commit a96c27fe authored by Olivier Sauter's avatar Olivier Sauter
Browse files

fix present problem of EQI, EQR, etc from within loadAUG (time base not...

fix present problem of EQI, EQR, etc from within loadAUG (time base not consistent) with a method which should be OK when OK, taking NTIME first points

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@4468 d63d8f72-b253-0410-a779-e742ad2e26cf
parent 2acfa533
Branches
Tags
No related merge requests found
......@@ -460,10 +460,23 @@ if size(data_type_eff,1)==2
AUGsiglocation(1:2,end+1)=[data_type_eff(1) ; {data_type_eff_noext}];
AUGexplocation{end+1}=shotfile_exp;
AUGsigtimeindx(end+1)=0;
elseif ~strcmp(AUGkeywrdcase{index},'simplereaddata')
msgbox(['Problem in loadAUGdata with data_type_eff = ' char(data_type_eff(end)) ...
'. Full paths of nodes should only be for case simplereaddata'],'in loadAUGdata','error')
error('in loadAUGdata')
AUGkeywrdcase{index}
elseif length(index)>1 || ~strcmp(AUGkeywrdcase{index},'simplereaddata')
% $$$ msgbox(['Problem in loadAUGdata with data_type_eff = ' char(data_type_eff(end)) ...
% $$$ '. Full paths of nodes should only be for case simplereaddata'],'in loadAUGdata','error')
% $$$ error('in loadAUGdata')
% chose a full path which coincides with a keyword
disp('the node path seems to coincide with the following keywords, which might be better to use:')
disp(' ')
disp(AUGkeywrdcase(index))
disp('proceed as is in any case')
index=length(AUGkeywrdall)+1;
AUGkeywrdall(end+1)={'new'};
AUGkeywrdcase(end+1)={'simplereaddata'};
AUGsiglocation(1:2,end+1)=[data_type_eff(1) ; {data_type_eff_noext}];
AUGexplocation{end+1}=shotfile_exp;
AUGsigtimeindx(end+1)=0;
AUGkeywrdcase{index}
end
else
index=strmatch(data_type_eff_noext,AUGkeywrdall,'exact');
......@@ -1059,6 +1072,9 @@ switch AUGkeywrdcase{index}
case {'equil', 'equil_fpp', 'equil_eqm', 'equil_eqr', 'equil_eqh', 'qrho', 'qrho_fpp'}
% In shotfiles time is 1st index in 2D arrays while gdat convention is to have time as last dimension
% so change equil.xxx to (rho,time) on the way...
shotfile_exp_eff = AUGexplocation{index};
if strcmp(AUGkeywrdcase{index},'equil_fpp') || strcmp(AUGkeywrdcase{index},'qrho_fpp')
......@@ -1079,90 +1095,110 @@ switch AUGkeywrdcase{index}
NTIME_par = sf2par(DIAG,shot,'NTIME','PARMV');
NTIME = NTIME_par.value; % nb of points
Lpf_par = rdaAUG_eff(shot,DIAG,'Lpf',shotfile_exp_eff);
% since June, nb of time points in EQ results is not consistent with NTIME and time
% It seems the first NTIME points are correct, so use this explicitely
NTIME_Lpf = length(Lpf_par.value);
Lpf_tot = Lpf_par.value; % nb of points: 100000*nb_SOL points + nb_core
if (NTIME < NTIME_Lpf)
disp('WARNING: nb of times points smaller then equil results, use first NTIME points')
elseif (NTIME > NTIME_Lpf)
disp('ERROR: nb of times points LARGER then equil results')
disp('this is unexpected, so stop there and ask Olivier.Sauter@epfl.ch')
return
end
Lpf_tot = Lpf_par.value(1:NTIME); % nb of points: 100000*nb_SOL points + nb_core
Lpf_SOL = fix(Lpf_tot/100000);
Lpf1_t = mod(Lpf_tot,100000)+1; % nb of Lpf points
% since Lpf depends on time, need to load all first and then loop over time for easier mapping
% since Lpf depends on time, need to load all first and then loop over time for easier mapping
[qpsi,e]=rdaAUG_eff(shot,DIAG,'Qpsi',shotfile_exp_eff);
ndimrho = size(qpsi.data,2);
if ndimrho==NTIME_Lpf
% data seems to be transposed
ndimrho = size(qpsi.data,1);
itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
else
itotransposeback = 0;
end
qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback);
[psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',shotfile_exp);
psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
[phi_tree,e]=rdaAUG_eff(shot,DIAG,'TFLx',shotfile_exp);
phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback);
[Vol,e]=rdaAUG_eff(shot,DIAG,'Vol',shotfile_exp);
Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback);
[Area,e]=rdaAUG_eff(shot,DIAG,'Area',shotfile_exp);
Area=adapt_rda(Area,NTIME,2*ndimrho,itotransposeback);
[Ri,e]=rdaAUG_eff(shot,DIAG,'Ri',shotfile_exp);
Ri=adapt_rda(Ri,NTIME,M_Rmesh,itotransposeback);
[Zj,e]=rdaAUG_eff(shot,DIAG,'Zj',shotfile_exp);
Zj=adapt_rda(Zj,NTIME,N_Zmesh,itotransposeback);
[PFM_tree,e]=rdaAUG_eff(shot,DIAG,'PFM',shotfile_exp);
PFM_tree=adaptPFM_rda(PFM_tree,M_Rmesh,N_Zmesh,NTIME);
[Pres,e]=rdaAUG_eff(shot,DIAG,'Pres',shotfile_exp);
Pres=adapt_rda(Pres,NTIME,2*ndimrho,itotransposeback);
[FFP,e]=rdaAUG_eff(shot,DIAG,'FFP',shotfile_exp);
if isempty(FFP.value); FFP.value=NaN*ones(NTIME,max(Lpf1_t)); end
if ~isempty(FFP.value)
FFP=adapt_rda(FFP,NTIME,ndimrho,itotransposeback);
else
FFP.value=NaN*ones(NTIME,max(Lpf1_t));
end
[LPFx,e]=rdaAUG_eff(shot,DIAG,'LPFx',shotfile_exp);
LPFx.value=LPFx.value(1:NTIME); LPFx.data=LPFx.value; LPFx.t=LPFx.t(1:NTIME);
[PFxx,e]=rdaAUG_eff(shot,DIAG,'PFxx',shotfile_exp);
PFxx=adapt_rda(PFxx,NTIME,max(LPFx.value)+1,itotransposeback);
[RPFx,e]=rdaAUG_eff(shot,DIAG,'RPFx',shotfile_exp);
RPFx=adapt_rda(RPFx,NTIME,max(LPFx.value)+1,itotransposeback);
[zPFx,e]=rdaAUG_eff(shot,DIAG,'zPFx',shotfile_exp);
zPFx=adapt_rda(zPFx,NTIME,max(LPFx.value)+1,itotransposeback);
% seems "LCFS" q-value is far too large, limit to some max (when diverted)
max_qValue = 40; % Note could just put a NaN on LCFS value since ill-defined when diverted
% check NTIME ok with qvalue
if NTIME== size(qpsi.value,2)
for it=1:NTIME
Lpf1 = Lpf1_t(it);
% Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part
% change it to (radial,time) and use only Lpf+1 points up to LCFS
if qpsi.value(1,2)<0
equil.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
max_qValue = 40.0; % Note could just put a NaN on LCFS value since ill-defined when diverted
for it=1:NTIME
Lpf1 = Lpf1_t(it);
% Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part
% change it to (radial,time) and use only Lpf+1 points up to LCFS
if qpsi.value(1,2)<0
equil.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
else
equil.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
end
% get x values
equil.psi(:,it)=psi_tree.value(it,Lpf1:-1:1)';
equil.psi_axis(it)= equil.psi(1,it);
equil.psi_lcfs(it)= equil.psi(end,it);
equil.rhopolnorm(:,it) = sqrt(abs((equil.psi(:,it)-equil.psi_axis(it)) ./(equil.psi_lcfs(it)-equil.psi_axis(it))));
if strcmp(DIAG,'EQR');
% q value has only a few values and from center to edge, assume they are from central rhopol values on
% But they are every other point starting from 3rd
ijk=find(equil.qvalue(:,it)~=0);
qeff=equil.qvalue(ijk(end:-1:1),it);
rhoeff=equil.rhopolnorm(3:2:2*length(ijk)+1,it);
if length(ijk)>2
ij_nonan=find(~isnan(equil.rhopolnorm(:,it)));
qfit = zeros(size(equil.rhopolnorm(:,it)));
qfit(ij_nonan)=interpos([0.; rhoeff],[qeff(1) ;qeff],equil.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff))]);
else
equil.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
qfit = zeros(size(equil.rhopolnorm(:,it)));
end
% get x values
equil.psi(:,it)=psi_tree.value(it,Lpf1:-1:1)';
equil.psi_axis(it)= equil.psi(1,it);
equil.psi_lcfs(it)= equil.psi(end,it);
equil.rhopolnorm(:,it) = sqrt(abs((equil.psi(:,it)-equil.psi_axis(it)) ./(equil.psi_lcfs(it)-equil.psi_axis(it))));
if strcmp(DIAG,'EQR');
% q value has only a few values and from center to edge, assume they are from central rhopol values on
% But they are every other point starting from 3rd
ijk=find(equil.qvalue(:,it)~=0);
qeff=equil.qvalue(ijk(end:-1:1),it);
rhoeff=equil.rhopolnorm(3:2:2*length(ijk)+1,it);
if length(ijk)>2
ij_nonan=find(~isnan(equil.rhopolnorm(:,it)));
qfit = zeros(size(equil.rhopolnorm(:,it)));
qfit(ij_nonan)=interpos([0.; rhoeff],[qeff(1) ;qeff],equil.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff))]);
else
qfit = zeros(size(equil.rhopolnorm(:,it)));
end
equil.qvalue(:,it) = qfit;
end
% get rhotor values
equil.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)';
equil.rhotornorm(:,it) = sqrt(abs(equil.phi(:,it) ./ equil.phi(end,it)));
% get rhovol values
equil.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)';
equil.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
equil.rhovolnorm(:,it) = sqrt(abs(equil.vol(:,it) ./ equil.vol(end,it)));
equil.area(:,it)=Area.value(it,2*Lpf1-1:-2:1)';
equil.dareadpsi(:,it)=Area.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
equil.Rmesh(:,it) = Ri.value(it,1:M_Rmesh);
equil.Zmesh(:,it) = Zj.value(it,1:N_Zmesh);
equil.psi2D(1:M_Rmesh,1:N_Zmesh,it) = PFM_tree.value(1:M_Rmesh,1:N_Zmesh,it);
equil.pressure(:,it)=Pres.value(it,2*Lpf1-1:-2:1)';
equil.dpressuredpsi(:,it)=Pres.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
equil.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)';
equil.Xpoints.psi(1:LPFx.value(it)+1,it) = PFxx.value(it,1:LPFx.value(it)+1);
equil.Xpoints.Rvalue(1:LPFx.value(it)+1,it) = RPFx.value(it,1:LPFx.value(it)+1);
equil.Xpoints.Zvalue(1:LPFx.value(it)+1,it) = zPFx.value(it,1:LPFx.value(it)+1);
%
end
else
disp('***********************************************************************')
disp('***********************************************************************')
disp(['problem with nb of time points: NTIME for parameter is: ' num2str(NTIME) ' while qpsi dim= ' num2str(size(qpsi.value))])
disp('***********************************************************************')
disp('***********************************************************************')
equil.qvalue = [];
equil.qpsi=qpsi;
equil.PFL=psi_tree;
equil.TFLx=phi_tree;
equil.rhopolnorm=[];
equil.qvalue(:,it) = qfit;
end
% get rhotor values
equil.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)';
equil.rhotornorm(:,it) = sqrt(abs(equil.phi(:,it) ./ equil.phi(end,it)));
% get rhovol values
equil.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)';
equil.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
equil.rhovolnorm(:,it) = sqrt(abs(equil.vol(:,it) ./ equil.vol(end,it)));
equil.area(:,it)=Area.value(it,2*Lpf1-1:-2:1)';
equil.dareadpsi(:,it)=Area.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
equil.Rmesh(:,it) = Ri.value(it,1:M_Rmesh);
equil.Zmesh(:,it) = Zj.value(it,1:N_Zmesh);
equil.psi2D(1:M_Rmesh,1:N_Zmesh,it) = PFM_tree.value(1:M_Rmesh,1:N_Zmesh,it);
equil.pressure(:,it)=Pres.value(it,2*Lpf1-1:-2:1)';
equil.dpressuredpsi(:,it)=Pres.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
equil.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)';
equil.Xpoints.psi(1:LPFx.value(it)+1,it) = PFxx.value(it,1:LPFx.value(it)+1);
equil.Xpoints.Rvalue(1:LPFx.value(it)+1,it) = RPFx.value(it,1:LPFx.value(it)+1);
equil.Xpoints.Zvalue(1:LPFx.value(it)+1,it) = zPFx.value(it,1:LPFx.value(it)+1);
%
end
equil.x = equil.rhopolnorm;
%
......@@ -1175,7 +1211,7 @@ switch AUGkeywrdcase{index}
equil.t = equil_time.value(1:NTIME);
%
[IpiPSI,e]=rdaAUG_eff(shot,DIAG,'IpiPSI',shotfile_exp);
equil.Ip = IpiPSI.value';
equil.Ip = IpiPSI.value(1:NTIME);
%
equil.data = equil.qvalue; % put q in data
equil.unit=[]; % not applicable
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment