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

add 1st version of ids_summary

parent 6e995a22
No related branches found
No related tags found
1 merge request!58Add ids summary
Pipeline #40383 passed
Stage: test
Hide whitespace changes
Inline Side-by-side
matlab/TCV_IMAS/tcv_get_ids_summary.m
+ 144
− 522
View file @ d9fe50d1
......@@ -39,10 +39,20 @@ ids_summary.code.repository = git_url;
% use equilibrium time as reference, so load first elongation(time)
params_eff = params_eff_ref;
params_eff.data_request='kappa'; params_eff.source='liuqe'; % to get liuqe time array
params_eff.data_request='ip_trapeze';
ip_trapeze = gdat(params_summary.shot,params_eff);
ipsign = mdsvalue('\pcs::data:iohfb');
i_t_end = find(ipsign.*ip_trapeze.data(end:-1:1) > 5e3); % an offset might cause a problem...
i_t_end = length(ip_trapeze.t)-i_t_end(1)+2;
t_end = ip_trapeze.t(i_t_end);
b0sign = mdsvalue('\pcs::data:if36fb');
params_eff = params_eff_ref;
params_eff.data_request='kappa';
kappa = gdat(params_summary.shot,params_eff);
ids_summary.time = kappa.t;
ids_summary_description.time = 'time from gdat(shot,''kappa'')';
% to get liuqe time array for main time, but add points before and after according to ip_trapeze, to have other data when available
dt_caract = 1e-3;
ids_summary.time = unique([[0:dt_caract:kappa.t(1)]'; kappa.t; [kappa.t(end):dt_caract:t_end]']);
ids_summary_description.time = 'time from gdat(shot,''kappa'') and adding from 0 up to t(ip_trapeze(end))';
params_eff = params_eff_ref;
% boundary:
......@@ -81,11 +91,12 @@ boundary_desc.triangularity_upper = params_eff.data_request;
boundary_fieldnames = fieldnames(boundary);
special_fields = {}; % fields needing non-automatic treatments
% assume all boundary have same time base since called with same gdat_params for liuqe
for i=1:numel(boundary_fieldnames)
if ~any(strcmp(boundary_fieldnames{i},special_fields))
if ~isstruct(ids_summary.boundary.(boundary_fieldnames{i}).value)
ids_summary.boundary.(boundary_fieldnames{i}).value = boundary.(boundary_fieldnames{i}).data;
ids_summary.boundary.(boundary_fieldnames{i}).value = interp1(boundary.(boundary_fieldnames{i}).t, ...
boundary.(boundary_fieldnames{i}).data,ids_summary.time,'linear',NaN);
ids_summary.boundary.(boundary_fieldnames{i}).source = ['gdat request: ' boundary_desc.(boundary_fieldnames{i})];
else
special_fields{end+1} = boundary_fieldnames{i};
......@@ -112,549 +123,160 @@ end
[ids_summary_out,ids_summary_out_description] = ids_summary_disruption(shot, ids_summary, gdat_params);
return
% to continue later...
params_eff = params_eff_ref;
params_eff.data_request='b0'; params_eff.source='liuqe'; % to get liuqe time array
bb = gdat(params_summary.shot,params_eff);
params_eff = rmfield(params_eff,'source'); % to get original magnetics data
vacuum_toroidal_field.b0=gdat(params_summary.shot,params_eff);
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};
vacuum_toroidal_field_desc.b0 = ['''b0'',''timing source'',''liuqe=' num2str(gdat_params.liuqe) ''''];
vacuum_toroidal_field_desc.r0 = '.r0 subfield from: [''b0'',''source'',''liuqe'']';
ids_summary.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
ids_summary.vacuum_toroidal_field.b0 = vacuum_toroidal_field.b0.data;
ids_summary_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;
ids_summary.time = vacuum_toroidal_field.b0.t;
ids_summary_description.time = '.t subfield from: [''b0'',''source'',''liuqe'']';
ntime = numel(ids_summary.time);
ids_summary.time_slice(1:ntime) = ids_summary.time_slice(1);
% gas: bottom/top?, deuterium, helium_xx, ?
params_eff.data_request = 'gas_flux';
gas_injection_rates.deuterium = gdat(params_summary.shot,params_eff);
gas_injection_rates_desc.deuterium = params_eff.data_request;
% 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:ntime
% $$$ ids_summary.time(it)
% $$$ temp.eqdsks{it}=gdat(params_summary.shot,'eqdsk','time',ids_summary.time(it),'write',0,'machine',gdat_params.machine);
% $$$ end
% $$$ temp_desc.eqdsks{1} = '''eqdsk'',''time'',ids_summary.time(it)';
params_eff = params_eff_ref;
params_eff.data_request = 'area_edge';
global_quantities.area = gdat(params_summary.shot,params_eff);
global_quantities_desc.area = params_eff.data_request;
params_eff.data_request = 'betan';
global_quantities.beta_normal = gdat(params_summary.shot,params_eff);
global_quantities_desc.beta_normal = params_eff.data_request;
params_eff.data_request = 'betap';
gas_injection_rates_fieldnames = fieldnames(gas_injection_rates);
special_fields = {}; % fields needing non-automatic treatments
% assume all gas_injection_rates have same time base since called with same gdat_params for liuqe
for i=1:numel(gas_injection_rates_fieldnames)
if ~any(strcmp(gas_injection_rates_fieldnames{i},special_fields))
if ~isstruct(ids_summary.gas_injection_rates.(gas_injection_rates_fieldnames{i}).value)
ids_summary.gas_injection_rates.(gas_injection_rates_fieldnames{i}).value = interp1( ...
gas_injection_rates.(gas_injection_rates_fieldnames{i}).t,gas_injection_rates.(gas_injection_rates_fieldnames{i}).data, ...
ids_summary.time,'linear',NaN);
ids_summary.gas_injection_rates.(gas_injection_rates_fieldnames{i}).source = ['gdat request: ' gas_injection_rates_desc.(gas_injection_rates_fieldnames{i})];
else
special_fields{end+1} = gas_injection_rates_fieldnames{i};
end
end
end
%params_eff.doplot=1;
% global_quantities:
params_eff.data_request = 'b0';
global_quantities.b0 = gdat(params_summary.shot,params_eff); % do not add 'source','liuqe' to get all times
global_quantities_desc.b0 = [params_eff.data_request ' ; ' global_quantities.b0.data_fullpath];
params_eff.data_request = 'beta_pol';
global_quantities.beta_pol = gdat(params_summary.shot,params_eff);
global_quantities_desc.beta_pol = params_eff.data_request;
params_eff.data_request = 'beta';
global_quantities_desc.beta_pol = [params_eff.data_request ' ; ' global_quantities.beta_pol.data_fullpath];
params_eff.data_request = 'beta_pol';
global_quantities.beta_pol_mhd = gdat(params_summary.shot,params_eff);
global_quantities_desc.beta_pol_mhd = [params_eff.data_request ' ; ' global_quantities.beta_pol_mhd.data_fullpath];
params_eff.data_request = 'beta_tor';
global_quantities.beta_tor = gdat(params_summary.shot,params_eff);
global_quantities_desc.beta_tor = params_eff.data_request;
global_quantities_desc.beta_tor = [params_eff.data_request ' ; ' global_quantities.beta_tor.data_fullpath];
params_eff.data_request = 'beta_tor';
global_quantities.beta_tor_mhd = gdat(params_summary.shot,params_eff);
global_quantities_desc.beta_tor_mhd = [params_eff.data_request ' ; ' global_quantities.beta_tor_mhd.data_fullpath];
params_eff.data_request = 'betan';
global_quantities.beta_tor_norm = gdat(params_summary.shot,params_eff);
global_quantities_desc.beta_tor_norm = [params_eff.data_request ' ; ' global_quantities.beta_tor_norm.data_fullpath];
params_eff.data_request = 'betan';
global_quantities.beta_tor_norm_mhd = gdat(params_summary.shot,params_eff);
global_quantities_desc.beta_tor_norm_mhd = [params_eff.data_request ' ; ' global_quantities.beta_tor_norm_mhd.data_fullpath];
params_eff.data_request = 'betan';
global_quantities.beta_tor_thermal_norm = gdat(params_summary.shot,params_eff);
global_quantities_desc.beta_tor_thermal_norm = [params_eff.data_request ' ; ' global_quantities.beta_tor_thermal_norm.data_fullpath];
params_eff.data_request = '\tcv_shot::top.results.ibs:ibs';
global_quantities.current_bootstrap = gdat(params_summary.shot,params_eff);
global_quantities_desc.current_bootstrap = [params_eff.data_request ' ; ' global_quantities.current_bootstrap.data_fullpath];
params_eff.data_request = '\tcv_shot::top.results.ibs:icd';
global_quantities.current_non_inductive = gdat(params_summary.shot,params_eff);
global_quantities_desc.current_non_inductive = [params_eff.data_request ' ; ' global_quantities.current_non_inductive.data_fullpath];
params_eff.data_request = 'w_mhd';
global_quantities.energy_diamagnetic = gdat(params_summary.shot,params_eff);
global_quantities_desc.energy_diamagnetic = [params_eff.data_request ' ; ' global_quantities.energy_diamagnetic.data_fullpath];
params_eff.data_request = '\tcv_shot::top.results.conf:we';
global_quantities.energy_electrons_thermal = gdat(params_summary.shot,params_eff);
global_quantities_desc.energy_electrons_thermal = [params_eff.data_request ' ; ' global_quantities.energy_electrons_thermal.data_fullpath];
params_eff.data_request = '\tcv_shot::top.results.conf:wi';
global_quantities.energy_ion_total_thermal = gdat(params_summary.shot,params_eff);
global_quantities_desc.energy_ion_total_thermal = [params_eff.data_request ' ; ' global_quantities.energy_ion_total_thermal.data_fullpath];
params_eff.data_request = 'w_mhd';
global_quantities.energy_mhd = gdat(params_summary.shot,params_eff);
global_quantities_desc.energy_mhd = params_eff.data_request;
global_quantities_desc.energy_mhd = [params_eff.data_request ' ; ' global_quantities.energy_mhd.data_fullpath];
params_eff.data_request = '.conf:we + .conf:wi';
global_quantities.energy_thermal = global_quantities.energy_electrons_thermal;
global_quantities.energy_thermal.data = global_quantities.energy_thermal.data + global_quantities.energy_ion_total_thermal.data;
global_quantities.energy_thermal.data_fullpath = [global_quantities.energy_thermal.data_fullpath ' + wi'];
global_quantities_desc.energy_thermal = [params_eff.data_request ' ; ' global_quantities.energy_thermal.data_fullpath];
params_eff.data_request = 'ngf';
global_quantities.greenwald_fraction = gdat(params_summary.shot,params_eff);
global_quantities_desc.greenwald_fraction = [params_eff.data_request ' ; ' global_quantities.greenwald_fraction.data_fullpath];
params_eff.data_request = 'h98y2';
global_quantities.h_98 = gdat(params_summary.shot,params_eff);
global_quantities_desc.h_98 = [params_eff.data_request ' ; ' global_quantities.h_98.data_fullpath];
params_eff.data_request = 'ip';
global_quantities.ip = gdat(params_summary.shot,params_eff);
global_quantities_desc.ip = params_eff.data_request;
% length_pol = gdat(params_summary.shot,'length_pol','machine',gdat_params.machine); % to be added
global_quantities_desc.ip = [params_eff.data_request ' ; ' global_quantities.ip.data_fullpath];
params_eff.data_request = 'li';
global_quantities.li_3 = gdat(params_summary.shot,params_eff);
global_quantities_desc.li_3 = params_eff.data_request;
params_eff.data_request = 'r_axis';
temp.r_magnetic_axis = gdat(params_summary.shot,params_eff);
temp_desc.r_magnetic_axis = params_eff.data_request;
params_eff.data_request = 'z_axis';
temp.z_magnetic_axis = gdat(params_summary.shot,params_eff);
temp_desc.z_magnetic_axis = params_eff.data_request;
params_eff.data_request = 'psi_axis';
temp.psi_axis = gdat(params_summary.shot,params_eff); % needs to add psi_edge sincepsi_axis liuqe assuming 0 dege value
temp_desc.psi_axis = params_eff.data_request;
params_eff.data_request = 'psi_edge';
global_quantities.psi_boundary = gdat(params_summary.shot,params_eff);
global_quantities_desc.psi_boundary = params_eff.data_request;
global_quantities.li = gdat(params_summary.shot,params_eff);
global_quantities_desc.li = [params_eff.data_request ' ; ' global_quantities.li.data_fullpath];
params_eff.data_request = '\tcv_shot::top.results.conf:ptot_ohm';
global_quantities.power_ohm = gdat(params_summary.shot,params_eff);
global_quantities_desc.power_ohm = [params_eff.data_request ' ; ' global_quantities.power_ohm.data_fullpath];
params_eff.data_request = '\results::bolo:prad:total';
global_quantities.power_radiated = gdat(params_summary.shot,params_eff);
if strcmp(lower(global_quantities.power_radiated.units),'kw')
global_quantities.power_radiated.data = global_quantities.power_radiated.data*1e3; %in kW
end
global_quantities_desc.power_radiated = [params_eff.data_request ' ; ' global_quantities.power_radiated.data_fullpath];
params_eff.data_request = 'q95';
global_quantities.q_95 = gdat(params_summary.shot,params_eff);
global_quantities_desc.q_95 = params_eff.data_request;
params_eff.data_request = 'q0';
global_quantities.q_axis = gdat(params_summary.shot,params_eff); % will be checked with q_rho?
global_quantities_desc.q_axis = params_eff.data_request;
params_eff.data_request = 'q_rho';
temp.q_rho = gdat(params_summary.shot,params_eff);
temp_desc.q_rho = params_eff.data_request;
% surface = gdat(params_summary.shot,'surface','machine',gdat_params.machine); % to be added
global_quantities_desc.q_95 = [params_eff.data_request ' ; ' global_quantities.q_95.data_fullpath];
params_eff.data_request = 'b0';
global_quantities.b0 = gdat(params_summary.shot,params_eff);
global_quantities_desc.b0 = [params_eff.data_request ' ; ' global_quantities.b0.data_fullpath];
params_eff.data_request = 'r0';
global_quantities.r0 = global_quantities.b0;
global_quantities.r0.data = global_quantities.r0.r0;
global_quantities.r0.dim = [];global_quantities.r0.t=[]; global_quantities.r0.dimunits = [];
global_quantities.r0.units = 'm';
global_quantities.r0.data_fullpath = ['gdat(b0).r0'];
global_quantities_desc.r0 = [params_eff.data_request ' ; ' global_quantities.r0.data_fullpath];
% $$$ params_eff.data_request = 'tau_tot';
% $$$ global_quantities.tau_energy = gdat(params_summary.shot,params_eff);
% $$$ global_quantities_desc.tau_energy = [params_eff.data_request ' ; ' global_quantities.tau_energy.data_fullpath];
params_eff.data_request = 'vloop';
global_quantities.v_loop = gdat(params_summary.shot,params_eff);
tensvloop=-1e3;
global_quantities.v_loop.data = interpos(global_quantities.v_loop.t,-global_quantities.v_loop.data,tensvloop);
global_quantities_desc.v_loop = [params_eff.data_request ' ; -1* interpos of ' global_quantities.v_loop.data_fullpath ...
' with tens=' num2str(tensvloop)];
params_eff.data_request = 'volume';
global_quantities.volume = gdat(params_summary.shot,params_eff);
global_quantities_desc.volume = params_eff.data_request;
params_eff.data_request = 'w_mhd';
global_quantities.w_mhd = gdat(params_summary.shot,params_eff);
global_quantities_desc.w_mhd = params_eff.data_request;
global_quantities_desc.volume = [params_eff.data_request ' ; ' global_quantities.volume.data_fullpath];
special_fields = {'r0'}; % fields needing non-automatic treatments
global_quantities_fieldnames = fieldnames(global_quantities);
special_fields = {'magnetic_axis', 'psi_axis', 'q_min'}; % fields needing non-automatic treatments
for it=1:ntime
for i=1:numel(global_quantities_fieldnames)
if ~any(strcmp(global_quantities_fieldnames{i},special_fields))
if ~isstruct(ids_summary.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}))
ids_summary.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}) = ...
global_quantities.(global_quantities_fieldnames{i}).data(it);
else
special_fields{end+1} = global_quantities_fieldnames{i};
end
for i=1:numel(global_quantities_fieldnames)
if ~any(strcmp(global_quantities_fieldnames{i},special_fields))
if ~isstruct(ids_summary.global_quantities.(global_quantities_fieldnames{i}).value)
ids_summary.global_quantities.(global_quantities_fieldnames{i}).value = interp1( ...
global_quantities.(global_quantities_fieldnames{i}).t,global_quantities.(global_quantities_fieldnames{i}).data, ...
ids_summary.time,'linear',NaN);
ids_summary.global_quantities.(global_quantities_fieldnames{i}).source = ['gdat request: ' global_quantities_desc.(global_quantities_fieldnames{i})];
else
special_fields{end+1} = global_quantities_fieldnames{i};
end
end
end
% special case
for it=1:ntime
ids_summary.time_slice{it}.global_quantities.magnetic_axis.r = temp.r_magnetic_axis.data(it);
ids_summary.time_slice{it}.global_quantities.magnetic_axis.z = temp.z_magnetic_axis.data(it);
ids_summary.time_slice{it}.global_quantities.psi_axis = temp.psi_axis.data(it) + ...
ids_summary.time_slice{it}.global_quantities.psi_boundary;
[zz,izz]=min(temp.q_rho.data(:,it));
ids_summary.time_slice{it}.global_quantities.q_min.value = zz;
ids_summary.time_slice{it}.global_quantities.q_min.rho_tor_norm = temp.q_rho.grids_1d.rhotornorm(izz);
end
% for boundary in addition to lcfs
% active_limiter_point = gdat(params_summary.shot,'active_limiter_point','machine',gdat_params.machine);
params_eff.data_request = 'kappa';
boundary.elongation = gdat(params_summary.shot,params_eff);
boundary_desc.elongation = params_eff.data_request;
% elongation_lower = gdat(params_summary.shot,'elongation_lower','machine',gdat_params.machine);
% elongation_upper = gdat(params_summary.shot,'elongation_upper','machine',gdat_params.machine);
params_eff.data_request = 'a_minor';
boundary.minor_radius = gdat(params_summary.shot,params_eff);
boundary_desc.minor_radius = params_eff.data_request;
% squareness_lower_inner = gdat(params_summary.shot,'squareness_lower_inner','machine',gdat_params.machine);
% squareness_lower_outer = gdat(params_summary.shot,'squareness_lower_outer','machine',gdat_params.machine);
% squareness_upper_inner = gdat(params_summary.shot,'squareness_upper_inner','machine',gdat_params.machine);
% squareness_upper_outer = gdat(params_summary.shot,'squareness_upper_outer','machine',gdat_params.machine);
% strike_point = gdat(params_summary.shot,'strike_point','machine',gdat_params.machine);
params_eff.data_request = 'delta';
boundary.triangularity = gdat(params_summary.shot,params_eff);
boundary_desc.triangularity = params_eff.data_request;
params_eff.data_request = 'delta_bottom';
boundary.triangularity_lower = gdat(params_summary.shot,params_eff);
boundary_desc.triangularity_lower = params_eff.data_request;
params_eff.data_request = 'delta_top';
boundary.triangularity_upper = gdat(params_summary.shot,params_eff);
boundary_desc.triangularity_upper = params_eff.data_request;
params_eff.data_request = 'tcv_eq(''''n_xpts'''',''''liuqe.m'''')';
temp.n_x_point = gdat(params_summary.shot,params_eff);
temp_desc.n_x_point = params_eff.data_request;
params_eff.data_request = 'tcv_eq(''''r_xpts'''',''''liuqe.m'''')';
temp.r_x_point = gdat(params_summary.shot,params_eff);
temp_desc.r_x_point = params_eff.data_request;
params_eff.data_request = 'tcv_eq(''''z_xpts'''',''''liuqe.m'''')';
temp.z_x_point = gdat(params_summary.shot,params_eff);
temp_desc.z_x_point = params_eff.data_request;
params_eff.data_request = 'rgeom';
temp.rgeom = gdat(params_summary.shot,params_eff);
temp_desc.rgeom = params_eff.data_request;
params_eff.data_request = 'zgeom';
temp.zgeom = gdat(params_summary.shot,params_eff);
temp_desc.zgeom = params_eff.data_request;
params_eff.data_request = 'r_contour_edge';
temp.r_lcfs = gdat(params_summary.shot,params_eff);
temp_desc.r_lcfs = params_eff.data_request;
params_eff.data_request = 'z_contour_edge';
temp.z_lcfs = gdat(params_summary.shot,params_eff);
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:ntime
for i=1:numel(boundary_fieldnames)
if ~any(strcmp(boundary_fieldnames{i},special_fields))
if ~isstruct(ids_summary.time_slice{it}.boundary.(boundary_fieldnames{i}))
ids_summary.time_slice{it}.boundary.(boundary_fieldnames{i}) = ...
boundary.(boundary_fieldnames{i}).data(it);
else
special_fields{end+1} = boundary_fieldnames{i};
end
end
end
end
% special cases
for it=1:ntime
ids_summary.time_slice{it}.boundary.outline.r = temp.r_lcfs.data(:,it);
ids_summary.time_slice{it}.boundary.outline.z = temp.z_lcfs.data(:,it);
ids_summary.time_slice{it}.boundary.lcfs.r = ids_summary.time_slice{it}.boundary.outline.r;
ids_summary.time_slice{it}.boundary.lcfs.z = ids_summary.time_slice{it}.boundary.outline.z;
ids_summary.time_slice{it}.boundary.geometric_axis.r = temp.rgeom.data(it);
ids_summary.time_slice{it}.boundary.geometric_axis.z = temp.zgeom.data(it);
if temp.n_x_point.data(it) > 0
% not that asking for time_out may lead to interpolated nb of X-points non integer, should included piece-wise constant interpolation for equil quantities...
ids_summary.time_slice{it}.boundary.x_point(1:fix(temp.n_x_point.data(it))) = ids_summary.time_slice{it}.boundary.x_point(1);
for i=1:fix(temp.n_x_point.data(it))
ids_summary.time_slice{it}.boundary.x_point{i}.r = temp.r_x_point.data(i,it);
ids_summary.time_slice{it}.boundary.x_point{i}.z = temp.z_x_point.data(i,it);
end
else
ids_summary.time_slice{it}.boundary.x_point = {};
end
end
%% constraints
% TODO: Add description
% Measured values
liuqe_time = ids_summary.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('tcv_idealloop("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);
ipol(:,27:29) = -ipol(:,27:29); % Bottom G-coils
dim_pol(30:32) = {'TOR_001'};
ipol(:,30:32) = 0; % TOR_001 is not used in LIUQE
ids_summary.global_quantities.global_quantities.r0.value = global_quantities.r0.data;
ids_summary.global_quantities.global_quantities.r0.source = ['gdat request: ' global_quantities_desc.r0];
% 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(27:29,:) = -ipol_liuqe(27:29,:); % Bottom G-coils
ipol_liuqe(30:32,:) = 0; % ... TOR
% cst zeff
params_eff.data_request = '\tcv_shot::top.results.ibs:z_eff';
volume_average.zeff = gdat(params_summary.shot,params_eff);
volume_average_desc.zeff = [params_eff.data_request ' ; ' volume_average.zeff.data_fullpath ' (not vol avrg but to match Iohm)'];
% 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(30:32) = NaN;
if ntime > 0
constraints_orig = ids_summary.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;
special_fields = {}; % fields needing non-automatic treatments
volume_average_fieldnames = fieldnames(volume_average);
for i=1:numel(volume_average_fieldnames)
if ~any(strcmp(volume_average_fieldnames{i},special_fields))
if ~isstruct(ids_summary.volume_average.(volume_average_fieldnames{i}).value)
ids_summary.volume_average.(volume_average_fieldnames{i}).value = interp1( ...
volume_average.(volume_average_fieldnames{i}).t,volume_average.(volume_average_fieldnames{i}).data, ...
ids_summary.time,'linear',NaN);
ids_summary.volume_average.(volume_average_fieldnames{i}).source = ['gdat request: ' volume_average_desc.(volume_average_fieldnames{i})];
else
pf_current{ic}.exact = 0;
pf_current{ic}.weight = 1/(ipol_err(ic)).^2;
pf_current{ic}.reconstructed = ipol_liuqe(ic,it);
special_fields{end+1} = volume_average_fieldnames{i};
end
end
constraints.pf_current = pf_current;
ids_summary.time_slice{it}.constraints = constraints;
end
%
%% profiles_1d (cannot use eqdsk since not same radial mesh)
%
% area = gdat(params_summary.shot,'area','machine',gdat_params.machine);
% b_average = gdat(params_summary.shot,'b_average','machine',gdat_params.machine);
% beta_pol = gdat(params_summary.shot,'beta_pol','machine',gdat_params.machine);
% b_field_average = gdat(params_summary.shot,'b_field_average','machine',gdat_params.machine);
% b_field_max = gdat(params_summary.shot,'b_field_max','machine',gdat_params.machine);
% b_field_min = gdat(params_summary.shot,'b_field_min','machine',gdat_params.machine);
% b_max = gdat(params_summary.shot,'b_max','machine',gdat_params.machine);
% b_min = gdat(params_summary.shot,'b_min','machine',gdat_params.machine);
% darea_dpsi = gdat(params_summary.shot,'darea_dpsi','machine',gdat_params.machine);
% darea_drho_tor = gdat(params_summary.shot,'darea_drho_tor','machine',gdat_params.machine);
params_eff.data_request = 'pprime';
profiles_1d.dpressure_dpsi = gdat(params_summary.shot,params_eff);
profiles_1d_desc.dpressure_dpsi = params_eff.data_request;
% dpsi_drho_tor = gdat(params_summary.shot,'dpsi_drho_tor','machine',gdat_params.machine);
% dvolume_dpsi = gdat(params_summary.shot,'dvolume_dpsi','machine',gdat_params.machine);
% dvolume_drho_tor = gdat(params_summary.shot,'dvolume_drho_tor','machine',gdat_params.machine);
% elongation = gdat(params_summary.shot,'elongation','machine',gdat_params.machine);
params_eff.data_request = 'ttprime';
profiles_1d.f_df_dpsi = gdat(params_summary.shot,params_eff);
profiles_1d_desc.f_df_dpsi = [params_eff.data_request '* 0.996^2'];
params_eff.data_request = 'rbphi_rho';
profiles_1d.f = gdat(params_summary.shot,params_eff);
profiles_1d_desc.f = [params_eff.data_request '* 0.996'];
profiles_1d.f.data = 0.996 * profiles_1d.f.data;
profiles_1d.f_df_dpsi.data = 0.996.^2 .* profiles_1d.f_df_dpsi.data;
% geometric_axis = gdat(params_summary.shot,'geometric_axis','machine',gdat_params.machine);
% gm1 = gdat(params_summary.shot,'gm1','machine',gdat_params.machine);
% gm2 = gdat(params_summary.shot,'gm2','machine',gdat_params.machine);
% gm3 = gdat(params_summary.shot,'gm3','machine',gdat_params.machine);
% gm4 = gdat(params_summary.shot,'gm4','machine',gdat_params.machine);
% gm5 = gdat(params_summary.shot,'gm5','machine',gdat_params.machine);
% gm6 = gdat(params_summary.shot,'gm6','machine',gdat_params.machine);
% gm7 = gdat(params_summary.shot,'gm7','machine',gdat_params.machine);
% gm8 = gdat(params_summary.shot,'gm8','machine',gdat_params.machine);
% gm9 = gdat(params_summary.shot,'gm9','machine',gdat_params.machine);
% j_parallel = gdat(params_summary.shot,'j_parallel','machine',gdat_params.machine);
% j_tor = gdat(params_summary.shot,'j_tor','machine',gdat_params.machine);
% magnetic_shear = gdat(params_summary.shot,'magnetic_shear','machine',gdat_params.machine);
% mass_density = gdat(params_summary.shot,'mass_density','machine',gdat_params.machine);
params_eff.data_request = 'phi_tor';
profiles_1d.phi = gdat(params_summary.shot,params_eff);
profiles_1d.phi.data = 0.996 * profiles_1d.phi.data;
profiles_1d_desc.phi = [params_eff.data_request '* 0.996'];
params_eff.data_request = 'pressure';
profiles_1d.pressure = gdat(params_summary.shot,params_eff);
profiles_1d_desc.pressure = params_eff.data_request;
% psi = gdat(params_summary.shot,'psi_rho','machine',gdat_params.machine); % (could take from .x of any like rhotor and psi_axis, psi_edge from global_quantities)
params_eff.data_request = 'q_rho';
profiles_1d.q = gdat(params_summary.shot,params_eff);
profiles_1d_desc.q = params_eff.data_request;
params_eff.data_request = 'rhotor';
profiles_1d.rho_tor = gdat(params_summary.shot,params_eff);
profiles_1d_desc.rho_tor = params_eff.data_request;
%rho_tor_norm = gdat(params_summary.shot,'rhotor_norm','machine',gdat_params.machine); % from rho_tor
params_eff.data_request = 'rhovol';
profiles_1d.rho_volume_norm = gdat(params_summary.shot,params_eff);
profiles_1d_desc.rho_volume_norm = params_eff.data_request;
% r_inboard = gdat(params_summary.shot,'r_inboard','machine',gdat_params.machine);
% r_outboard = gdat(params_summary.shot,'r_outboard','machine',gdat_params.machine);
% squareness_lower_inner = gdat(params_summary.shot,'squareness_lower_inner','machine',gdat_params.machine);
% squareness_lower_outer = gdat(params_summary.shot,'squareness_lower_outer','machine',gdat_params.machine);
% squareness_upper_inner = gdat(params_summary.shot,'squareness_upper_inner','machine',gdat_params.machine);
% squareness_upper_outer = gdat(params_summary.shot,'squareness_upper_outer','machine',gdat_params.machine);
% surface = gdat(params_summary.shot,'surface','machine',gdat_params.machine);
% trapped_fraction = gdat(params_summary.shot,'trapped_fraction','machine',gdat_params.machine);
% triangularity_lower = gdat(params_summary.shot,'triangularity_lower','machine',gdat_params.machine);
% triangularity_upper = gdat(params_summary.shot,'triangularity_upper','machine',gdat_params.machine);
params_eff.data_request = 'volume_rho';
profiles_1d.volume = gdat(params_summary.shot,params_eff);
profiles_1d_desc.volume = params_eff.data_request;
liuqe_opt = gdat_params.liuqe; % default at this stage but could take from gdat params like error bar
switch liuqe_opt
case {-1}, psitbx_str='FBTE';
case {1,21}, psitbx_str='LIUQE.M';
case {11}, psitbx_str='LIUQE';
case {2, 3, 22, 23}, psitbx_str=['LIUQE.M' num2str(mod(liuqe_opt,10))];
case {12,13}, psitbx_str=['LIUQE' num2str(mod(liuqe_opt,10))];
otherwise, error(['Unknown LIUQE version, liuqe = ' num2str(liuqe_opt)]);
end
fsd = psitbxtcv2(shot,profiles_1d.volume.t,'FS',psitbx_str,false); % will get automatically the correct time interval
grho_metric_3D = metric(fsd,-1);
% Introduced new anonymous function to compute FS average ...
metric_FS = metric(grho_metric_3D.grid,[2,3]);
denom=sum(metric_FS./grho_metric_3D,[2,3]);
FS_av = @(x) sum(x.*metric_FS./grho_metric_3D,[2,3])./denom;
R=metric(fsd,3);
Rm2av=FS_av(1./R.^2);
profiles_1d.gm1.data = Rm2av.x;
profiles_1d_desc.gm1 = ['psitbxtcv2 with ' psitbx_str ' then Rm2av=FS_av(1./R.^2)'];
%tmp_gm = FS_av(grho_metric_3D.^2./R.^2); % this gives (grad rhopol/R)^2 not gm2 which is grad rhotor^2
%profiles_1d.gm2.data = tmp_gm.x;
tmp_gm = FS_av(1./R.^1);
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>
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_summary.time_slice{it}.global_quantities.psi_boundary-ids_summary.time_slice{it}.global_quantities.psi_axis).^2;
end
mu0 = 4.e-7 * pi;
% Eq. (30) cocos paper cocos=17
% j_tor=<jphi/R>/<1/R>=-sigma_Bp (2pi)^e_Bp dp/dpsi / <1/R> - sigma_Bp (2pi)^e_Bp F dF/dpsi / mu0 <1/R^2> / <1/R>
% simaBp=-1 and eBp=1 for cocos=17 from TCV LIUQE
profiles_1d.j_tor.data = - (-1.) .* 2.*pi .* profiles_1d.dpressure_dpsi.data ./ profiles_1d.gm9.data ...
- (-1.) .* 2.*pi .* profiles_1d.gm1.data ./ profiles_1d.gm9.data .* profiles_1d.f_df_dpsi.data ./ mu0;
%
% <j.B> = - sigma_Bp (2pi)^e_Bp dp/dpsi F - sigma_Bp F dF/dpsi / mu0 [ (2pi)^e_Bp F <1/R^2> + 1/(2pi)^e_Bp * <|grad psi|^2/R^2> / F ]
% simaBp=-1 and eBp=1 for cocos=17 from TCV LIUQE
%
j_par = - (-1.) .* 2*pi .* profiles_1d.dpressure_dpsi.data .* profiles_1d.f.data ...
- (-1.) .* profiles_1d.f_df_dpsi.data ./ mu0 .* ...
( (2.*pi) .* profiles_1d.f.data .* profiles_1d.gm1.data + 1./(2.*pi) .* gradpsi_over_R_sq ./ profiles_1d.f.data);
profiles_1d.j_parallel.data = j_par./repmat(ids_summary.vacuum_toroidal_field.b0',size(profiles_1d.f.data,1),1);
profiles_1d_fieldnames = fieldnames(profiles_1d);
special_fields = {'geometric_axis', 'rho_tor_norm', 'psi'}; % fields needing non-automatic treatments
for it=1:ntime
for i=1:numel(profiles_1d_fieldnames)
if ~any(strcmp(profiles_1d_fieldnames{i},special_fields))
if ~isstruct(ids_summary.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}))
if ~ischar(profiles_1d.(profiles_1d_fieldnames{i}).data) && ~isempty(profiles_1d.(profiles_1d_fieldnames{i}).data) ...
&& size(profiles_1d.(profiles_1d_fieldnames{i}).data,2)>=it
ids_summary.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}) = ...
profiles_1d.(profiles_1d_fieldnames{i}).data(:,it);
end
else
special_fields{end+1} = profiles_1d_fieldnames{i};
end
end
end
end
% special cases
for it=1:ntime
ids_summary.time_slice{it}.global_quantities.magnetic_axis.b_field_tor = ids_summary.time_slice{it}.profiles_1d.f(1) ...
./ids_summary.time_slice{it}.global_quantities.magnetic_axis.r;
ids_summary.time_slice{it}.global_quantities.magnetic_axis.b_tor = ids_summary.time_slice{it}.global_quantities.magnetic_axis.b_field_tor;
ids_summary.time_slice{it}.profiles_1d.rho_tor_norm = ids_summary.time_slice{it}.profiles_1d.rho_tor./ ...
ids_summary.time_slice{it}.profiles_1d.rho_tor(end);
ids_summary.time_slice{it}.profiles_1d.psi = ids_summary.time_slice{it}.global_quantities.psi_axis + ...
profiles_1d.rho_tor.x.^2 .* ...
(ids_summary.time_slice{it}.global_quantities.psi_boundary- ids_summary.time_slice{it}.global_quantities.psi_axis);
end
%
%% profiles_2d{1} ala eqdsk, only this one thus grid_type=1
%
% b_field_r = gdat(params_summary.shot,'b_field_r','machine',gdat_params.machine);
% b_field_tor = gdat(params_summary.shot,'b_field_tor','machine',gdat_params.machine);
% b_field_z = gdat(params_summary.shot,'b_field_z','machine',gdat_params.machine);
% b_r = gdat(params_summary.shot,'b_r','machine',gdat_params.machine);
% b_tor = gdat(params_summary.shot,'b_tor','machine',gdat_params.machine);
% b_z = gdat(params_summary.shot,'b_z','machine',gdat_params.machine);
% grid = gdat(params_summary.shot,'grid','machine',gdat_params.machine); % special
profiles_2d.grid_type.name = 'rectangular';
profiles_2d.grid_type.index = 1;
profiles_2d.grid_type.description = 'Cylindrical R,Z ala eqdsk';
% j_parallel = gdat(params_summary.shot,'j_parallel','machine',gdat_params.machine);
% j_tor = gdat(params_summary.shot,'j_tor','machine',gdat_params.machine);
% phi = gdat(params_summary.shot,'phi','machine',gdat_params.machine);
params_eff.data_request = 'psi';
profiles_2d.psi = gdat(params_summary.shot,params_eff); % add psi_bound in a second step in special cases
profiles_2d_desc.psi = [params_eff.data_request ' adding psi_bound in a 2nd step'];
% r = gdat(params_summary.shot,'r','machine',gdat_params.machine); % not to be filled since in grid.dim1
% theta = gdat(params_summary.shot,'theta','machine',gdat_params.machine);
% z = gdat(params_summary.shot,'z','machine',gdat_params.machine); % not to be filled since in grid.dim2
profiles_2d_fieldnames = fieldnames(profiles_2d);
special_fields = {'grid', 'grid_type'}; % fields needing non-automatic treatments
for it=1:ntime
for i=1:numel(profiles_2d_fieldnames)
if ~any(strcmp(profiles_2d_fieldnames{i},special_fields))
if ~isstruct(ids_summary.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}))
if ~ischar(profiles_2d.(profiles_2d_fieldnames{i}).data) && ~isempty(profiles_2d.(profiles_2d_fieldnames{i}).data) ...
&& size(profiles_2d.(profiles_2d_fieldnames{i}).data,3)>=it
ids_summary.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}) = ...
profiles_2d.(profiles_2d_fieldnames{i}).data(:,:,it);
end
else
special_fields{end+1} = profiles_2d_fieldnames{i};
end
end
end
end
% special cases
for it=1:ntime
ids_summary.time_slice{it}.profiles_2d{1}.grid_type.name = profiles_2d.grid_type.name;
ids_summary.time_slice{it}.profiles_2d{1}.grid_type.index = profiles_2d.grid_type.index;
ids_summary.time_slice{it}.profiles_2d{1}.grid_type.description = profiles_2d.grid_type.description;
ids_summary.time_slice{it}.profiles_2d{1}.grid.dim1 = profiles_2d.psi.dim{1};
ids_summary.time_slice{it}.profiles_2d{1}.grid.dim2 = profiles_2d.psi.dim{2};
ids_summary.time_slice{it}.profiles_2d{1}.psi(:,:) = ids_summary.time_slice{it}.profiles_2d{1}.psi(:,:) + ...
global_quantities.psi_boundary.data(it);
end
% make arrays not filled in empty:
ids_summary.grids_ggd = {};
for it=1:numel(ids_summary.time_slice)
ids_summary.time_slice{it}.ggd = {};
ids_summary.time_slice{it}.boundary.strike_point = {};
ids_summary.time_slice{it}.boundary_separatrix.x_point = {};
ids_summary.time_slice{it}.boundary_separatrix.strike_point = {};
end
% special test matrix cocos transform
% $$$ ldim1=129;
% $$$ ldim2=257;
% $$$ it=1;
% $$$ ids_summary.time_slice{it}.coordinate_system.grid_type.index = 13;
% $$$ ids_summary.time_slice{it}.coordinate_system.grid.dim1 = linspace(0,1,ldim1)';
% $$$ ids_summary.time_slice{it}.coordinate_system.grid.dim2 = linspace(0,2*pi,ldim2);
% $$$ ids_summary.time_slice{it}.coordinate_system.tensor_contravariant = 2.*ones(ldim1,ldim2,3,3);
% $$$ ids_summary.time_slice{it}.coordinate_system.tensor_covariant = 0.5*ones(ldim1,ldim2,3,3);
% $$$ ids_summary.time_slice{it}.coordinate_system.g13_contravariant = 13.*ones(ldim1,ldim2,3,3);
% $$$ ids_summary.time_slice{it}.coordinate_system.g13_contravariant_error_upper = 14.*ones(ldim1,ldim2,3,3);
% $$$ ids_summary.time_slice{it}.coordinate_system.g13_contravariant_error_lower = 12.*ones(ldim1,ldim2,3,3);
% $$$ for it=1:2100
% $$$ ids_summary.time_slice{it}.coordinate_system.g11_contravariant = 11.*ones(ldim1,ldim2,3,3);
% $$$ ids_summary.time_slice{it}.coordinate_system.tensor_covariant = 0.5*ones(ldim1,ldim2,3,3);
% $$$ ids_summary.time_slice{it}.profiles_2d{2}.grid_type.name = profiles_2d.grid_type.name;
% $$$ ids_summary.time_slice{it}.profiles_2d{2}.grid_type.index = 11;
% $$$ ids_summary.time_slice{it}.profiles_2d{2}.grid_type.description = profiles_2d.grid_type.description;
% $$$ ids_summary.time_slice{it}.profiles_2d{2}.grid.dim1 = linspace(0,1,ldim1)';
% $$$ ids_summary.time_slice{it}.profiles_2d{2}.grid.dim1_error_upper = 1.2.*linspace(0,1,ldim1)';
% $$$ ids_summary.time_slice{it}.profiles_2d{2}.grid.dim1_error_lower = 0.8.*linspace(0,1,ldim1)';
% $$$ ids_summary.time_slice{it}.profiles_2d{2}.grid.dim2 = linspace(0,2*pi,ldim2);
% $$$ ids_summary.time_slice{it}.profiles_2d{2}.psi(:,:) = 11.*ones(ldim1,ldim2);
% $$$ ids_summary.time_slice{it}.profiles_2d{2}.psi_error_upper(:,:) = 12.*ones(ldim1,ldim2);
% $$$ ids_summary.time_slice{it}.profiles_2d{2}.psi_error_lower(:,:) = 10.*ones(ldim1,ldim2);
% $$$ end
% cocos automatic transform
if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic'))
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment