From 87bccfb3e5b47a8730fd3c5f77345332ceb331f4 Mon Sep 17 00:00:00 2001
From: Antonia Frank <antonia.frank@epfl.ch>
Date: Thu, 8 Aug 2024 17:36:27 +0200
Subject: [PATCH] Finalize the EC section with correct conversion factor &
integration, reduce gdat_powers grid to times when waves are injected to the
plasma
---
matlab/TCV_IMAS/tcv_get_ids_core_sources.m | 203 ++++++++++++---------
1 file changed, 114 insertions(+), 89 deletions(-)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index 9f81305d..6a25f91d 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -52,7 +52,8 @@ powers_gdat = gdat(params_core_sources.shot,params_eff);
% ohm
ohm_t_grid = powers_gdat.ohm.t; ohm_n_t = numel(ohm_t_grid);
-id_ohm.description = 'Source from ohmic heating';
+main_desc = 'Source from ohmic heating'; production = 'IBS nodes';
+id_ohm.description = sprintf('%s from %s',main_desc,production);
id_ohm.index = 7; id_ohm.name = 'ohmic';
ids_core_sources.source{last_index+1} = source_template;
ids_core_sources.source{last_index+1}.identifier = id_ohm;
@@ -106,7 +107,8 @@ if any(bs_t_grid-ohm_t_grid')>0
warning('Bootstrap and ohmic time grids are not the same! Interpolation needed.')
end
-id_bs.description = 'Bootstrap current';
+main_desc = 'Bootstrap current'; production = 'IBS nodes';
+id_bs.description = sprintf('%s from %s',main_desc,production);
id_bs.index = 13; id_bs.name = 'bootstrap_current';
ids_core_sources.source{last_index+1} = source_template;
ids_core_sources.source{last_index+1}.identifier = id_bs;
@@ -136,86 +138,90 @@ ids_core_sources.source{last_index+1}.global_quantities.current_parallel = curre
last_index = last_index+1; % add if statement to only increment if bs source has been added
%% ec
-ec_gdat = gdat(shot,'ec_data','ec_inputs',1);
+params_eff = params_eff_ref;
+params_eff.data_request='ec_data';
+params_eff.ec_inputs = 1; % load EC input information
+ec_gdat = gdat(params_core_sources.shot,params_eff);
if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
ec_data = ec_gdat.ec.ec_data; ec_inputs = ec_gdat.ec.ec_inputs;
% get tgrid from gdat ec_data
- ec_data_tgrid = ec_data.p_dens.t; nt_ec_data = numel(ec_data_tgrid);
-
- % get tgrid from gdat powers
- ec_powers_tgrid = powers_gdat.ec.t; nt_ec_powers = numel(ec_powers_tgrid);
+ ec_tgrid = ec_data.p_dens.t; nt_ec_data = numel(ec_tgrid);
% retrieve active launcher information from ec_inputs
nb_launchers = numel(ec_inputs.launchers_active.data);
active_launchers = find(ec_inputs.launchers_active.data==1)';
+ % find times of injected ec powers, for setting up source and
+ % interpolating the power & current densities on p_inj time grid
+ ec_powers_tgrid = powers_gdat.ec.t; %nt_ec_powers = numel(ec_powers_tgrid);
+ itime_ec = powers_gdat.ec.data(:,end)>0;
+ interp_tgrid = ec_powers_tgrid(itime_ec); nt_interp = numel(interp_tgrid);
+ interp_p_ec_injected = powers_gdat.ec.data(itime_ec,end);
+
% Setup structures for active launchers from template
main_desc = 'Source from electron cyclotron heating and current drive';
- id_ec.index = 3; id_ec.name = 'ec';
+ production = 'TORAY';
+ id_ec.index = 3; id_ec.name = 'ec';
for i_lau = active_launchers
- id_ec.description = sprintf('L%i/G%i, %s',i_lau,ec_inputs.gyro2launcher.data(i_lau),main_desc);
+ id_ec.description = sprintf('L%i/G%i, %s from %s double width CD profiles',i_lau,ec_inputs.gyro2launcher.data(i_lau),main_desc,production);
ids_core_sources.source{last_index+i_lau} = source_template;
ids_core_sources.source{last_index+i_lau}.identifier = id_ec;
- ids_core_sources.source{last_index+i_lau}.profiles_1d(1:nt_ec_powers) = {profiles_template};
+ ids_core_sources.source{last_index+i_lau}.profiles_1d(1:nt_interp) = {profiles_template};
ids_core_sources.source{last_index+i_lau}.global_quantities = globals_template;
end
- % get data for globals from gdat powers and fill in ids_structure
- ec_total_pow = transpose(powers_gdat.ec.data(:,nb_launchers+1)); %use power from powers_gdat(injected ec power) instead of ec_data power
- ec_total_pow(isnan(ec_total_pow)) = 0;
- ec_total_cur = ec_data.cd_tot.data(nb_launchers+1,:);
- ec_total_cur(isnan(ec_total_cur)) = 0;
-
- for i_lau = active_launchers
- ids_core_sources.source{last_index+i_lau}.global_quantities.time = ec_powers_tgrid;
- ids_core_sources.source{last_index+i_lau}.global_quantities.power = ec_total_pow;
- ids_core_sources.source{last_index+i_lau}.global_quantities.current_parallel = ec_total_cur;
- end
-
% load LIUQE data to convert
% \tilde{j}_// = <jdotB>/(R0<Bphi/R>) to j_//0 = <jdotB>/B0; Bphi=T(psi)/R
- [L,~,LY] = liuqe(shot,ec_data_tgrid,'iterq',50,'ilim',3,'icsint',true);
+ [L,~,LY] = liuqe(shot,ec_tgrid,'iterq',50,'ilim',3,'icsint',true);
rho_pol_norm_liu = L.pQ; % normalized sqrt poloidal flux including axis
- mu0 = 4*pi*10^-7;
- Rm2_fs = LY.Q2Q;
- Vprime = -1./LY.Q1Q; % Q1Q: -dpsi/dV
- Ip_psi = LY.ItQ; % Ip(psi) = I_tor(psi)
- T = LY.TQ;
- vol = LY.VQ;
- B2_fs = mu0*Ip_psi.*2.*pi./Vprime + T.^2.*Rm2_fs;
- toray_to_parallel = 1./(T.*Rm2_fs).*B2_fs./R0./B0;
+ mu0 = 4*pi*10^-7;
+ Rm2_fs = LY.Q2Q;
+ mVprime = 1./LY.Q1Q; % Q1Q: -dpsi/dV
+ Ip_psi = LY.ItQ; % Ip(psi) = I_tor(psi)
+ T = LY.TQ;
+ vol = LY.VQ;
+ B2_fs = mu0*Ip_psi./mVprime + T.^2.*Rm2_fs;
+ R0 = ids_core_sources.vacuum_toroidal_field.r0;
+ B0 = interp1(ids_core_sources.time,ids_core_sources.vacuum_toroidal_field.b0,ec_tgrid);
+ jtoray_to_jpar0 = 2*pi./B0./(T.*Rm2_fs).*B2_fs;
% interpolate on ec_data rho_pol grid
- toray_to_parallel_ec = nan(size(toray_to_parallel));
- vol_ec = nan(size(toray_to_parallel));
+ jtoray_to_jpar0_mapped = nan(size(ec_data.cd_dens.x));
+ vol_mapped = nan(size(ec_data.cd_dens.x));
for ii = 1:nt_ec_data
- toray_to_parallel_ec(:,ii) = ...
- interp1(rho_pol_norm_liu,toray_to_parallel(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
- vol_ec(:,ii) = ...
+ jtoray_to_jpar0_mapped(:,ii) = ...
+ interp1(rho_pol_norm_liu,jtoray_to_jpar0(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+ vol_mapped(:,ii) = ...
interp1(rho_pol_norm_liu,vol(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
end
- % interpoating p_dens profiles from 'ec_data_tgrid' grid (toray tgrid) to 'ec_powers_tgrid' (powers tgrid)
+ % interpoating p_dens profiles from 'ec_tgrid' grid (toray tgrid) to 'ec_powers_tgrid' (powers tgrid)
p_dens = ec_data.p_dens.data;
p_integrated = ec_data.p_integrated.data;
- cd_dens = nan(size(ec_data.cd_dens.data));
- cd_integrated = nan(size(ec_data.cd_integrated.data));
- %%
+ cd_dens = ec_data.cd_dens_doublewidth.data; % use double width to have realistic deposition broadening
+ cd_integrated = ec_data.cd_integrated.data;
+ %
for i_lau = active_launchers
- cd_dens(:,i_lau,:) = squeeze(cd_dens(:,i_lau,:)).*toray_to_parallel_ec;
+ % convert to j// with conversion factor
+ cd_dens(:,i_lau,:) = squeeze(cd_dens(:,i_lau,:)).*jtoray_to_jpar0_mapped;
for ii = 1:nt_ec_data
- cd_integrated(:,i_lau,ii) = cumtrapz(vol_ec(:,ii),cd_dens(:,i_lau,ii))/(2*pi*0.88);
+ % integrate j//
+ cd_integrated(:,i_lau,ii) = cumtrapz(vol_mapped(:,ii),cd_dens(:,i_lau,ii))/(2*pi*R0);
end
end
-%%
- p_ec_injected = powers_gdat.ec.data;
- ij = iround_os(ec_powers_tgrid,ec_data_tgrid);
- sparse_p_ec_injected = p_ec_injected(ij,:); % injected ec power vals corresponding to ec_data_tgrid
+
+%% remap ec outputs on power time grid
+
+ ij = iround_os(interp_tgrid,ec_tgrid);
+ sparse_p_ec_injected = interp_p_ec_injected(ij,:); % injected ec power vals corresponding to ec_tgrid
n_rho = size(p_dens, 1);
+ %rho grids to be interpolated on power time grid
+ rho_pol_norm = ec_data.cd_dens.grids.rho_pol_norm;
+ rho_tor_norm = ec_data.cd_dens.grids.rho_tor_norm;
- % calculate normalised profiles on ec_data_tgrid grid
+ % calculate normalised profiles on ec_tgrid grid
norm_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_data);
norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_data);
norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_data);
@@ -228,10 +234,14 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
end
% interpolate normalised p_dens profiles on ec_powers_tgrid
- interp_norm_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers);
- interp_norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers);
- interp_norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers);
- interp_norm_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers);
+ interp_norm_p_dens = zeros(n_rho,nb_launchers+1,nt_interp);
+ interp_norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_interp);
+ interp_norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_interp);
+ interp_norm_cd_integrated = zeros(n_rho,nb_launchers+1,nt_interp);
+ % interpolate the rho_pol & rho_tor on ec_powers_tgrid
+ interp_rho_pol = zeros(n_rho,nt_interp);
+ interp_rho_tor = zeros(n_rho,nt_interp);
+ interp_rho_flag = 1;
for i_lau = active_launchers
for irho = 1:n_rho
% get power and current density at each rho
@@ -240,47 +250,64 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
trace_cd_dens = squeeze(norm_cd_dens(irho,i_lau,:));
trace_cd_integrated = squeeze(norm_cd_integrated(irho,i_lau,:));
% interpolate on gdat powers tgrid
- interp_norm_p_dens(irho,i_lau,:) = interp1(ec_data_tgrid,trace_p_dens,ec_powers_tgrid);
- interp_norm_p_integrated(irho,i_lau,:) = interp1(ec_data_tgrid,trace_p_integrated,ec_powers_tgrid);
- interp_norm_cd_dens(irho,i_lau,:) = interp1(ec_data_tgrid,trace_cd_dens,ec_powers_tgrid);
- interp_norm_cd_integrated(irho,i_lau,:) = interp1(ec_data_tgrid,trace_cd_integrated,ec_powers_tgrid);
+ interp_norm_p_dens(irho,i_lau,:) = interpos(ec_tgrid,trace_p_dens,interp_tgrid);
+ interp_norm_p_integrated(irho,i_lau,:) = interpos(ec_tgrid,trace_p_integrated,interp_tgrid);
+ interp_norm_cd_dens(irho,i_lau,:) = interpos(ec_tgrid,trace_cd_dens,interp_tgrid);
+ interp_norm_cd_integrated(irho,i_lau,:) = interpos(ec_tgrid,trace_cd_integrated,interp_tgrid);
+
+ % only interpolate rho grids once
+ if interp_rho_flag
+ interp_rho_pol(irho,:) = interpos(ec_tgrid,rho_pol_norm(irho,:),interp_tgrid');
+ interp_rho_tor(irho,:) = interpos(ec_tgrid,rho_tor_norm(irho,:),interp_tgrid');
+ end
end
+ interp_rho_flag = 0;
end
- % normalised & interpolated profiles * p_ec_injected on ec_powers_tgrid
- interp_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers);
- interp_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers);
- interp_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers);
- interp_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers);
- for it = 1:nt_ec_powers
- interp_p_dens(:,:,it) = interp_norm_p_dens(:,:,it).*p_ec_injected(it,:);
- interp_p_integrated(:,:,it) = interp_norm_p_integrated(:,:,it).*p_ec_injected(it,:);
- interp_cd_dens(:,:,it) = interp_norm_cd_dens(:,:,it).*p_ec_injected(it,:);
- interp_cd_integrated(:,:,it) = interp_norm_cd_integrated(:,:,it).*p_ec_injected(it,:);
+ % normalised & interpolated profiles * p_ec_injected on interp_tgrid
+ interp_p_dens = zeros(n_rho,nb_launchers+1,nt_interp);
+ interp_p_integrated = zeros(n_rho,nb_launchers+1,nt_interp);
+ interp_cd_dens = zeros(n_rho,nb_launchers+1,nt_interp);
+ interp_cd_integrated = zeros(n_rho,nb_launchers+1,nt_interp);
+ for it = 1:nt_interp
+ interp_p_dens(:,:,it) = interp_norm_p_dens(:,:,it).*interp_p_ec_injected(it,:);
+ interp_p_integrated(:,:,it) = interp_norm_p_integrated(:,:,it).*interp_p_ec_injected(it,:);
+ interp_cd_dens(:,:,it) = interp_norm_cd_dens(:,:,it).*interp_p_ec_injected(it,:);
+ interp_cd_integrated(:,:,it) = interp_norm_cd_integrated(:,:,it).*interp_p_ec_injected(it,:);
end
% fill profiles_1d with interpolated profiles
for i_lau = active_launchers
- for ii = 1:nt_ec_powers
+ % initialize var for current_parallel and coupled power which needs to be obatined from
+ % integrated j// profiles and pdens last entry
+ current_parallel_tmp = nan(1,nt_interp);
+ power_coupled_tmp = nan(1,nt_interp);
+
+ for ii = 1:nt_interp
% time
- ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.time = ec_powers_tgrid(ii);
- % rhotor grid
- ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_tor_norm = ...
- ec_data.p_dens.rhotor_norm(1,:); % constant in time
+ ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.time = interp_tgrid(ii);
+ % grids
+ ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_tor_norm = interp_rho_tor(:,ii);
+ ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_pol_norm = interp_rho_pol(:,ii);
% power density
ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.energy = ...
interp_p_dens(:,i_lau,ii);
% integrated power density
ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.power_inside = ...
interp_p_integrated(:,i_lau,ii);
- % current density (to adapt to <J.B>/B0)
+ power_coupled_tmp(ii) = interp_p_integrated(end,i_lau,ii);
+ % current density
ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.j_parallel = ...
interp_cd_dens(:,i_lau,ii);
- % integrated current density (to adapt to INTEGRAL(<J.B>/B0)*ds_phi)
+ % integrated current density
ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.current_parallel_inside = ...
interp_cd_integrated(:,i_lau,ii);
+ current_parallel_tmp(ii) = interp_cd_integrated(end,i_lau,ii);
end
- end
+ ids_core_sources.source{last_index+i_lau}.globa_quantities.time = interp_tgrid;
+ ids_core_sources.source{last_index+i_lau}.global_quantities.power = power_coupled_tmp;
+ ids_core_sources.source{last_index+i_lau}.global_quantities.current_parallel = current_parallel_tmp;
+ end
%add empty sources for rest of unsused launchers
if numel(ids_core_sources.source)-last_index ~= nb_launchers
@@ -292,22 +319,20 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
end
%% nbi
-id_nbi.description = 'Source from Neutral Beam Injection';
-id_nbi.index = 2; id_nbi.name = 'nbi';
-%ids_core_sources.source{} = source_template;
-%ids_core_sources.source{}.identifier = id_nbi;
-%ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template};
-%ids_core_sources.source{}.global_quantities(1:n_t) = {globals_template};
-
-
-% total
-id_total.description = 'Total source; combines all sources';
-id_total.index = 1; id_total.name = 'total';
-%ids_core_sources.source{} = source_template;
-%ids_core_sources.source{}.identifier = id_total;
-%ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template};
-%ids_core_sources.source{}.global_quantities(1:n_t) = {globals_template};
-
+% id_nbi.description = 'Source from Neutral Beam Injection';
+% id_nbi.index = 2; id_nbi.name = 'nbi';
+% ids_core_sources.source{} = source_template;
+% ids_core_sources.source{}.identifier = id_nbi;
+% ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template};
+% ids_core_sources.source{}.global_quantities(1:n_t) = {globals_template};
+
+%% total
+% id_total.description = 'Total source; combines all sources';
+% id_total.index = 1; id_total.name = 'total';
+% ids_core_sources.source{} = source_template;
+% ids_core_sources.source{}.identifier = id_total;
+% ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template};
+% ids_core_sources.source{}.global_quantities(1:n_t) = {globals_template};
%% add descriptions for profiles_1d
--
GitLab