Skip to content
Snippets Groups Projects
Commit 87bccfb3 authored by Antonia Frank's avatar Antonia Frank
Browse files

Finalize the EC section with correct conversion factor & integration, reduce... 

Finalize the EC section with correct conversion factor & integration, reduce gdat_powers grid to times when waves are injected to the plasma
parent c173fa38
No related branches found
No related tags found
1 merge request!137Add quantities to ids for MRE
Hide whitespace changes
Inline Side-by-side
matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+ 114
− 89
View file @ 87bccfb3
...@@ -52,7 +52,8 @@ powers_gdat = gdat(params_core_sources.shot,params_eff); ...@@ -52,7 +52,8 @@ powers_gdat = gdat(params_core_sources.shot,params_eff);
% ohm % ohm
ohm_t_grid = powers_gdat.ohm.t; ohm_n_t = numel(ohm_t_grid); 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'; id_ohm.index = 7; id_ohm.name = 'ohmic';
ids_core_sources.source{last_index+1} = source_template; ids_core_sources.source{last_index+1} = source_template;
ids_core_sources.source{last_index+1}.identifier = id_ohm; ids_core_sources.source{last_index+1}.identifier = id_ohm;
...@@ -106,7 +107,8 @@ if any(bs_t_grid-ohm_t_grid')>0 ...@@ -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.') warning('Bootstrap and ohmic time grids are not the same! Interpolation needed.')
end 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'; 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} = source_template;
ids_core_sources.source{last_index+1}.identifier = id_bs; 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 ...@@ -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 last_index = last_index+1; % add if statement to only increment if bs source has been added
%% ec %% 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 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; ec_data = ec_gdat.ec.ec_data; ec_inputs = ec_gdat.ec.ec_inputs;
% get tgrid from gdat ec_data % get tgrid from gdat ec_data
ec_data_tgrid = ec_data.p_dens.t; nt_ec_data = numel(ec_data_tgrid); ec_tgrid = ec_data.p_dens.t; nt_ec_data = numel(ec_tgrid);
% get tgrid from gdat powers
ec_powers_tgrid = powers_gdat.ec.t; nt_ec_powers = numel(ec_powers_tgrid);
% retrieve active launcher information from ec_inputs % retrieve active launcher information from ec_inputs
nb_launchers = numel(ec_inputs.launchers_active.data); nb_launchers = numel(ec_inputs.launchers_active.data);
active_launchers = find(ec_inputs.launchers_active.data==1)'; 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 % Setup structures for active launchers from template
main_desc = 'Source from electron cyclotron heating and current drive'; 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 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} = source_template;
ids_core_sources.source{last_index+i_lau}.identifier = id_ec; 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; ids_core_sources.source{last_index+i_lau}.global_quantities = globals_template;
end 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 % load LIUQE data to convert
% \tilde{j}_// = <jdotB>/(R0<Bphi/R>) to j_//0 = <jdotB>/B0; Bphi=T(psi)/R % \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 rho_pol_norm_liu = L.pQ; % normalized sqrt poloidal flux including axis
mu0 = 4*pi*10^-7; mu0 = 4*pi*10^-7;
Rm2_fs = LY.Q2Q; Rm2_fs = LY.Q2Q;
Vprime = -1./LY.Q1Q; % Q1Q: -dpsi/dV mVprime = 1./LY.Q1Q; % Q1Q: -dpsi/dV
Ip_psi = LY.ItQ; % Ip(psi) = I_tor(psi) Ip_psi = LY.ItQ; % Ip(psi) = I_tor(psi)
T = LY.TQ; T = LY.TQ;
vol = LY.VQ; vol = LY.VQ;
B2_fs = mu0*Ip_psi.*2.*pi./Vprime + T.^2.*Rm2_fs; B2_fs = mu0*Ip_psi./mVprime + T.^2.*Rm2_fs;
toray_to_parallel = 1./(T.*Rm2_fs).*B2_fs./R0./B0; 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 % interpolate on ec_data rho_pol grid
toray_to_parallel_ec = nan(size(toray_to_parallel)); jtoray_to_jpar0_mapped = nan(size(ec_data.cd_dens.x));
vol_ec = nan(size(toray_to_parallel)); vol_mapped = nan(size(ec_data.cd_dens.x));
for ii = 1:nt_ec_data for ii = 1:nt_ec_data
toray_to_parallel_ec(:,ii) = ... jtoray_to_jpar0_mapped(:,ii) = ...
interp1(rho_pol_norm_liu,toray_to_parallel(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii)); interp1(rho_pol_norm_liu,jtoray_to_jpar0(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
vol_ec(:,ii) = ... vol_mapped(:,ii) = ...
interp1(rho_pol_norm_liu,vol(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii)); interp1(rho_pol_norm_liu,vol(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
end 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_dens = ec_data.p_dens.data;
p_integrated = ec_data.p_integrated.data; p_integrated = ec_data.p_integrated.data;
cd_dens = nan(size(ec_data.cd_dens.data)); cd_dens = ec_data.cd_dens_doublewidth.data; % use double width to have realistic deposition broadening
cd_integrated = nan(size(ec_data.cd_integrated.data)); cd_integrated = ec_data.cd_integrated.data;
%% %
for i_lau = active_launchers 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 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
end end
%%
p_ec_injected = powers_gdat.ec.data; %% remap ec outputs on power time grid
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 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); 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_dens = zeros(n_rho,nb_launchers+1,nt_ec_data);
norm_p_integrated = 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); 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 ...@@ -228,10 +234,14 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
end end
% interpolate normalised p_dens profiles on ec_powers_tgrid % 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_dens = zeros(n_rho,nb_launchers+1,nt_interp);
interp_norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers); interp_norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_interp);
interp_norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers); interp_norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_interp);
interp_norm_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers); 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 i_lau = active_launchers
for irho = 1:n_rho for irho = 1:n_rho
% get power and current density at each 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 ...@@ -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_dens = squeeze(norm_cd_dens(irho,i_lau,:));
trace_cd_integrated = squeeze(norm_cd_integrated(irho,i_lau,:)); trace_cd_integrated = squeeze(norm_cd_integrated(irho,i_lau,:));
% interpolate on gdat powers tgrid % 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_dens(irho,i_lau,:) = interpos(ec_tgrid,trace_p_dens,interp_tgrid);
interp_norm_p_integrated(irho,i_lau,:) = interp1(ec_data_tgrid,trace_p_integrated,ec_powers_tgrid); interp_norm_p_integrated(irho,i_lau,:) = interpos(ec_tgrid,trace_p_integrated,interp_tgrid);
interp_norm_cd_dens(irho,i_lau,:) = interp1(ec_data_tgrid,trace_cd_dens,ec_powers_tgrid); interp_norm_cd_dens(irho,i_lau,:) = interpos(ec_tgrid,trace_cd_dens,interp_tgrid);
interp_norm_cd_integrated(irho,i_lau,:) = interp1(ec_data_tgrid,trace_cd_integrated,ec_powers_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 end
interp_rho_flag = 0;
end end
% normalised & interpolated profiles * p_ec_injected on ec_powers_tgrid % normalised & interpolated profiles * p_ec_injected on interp_tgrid
interp_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers); interp_p_dens = zeros(n_rho,nb_launchers+1,nt_interp);
interp_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers); interp_p_integrated = zeros(n_rho,nb_launchers+1,nt_interp);
interp_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers); interp_cd_dens = zeros(n_rho,nb_launchers+1,nt_interp);
interp_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers); interp_cd_integrated = zeros(n_rho,nb_launchers+1,nt_interp);
for it = 1:nt_ec_powers for it = 1:nt_interp
interp_p_dens(:,:,it) = interp_norm_p_dens(:,:,it).*p_ec_injected(it,:); interp_p_dens(:,:,it) = interp_norm_p_dens(:,:,it).*interp_p_ec_injected(it,:);
interp_p_integrated(:,:,it) = interp_norm_p_integrated(:,:,it).*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).*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).*p_ec_injected(it,:); interp_cd_integrated(:,:,it) = interp_norm_cd_integrated(:,:,it).*interp_p_ec_injected(it,:);
end end
% fill profiles_1d with interpolated profiles % fill profiles_1d with interpolated profiles
for i_lau = active_launchers 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 % time
ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.time = ec_powers_tgrid(ii); ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.time = interp_tgrid(ii);
% rhotor grid % grids
ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_tor_norm = ... ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_tor_norm = interp_rho_tor(:,ii);
ec_data.p_dens.rhotor_norm(1,:); % constant in time ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_pol_norm = interp_rho_pol(:,ii);
% power density % power density
ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.energy = ... ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.energy = ...
interp_p_dens(:,i_lau,ii); interp_p_dens(:,i_lau,ii);
% integrated power density % integrated power density
ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.power_inside = ... ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.power_inside = ...
interp_p_integrated(:,i_lau,ii); 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 = ... ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.j_parallel = ...
interp_cd_dens(:,i_lau,ii); 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 = ... ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.current_parallel_inside = ...
interp_cd_integrated(:,i_lau,ii); interp_cd_integrated(:,i_lau,ii);
current_parallel_tmp(ii) = interp_cd_integrated(end,i_lau,ii);
end 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 %add empty sources for rest of unsused launchers
if numel(ids_core_sources.source)-last_index ~= nb_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 ...@@ -292,22 +319,20 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
end end
%% nbi %% nbi
id_nbi.description = 'Source from Neutral Beam Injection'; % id_nbi.description = 'Source from Neutral Beam Injection';
id_nbi.index = 2; id_nbi.name = 'nbi'; % id_nbi.index = 2; id_nbi.name = 'nbi';
%ids_core_sources.source{} = source_template; % ids_core_sources.source{} = source_template;
%ids_core_sources.source{}.identifier = id_nbi; % ids_core_sources.source{}.identifier = id_nbi;
%ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template}; % ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template};
%ids_core_sources.source{}.global_quantities(1:n_t) = {globals_template}; % ids_core_sources.source{}.global_quantities(1:n_t) = {globals_template};
%% total
% total % id_total.description = 'Total source; combines all sources';
id_total.description = 'Total source; combines all sources'; % id_total.index = 1; id_total.name = 'total';
id_total.index = 1; id_total.name = 'total'; % ids_core_sources.source{} = source_template;
%ids_core_sources.source{} = source_template; % ids_core_sources.source{}.identifier = id_total;
%ids_core_sources.source{}.identifier = id_total; % ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template};
%ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template}; % ids_core_sources.source{}.global_quantities(1:n_t) = {globals_template};
%ids_core_sources.source{}.global_quantities(1:n_t) = {globals_template};
%% add descriptions for profiles_1d %% add descriptions for profiles_1d
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment