From 2770168e0ede3b9eabda784878d8dfdd4d0134ef Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Mon, 3 Jul 2017 10:32:15 +0000
Subject: [PATCH] fix phi_tor since not yet ready with liuqe.m
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@7733 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/TCV/gdat_tcv.m | 81 +++++++++++++++++++++++++-----
crpptbx/TCV/tcv_requests_mapping.m | 5 +-
2 files changed, 72 insertions(+), 14 deletions(-)
diff --git a/crpptbx/TCV/gdat_tcv.m b/crpptbx/TCV/gdat_tcv.m
index af483bcf..2d9cb68e 100644
--- a/crpptbx/TCV/gdat_tcv.m
+++ b/crpptbx/TCV/gdat_tcv.m
@@ -281,6 +281,17 @@ if shot==-1 || (shot>=100000 && shot < 200000) || liuqe_version==-1
substr_liuqe = '", "FBTE" )';
end
+% should replace all above by just psitbx_str...
+liuqe_matlab = 1;
+switch liuqe_version
+ case {-1}, liuqe_ext=''; psitbx_str='FBTE'; liuqe_matlab = 0;
+ case {1,21}, liuqe_ext=''; psitbx_str='LIUQE.M';
+ case {11}, liuqe_ext=''; psitbx_str='LIUQE';liuqe_matlab = 0;
+ case {2, 3, 22, 23}, liuqe_ext=['_' num2str(mod(liuqe_version,10))]; psitbx_str=['LIUQE.M' num2str(mod(liuqe_version,10))];
+ case {12,13}, liuqe_ext=['_' num2str(mod(liuqe_version,10))]; psitbx_str=['LIUQE' num2str(mod(liuqe_version,10))];liuqe_matlab = 0;
+ otherwise, error(['Unknown LIUQE version, liuqe = ' num2str(liuqe_version)]);
+end
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specifications on how to get the data provided in tcv_requests_mapping
@@ -692,11 +703,11 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
gdat_data.dim = beta.dim;
gdat_data.t = beta.dim{1};
gdat_data.data = beta.data;
- ij=find(~isnan(ip.data));
+ ij=find(isfinite(ip.data));
ip_t = interp1(ip.dim{1}(ij),ip.data(ij),gdat_data.t);
- ij=find(~isnan(b0.data));
+ ij=find(isfinite(b0.data));
b0_t = interp1(b0.dim{1}(ij),b0.data(ij),gdat_data.t);
- ij=find(~isnan(a_minor.data));
+ ij=find(isfinite(a_minor.data));
a_minor_t = interp1(a_minor.dim{1}(ij),a_minor.data(ij),gdat_data.t);
gdat_data.data = 100.*beta.data ./ abs(ip_t).*1.e6 .* abs(b0_t) .* a_minor_t;
gdat_data.data_fullpath='100*beta/ip*1e6*b0*a_minor, each from gdat_tcv';
@@ -1464,6 +1475,42 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
% add grids_1d to have rhotor, etc
gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose);
+ case {'phi_tor', 'phitor', 'toroidal_flux'}
+ % Phi(LCFS) = int(Bphi dSphi), can use Eq.(11) of "Tokamak coordinate conventions: COCOS" paper:
+ % O. Sauter, S.Yu. Medvedev, Comput. Phys. Commun. 184 (2013) 293–302
+ % since cocos=17 for LIUQE we get:
+ % q = -dPhi/dpsi => Phi = - int(q*dpsi) which should always have the sign of B0
+ % need to get q_rho but to avoid loop for rhotor in grids_1d, get q_rho explicitely here
+ params_eff = gdat_data.gdat_params;
+ if liuqe_matlab==1
+ nodenameeff = 'tcv_eq(''tor_flux_tot'',''LIUQE.M'')';
+ phi_tor = tdi(nodenameeff);
+ else
+ phi_tor.data = [];
+ phi_tor.units = 'Wb';
+ end
+ if ~any(any(isfinite(phi_tor.data)))
+ q_rho = tdi(['tcv_eq(''''''' liuqefortran2liuqematlab('q',1,liuqe_matlab) ''''''',''''''' psitbx_str ''''''')']);
+ if liuqe_matlab==0
+ q_rho.dim{1} = sqrt(linspace(0,1,numel(q_rho.dim{1})));
+ end
+ phi_tor.dim = q_rho.dim;
+ phi_tor.dimunits = q_rho.dimunits;
+ for it=1:size(q_rho.data,2)
+ ij=find(isfinite(q_rho.data(:,it)));
+ if ~isempty(ij) && numel(ij)>5
+ [qfit,~,~,phi]=interpos(q_rho.dim{1}(ij),q_rho.data(ij,it),phi_tor.dim{1});
+ phi_tor.data(:,it) = phi';
+ end
+ end
+ end
+ gdat_data.data = phi_tor.data;
+ gdat_data.dim = phi_tor.dim;
+ gdat_data.dimunits = phi_tor.dimunits;
+ gdat_data.units = phi_tor.units;
+ gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim};
+ gdat_data.x = gdat_data.dim{1};
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'psi_edge'}
% psi at edge, 0 by construction in Liuqe, thus not given
@@ -1473,7 +1520,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
nodenameeff=[begstr 'q_psi' substr_liuqe];
end
tracetdi=tdi(nodenameeff);
- if isempty(tracetdi.data) || isempty(tracetdi.dim) || ~any(~isnan(tracetdi.data)) % || ischar(tracetdi.data) (to add?)
+ if isempty(tracetdi.data) || isempty(tracetdi.dim) || ~any(isfinite(tracetdi.data)) % || ischar(tracetdi.data) (to add?)
if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
return
@@ -1550,7 +1597,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
return
end
for it=1:size(q_rho.data,2)
- ij=find(~isnan(q_rho.data(:,it)));
+ ij=find(isfinite(q_rho.data(:,it)));
if ~isempty(ij)
[qfit,~,~,phi]=interpos(q_rho.dim{1}(ij).^2,q_rho.data(ij,it),rhoequal.^2);
dataeff = sqrt(phi .* psi_axis.data(it) ./ b0tpsi(it) ./ pi) ; % Delta_psi = -psi_axis
@@ -1569,11 +1616,22 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
params_eff = gdat_data.gdat_params;
params_eff.data_request='phi_tor';
phi_tor=gdat_tcv([],params_eff);
+% $$$ if ~any(any(isfinite(phi_tor.data)))
+% $$$ % not yet implemented, compute phi_tor
+% $$$ q_rho = tdi(['tcv_eq(''''''' liuqefortran2liuqematlab('q',1,liuqe_matlab) ''''''',''''''' psitbx_str ''''''')']);
+% $$$ for it=1:size(q_rho.data,2)
+% $$$ ij=find(isfinite(q_rho.data(:,it)));
+% $$$ if ~isempty(ij) && numel(ij)>5
+% $$$ [qfit,~,~,phi]=interpos(q_rho.dim{1}(ij),q_rho.data(ij,it),phi_tor.dim{1});
+% $$$ phi_tor.data(:,it) = phi;
+% $$$ end
+% $$$ end
+% $$$ end
params_eff = gdat_data.gdat_params;
params_eff.data_request='b0';
b0=gdat_tcv([],params_eff);
gdat_data.t = phi_tor.t;
- ij=find(~isnan(b0.data));
+ ij=find(isfinite(b0.data));
gdat_data.b0 = interpos(b0.t(ij),b0.data(ij),gdat_data.t);
if strcmp(data_request_eff,'rhotor_edge')
gdat_data.data = sqrt(phi_tor.data(end,:)./pi./reshape(gdat_data.b0,1,numel(gdat_data.b0)));
@@ -1611,7 +1669,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
if liuqe_matlab==0
nodenameeff=['\results::psitbx:vol'];
if liuqe_version_eff > 1
- disp('needs to construct volume')
+ disp('needs to construct volume');
return
end
else
@@ -1667,7 +1725,6 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
return
end
end
- mapping_for_tcv.gdat_timedim
gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -2118,7 +2175,7 @@ else
return
end
-if ~isempty(tracefirrat.data) && ~ischar(tracefirrat.data) && any(~isnan(tracefirrat.data)) ...
+if ~isempty(tracefirrat.data) && ~ischar(tracefirrat.data) && any(isfinite(tracefirrat.data)) ...
&& ~isempty(tracefirrat.dim) && ~isempty(tracefirrat.dim{1})
firthomratio = interp1(tracefirrat.dim{1},tracefirrat.data,timebase);
end
@@ -2181,7 +2238,7 @@ if strcmp(data_request_eff(1:2),'ne')
gdat_data.error_bar_raw = gdat_data.error_bar;
gdat_data.error_bar = gdat_data.error_bar_raw * diag(tracefirrat_data);
gdat_data.data_fullpath=[gdat_data.data_fullpath '; fir_thom_rat ; _raw without firrat'];
- ij=find(~isnan(tracefirrat_data));
+ ij=find(isfinite(tracefirrat_data));
if isempty(ij)
gdat_data.data_fullpath=[gdat_data.data_fullpath ' WARNING use ne_raw since firrat has NaNs thus ne=NaNs; fir_thom_rat ; _raw without firrat'];
disp('***********************************************************************')
@@ -2246,9 +2303,9 @@ for it=1:length(gdat_data.t)
it_rt_eff = it_rt(it);
it_vol_eff = it_vol(it);
if (nbdim_x == 1)
- ii=find(~isnan(gdat_data.grids_1d.rhopolnorm));
+ ii=find(isfinite(gdat_data.grids_1d.rhopolnorm));
else
- ii=find(~isnan(gdat_data.grids_1d.rhopolnorm(:,it)));
+ ii=find(isfinite(gdat_data.grids_1d.rhopolnorm(:,it)));
end
if (nbdim_x == 1)
if length(ii)==length(gdat_data.grids_1d.rhopolnorm)
diff --git a/crpptbx/TCV/tcv_requests_mapping.m b/crpptbx/TCV/tcv_requests_mapping.m
index 07d4e746..b66cef30 100644
--- a/crpptbx/TCV/tcv_requests_mapping.m
+++ b/crpptbx/TCV/tcv_requests_mapping.m
@@ -219,9 +219,10 @@ switch lower(data_request)
case {'phi_tor', 'phitor', 'toroidal_flux'}
mapping.timedim = 2;
mapping.method = 'tdiliuqe';
- mapping.expression = '\results::psi_axis';
- mapping.expression = '\tcv_shot::top.results.equil_1.results:psi_axis';
mapping.expression = 'tcv_eq(''''tor_flux_tot'''',''''LIUQE.M'''')';
+ % node not filled in and in any case need special case for
+ mapping.label = 'toroidal_flux';
+ mapping.method = 'switchcase';
case 'powers'
mapping.timedim = 1;
mapping.label = 'various powers';
--
GitLab