From 4dffbd0fd036b0c6b7081c87062c72476d091a6e Mon Sep 17 00:00:00 2001
From: Antoine Merle <antoine.merle@epfl.ch>
Date: Fri, 24 Jan 2020 12:46:34 +0100
Subject: [PATCH] Add equilibrium constraints when fetching LIUQE data.
---
matlab/TCV_IMAS/tcv_get_ids_equilibrium.m | 157 ++++++++++++++++++----
matlab/TCV_IMAS/tcv_get_ids_magnetics.m | 2 +-
2 files changed, 129 insertions(+), 30 deletions(-)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
index 391d8e44..bb275385 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
@@ -1,4 +1,4 @@
-function [ids_equilibrium,ids_equilibrium_description,varargout] = tcv_get_ids_equilibrium(shot,ids_equil_empty, gdat_params,varargin);
+function [ids_equilibrium,ids_equilibrium_description,varargout] = tcv_get_ids_equilibrium(shot,ids_equil_empty, gdat_params,varargin)
%
% [ids_equilibrium,ids_equilibrium_description,varargout] = tcv_get_ids_equilibrium(shot,ids_equil_empty,varargin);
%
@@ -6,7 +6,7 @@ function [ids_equilibrium,ids_equilibrium_description,varargout] = tcv_get_ids_e
% gdat_params: gdat_data.gdat_params to get all params passed from original call, in particular error_bar and cocos_out options
%
-if exist('gdat_params')
+if exist('gdat_params','var')
[ids_equilibrium, params_equilibrium] = tcv_ids_headpart(shot,ids_equil_empty,'equilibrium','gdat_params',gdat_params,varargin{:});
else
[ids_equilibrium, params_equilibrium] = tcv_ids_headpart(shot,ids_equil_empty,'equilibrium',varargin{:});
@@ -22,14 +22,14 @@ end
%
%% liuqe.m at this stage is missing correction 0.996, so use std source by time of default liuqe to make sure
params_eff_ref = gdat_params; params_eff_ref.doplot=0;
-try;params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def)
+try params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def)
params_eff = params_eff_ref;
params_eff.data_request='b0'; params_eff.source='liuqe'; % to get liuqe time array
bb = gdat(params_equilibrium.shot,params_eff);
params_eff = rmfield(params_eff,'source'); % to get original magnetics data
vacuum_toroidal_field.b0=gdat(params_equilibrium.shot,params_eff);
-ij_ok = [isfinite(vacuum_toroidal_field.b0.data)];
+ij_ok = isfinite(vacuum_toroidal_field.b0.data);
vacuum_toroidal_field.b0.data = interpos(21,vacuum_toroidal_field.b0.t(ij_ok),vacuum_toroidal_field.b0.data(ij_ok),bb.t);
vacuum_toroidal_field.b0.t = bb.t;
vacuum_toroidal_field.b0.dim = {vacuum_toroidal_field.b0.t};
@@ -41,15 +41,16 @@ ids_equilibrium_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;
ids_equilibrium.time = vacuum_toroidal_field.b0.t;
ids_equilibrium_description.time = '.t subfield from: [''b0'',''source'',''liuqe'']';
+ntime = numel(ids_equilibrium.time);
-ids_equilibrium.time_slice(1:numel(ids_equilibrium.time)) = ids_equilibrium.time_slice(1);
+ids_equilibrium.time_slice(1:ntime) = ids_equilibrium.time_slice(1);
% load time array data to copy to time_slices
% global_quantities data into local global_quantities.* structure with correct end names and global_quantities_desc.* with description. Use temp.* and temp_desc.* structures for temporary data
% brute force solution load all eqdsks
-% $$$ for it=1:numel(ids_equilibrium.time)
+% $$$ for it=1:ntime
% $$$ ids_equilibrium.time(it)
% $$$ temp.eqdsks{it}=gdat(params_equilibrium.shot,'eqdsk','time',ids_equilibrium.time(it),'write',0,'machine',gdat_params.machine);
% $$$ end
@@ -109,7 +110,7 @@ global_quantities_desc.w_mhd = params_eff.data_request;
global_quantities_fieldnames = fieldnames(global_quantities);
special_fields = {'magnetic_axis', 'psi_axis', 'q_min'}; % fields needing non-automatic treatments
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
for i=1:numel(global_quantities_fieldnames)
if ~any(strcmp(global_quantities_fieldnames{i},special_fields))
if ~isstruct(ids_equilibrium.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}))
@@ -123,7 +124,7 @@ for it=1:numel(ids_equilibrium.time)
end
% special case
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.r = temp.r_magnetic_axis.data(it);
ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.z = temp.z_magnetic_axis.data(it);
ids_equilibrium.time_slice{it}.global_quantities.psi_axis = temp.psi_axis.data(it) + ...
@@ -181,7 +182,7 @@ temp_desc.z_lcfs = params_eff.data_request;
boundary_fieldnames = fieldnames(boundary);
special_fields = {'lcfs', 'geometric_axis', 'x_point'}; % fields needing non-automatic treatments
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
for i=1:numel(boundary_fieldnames)
if ~any(strcmp(boundary_fieldnames{i},special_fields))
if ~isstruct(ids_equilibrium.time_slice{it}.boundary.(boundary_fieldnames{i}))
@@ -195,7 +196,7 @@ for it=1:numel(ids_equilibrium.time)
end
% special cases
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
ids_equilibrium.time_slice{it}.boundary.outline.r = temp.r_lcfs.data(:,it);
ids_equilibrium.time_slice{it}.boundary.outline.z = temp.z_lcfs.data(:,it);
ids_equilibrium.time_slice{it}.boundary.lcfs.r = ids_equilibrium.time_slice{it}.boundary.outline.r;
@@ -214,6 +215,115 @@ for it=1:numel(ids_equilibrium.time)
end
end
+%% constraints
+% TODO: Add description
+
+% Measured values
+liuqe_time = ids_equilibrium.time; % Not necessarily magnetics time so far
+mag_time = mdsvalue('\magnetics::bpol_003:dim0');
+itime = iround_os(mag_time, liuqe_time);
+mag_time_req = mdscvt(mag_time(itime),'f');
+bpol = mdsvalue('\magnetics::bpol_003[$1,*]',mag_time_req);
+flux = mdsvalue('\magnetics::flux[$1,*]',mag_time_req);
+diam = mdsvalue('\results::dmlcor[$1]',mag_time_req);
+ip = mdsvalue('\magnetics::iplasma:trapeze[$1]',mag_time_req);
+% Coil currents since dim of constraints pf_current is IDS:pf_active/coil
+dim_pol = {'OH_001','OH_002','OH_002','OH_002','OH_002','OH_002','OH_002',...
+ 'E_001','E_002','E_003','E_004','E_005','E_006','E_007','E_008',...
+ 'F_001','F_002','F_003','F_004','F_005','F_006','F_007','F_008',...
+ 'G_001','G_001','G_001','G_001','G_001','G_001'};
+ipol = mdsvalue('\magnetics::ipol[$1,$2]',mag_time_req,dim_pol);
+dim_pol(30:32) = {'TOR_001'};
+ipol(:,30:32) = 0; % TOR_001 is not used in LIUQE
+
+% Reconstructed values
+ipol_liuqe_order = [18,19*ones(1,6),1:16,17*ones(1,6)]; % LIUQE order is E F G OH
+bpol_liuqe = mdsvalue('tcv_eq("b_probe","liuqe.m")');
+flux_liuqe = mdsvalue('tcv_eq("psi_loop","liuqe.m")');
+diam_liuqe = mdsvalue('tcv_eq("tor_flux_dml","liuqe.m")');
+ip_liuqe = mdsvalue('tcv_eq("i_pl","liuqe.m")');
+ipol_liuqe = mdsvalue('tcv_eq("i_pol","liuqe.m")');
+ipol_liuqe = ipol_liuqe(ipol_liuqe_order,:);
+ipol_liuqe(20,:) = 0; % ... TOR
+
+% Weights (using old parameters tree for now)
+bpol_err = mdsvalue('\results::parameters:berr')./mdsvalue('\results::parameters:vvv[0:37]');
+flux_err = mdsvalue('\results::parameters:ferr')./mdsvalue('\results::parameters:www[0:37]')*2*pi;
+diam_err = 0.13e-3./mdsvalue('\results::parameters:idml');
+ip_err = mdsvalue('\results::parameters:plcerr')*1e3;
+ipol_err = mdsvalue('\results::parameters:cerr')./mdsvalue('\results::parameters:uuu[0:18]')*1e3;
+ipol_err = ipol_err(ipol_liuqe_order);
+ipol_err(20) = NaN;
+
+if ntime > 0
+ constraints_orig = ids_equilibrium.time_slice{1}.constraints;
+ % Remove unused arrays
+ ununsed_constraints = {'faraday_angle','mse_polarisation_angle','iron_core_segment',...
+ 'n_e','n_e_line','pressure','q','x_point'};
+ for name = ununsed_constraints, constraints_orig.(name{1})={}; end
+end
+for it = 1:ntime
+ constraints = constraints_orig;
+ % bpol_probe
+ nbpol = size(bpol,2);
+ bpol_probe(1:nbpol) = constraints.bpol_probe(1);
+ for ib = 1:nbpol
+ bpol_probe{ib}.measured = bpol(it,ib);
+ bpol_probe{ib}.source = sprintf('IDS:magnetics/bpol_probe[%02d]/field',ib);
+ bpol_probe{ib}.time_measurement = mag_time(itime(it));
+ bpol_probe{ib}.exact = 0;
+ bpol_probe{ib}.weight = 1/(bpol_err(ib)).^2;
+ bpol_probe{ib}.reconstructed = bpol_liuqe(ib,it);
+ end
+ constraints.bpol_probe = bpol_probe;
+ % flux_loop
+ nflux = size(flux,2);
+ flux_loop(1:nflux) = constraints.flux_loop(1);
+ for il = 1:nflux
+ flux_loop{il}.measured = flux(it,il);
+ flux_loop{il}.source = sprintf('IDS:magnetics/flux_loop[%02d]/flux',il);
+ flux_loop{il}.time_measurement = mag_time(itime(it));
+ flux_loop{il}.exact = 0;
+ flux_loop{il}.weight = 1/(flux_err(il)).^2;
+ flux_loop{il}.reconstructed = flux_liuqe(il,it);
+ end
+ constraints.flux_loop = flux_loop;
+ % ip
+ constraints.ip.measured = ip(it);
+ constraints.ip.source = 'IDS:magnetics/method[1]/ip';
+ constraints.ip.time_measurement = mag_time(itime(it));
+ constraints.ip.exact = 0;
+ constraints.ip.weight = 1/(ip_err).^2;
+ constraints.ip.reconstructed = ip_liuqe(it);
+ % diamagnetic_flux
+ constraints.diamagnetic_flux.measured = diam(it);
+ constraints.diamagnetic_flux.source = 'IDS:magnetics/method[1]/diamagnetic_flux';
+ constraints.diamagnetic_flux.time_measurement = mag_time(itime(it));
+ constraints.diamagnetic_flux.exact = 0;
+ constraints.diamagnetic_flux.weight = 1/(diam_err).^2;
+ constraints.diamagnetic_flux.reconstructed = diam_liuqe(it);
+ % pf_current
+ nipol = size(ipol,2);
+ pf_current(1:nipol) = constraints.pf_current(1);
+ for ic = 1:nipol
+ pf_current{ic}.measured = ipol(it,ic);
+ pf_current{ic}.source = sprintf('IDS:pf_active/coil[%02d]/current',ic);
+ pf_current{ic}.time_measurement = mag_time(itime(it));
+ if strcmp(dim_pol{ic},'TOR_001')
+ pf_current{ic}.source = [pf_current{ic}.source,' replaced with 0'];
+ pf_current{ic}.exact = 1;
+ else
+ pf_current{ic}.exact = 0;
+ pf_current{ic}.weight = 1/(ipol_err(ic)).^2;
+ pf_current{ic}.reconstructed = ipol_liuqe(ic,it);
+ end
+ end
+ constraints.pf_current = pf_current;
+
+ ids_equilibrium.time_slice{it}.constraints = constraints;
+end
+
+
%
%% profiles_1d (cannot use eqdsk since not same radial mesh)
%
@@ -314,7 +424,9 @@ profiles_1d.gm9.data = tmp_gm.x;
profiles_1d_desc.gm9 = 'FS_av(1./R.^1)';
tmp_gm = FS_av(grho_metric_3D.^2./R.^2); % grad rhopol^2 to get <grad psi^2>
-for it=1:numel(ids_equilibrium.time)
+nrho = numel(profiles_1d.rho_tor.x);
+gradpsi_over_R_sq = NaN(nrho,ntime);
+for it=1:ntime
gradpsi_over_R_sq(:,it) = tmp_gm.x(:,it) .* 4 .* profiles_1d.volume.x.^2 .* ...
(ids_equilibrium.time_slice{it}.global_quantities.psi_boundary-ids_equilibrium.time_slice{it}.global_quantities.psi_axis).^2;
end
@@ -335,7 +447,7 @@ profiles_1d.j_parallel.data = j_par./repmat(ids_equilibrium.vacuum_toroidal_fiel
profiles_1d_fieldnames = fieldnames(profiles_1d);
special_fields = {'geometric_axis', 'rho_tor_norm', 'psi'}; % fields needing non-automatic treatments
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
for i=1:numel(profiles_1d_fieldnames)
if ~any(strcmp(profiles_1d_fieldnames{i},special_fields))
if ~isstruct(ids_equilibrium.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}))
@@ -352,9 +464,7 @@ for it=1:numel(ids_equilibrium.time)
end
% special cases
-nrho = numel(profiles_1d.rho_tor.x);
-ntime = numel(temp.psi_axis.data);
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.b_field_tor = ids_equilibrium.time_slice{it}.profiles_1d.f(1) ...
./ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.r;
ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.b_tor = ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.b_field_tor;
@@ -390,7 +500,7 @@ profiles_2d_desc.psi = [params_eff.data_request ' adding psi_bound in a 2nd step
profiles_2d_fieldnames = fieldnames(profiles_2d);
special_fields = {'grid', 'grid_type'}; % fields needing non-automatic treatments
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
for i=1:numel(profiles_2d_fieldnames)
if ~any(strcmp(profiles_2d_fieldnames{i},special_fields))
if ~isstruct(ids_equilibrium.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}))
@@ -407,7 +517,7 @@ for it=1:numel(ids_equilibrium.time)
end
% special cases
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
ids_equilibrium.time_slice{it}.profiles_2d{1}.grid_type.name = profiles_2d.grid_type.name;
ids_equilibrium.time_slice{it}.profiles_2d{1}.grid_type.index = profiles_2d.grid_type.index;
ids_equilibrium.time_slice{it}.profiles_2d{1}.grid_type.description = profiles_2d.grid_type.description;
@@ -424,17 +534,6 @@ for it=1:numel(ids_equilibrium.time_slice)
ids_equilibrium.time_slice{it}.boundary.strike_point = {};
ids_equilibrium.time_slice{it}.boundary_separatrix.x_point = {};
ids_equilibrium.time_slice{it}.boundary_separatrix.strike_point = {};
- ids_equilibrium.time_slice{it}.constraints.bpol_probe = {};
- ids_equilibrium.time_slice{it}.constraints.faraday_angle = {};
- ids_equilibrium.time_slice{it}.constraints.mse_polarisation_angle = {};
- ids_equilibrium.time_slice{it}.constraints.flux_loop = {};
- ids_equilibrium.time_slice{it}.constraints.iron_core_segment = {};
- ids_equilibrium.time_slice{it}.constraints.n_e = {};
- ids_equilibrium.time_slice{it}.constraints.n_e_line = {};
- ids_equilibrium.time_slice{it}.constraints.pf_current = {};
- ids_equilibrium.time_slice{it}.constraints.pressure = {};
- ids_equilibrium.time_slice{it}.constraints.q = {};
- ids_equilibrium.time_slice{it}.constraints.x_point = {};
end
% special test matrix cocos transform
@@ -465,7 +564,7 @@ end
% $$$ end
% cocos automatic transform
-if exist('ids_generic_cocos_nodes_transformation_symbolic') == 2
+if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic'))
[ids_equilibrium,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_equilibrium,'equilibrium',gdat_params.cocos_in, ...
gdat_params.cocos_out,gdat_params.ipsign_out,gdat_params.b0sign_out,gdat_params.ipsign_in,gdat_params.b0sign_in, ...
gdat_params.error_bar,gdat_params.nverbose);
diff --git a/matlab/TCV_IMAS/tcv_get_ids_magnetics.m b/matlab/TCV_IMAS/tcv_get_ids_magnetics.m
index 8f405901..b20f74f1 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_magnetics.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_magnetics.m
@@ -24,7 +24,7 @@ ids_magnetics.b_field_tor_probe = {};
ids_magnetics.rogowski_coil = {};
% cocos automatic transform
-if exist('ids_generic_cocos_nodes_transformation_symbolic','file')
+if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic'))
[ids_magnetics,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_magnetics,'magnetics',gdat_params.cocos_in, ...
gdat_params.cocos_out,gdat_params.ipsign_out,gdat_params.b0sign_out,gdat_params.ipsign_in,gdat_params.b0sign_in, ...
gdat_params.error_bar,gdat_params.nverbose);
--
GitLab