From 73f02bb33ccca13c4ada3dfc130a59bd1b9ec72d Mon Sep 17 00:00:00 2001
From: Antonia Frank <antonia.frank@epfl.ch>
Date: Tue, 13 Aug 2024 12:07:58 +0200
Subject: [PATCH] Add implementation and check for ASTAR file in Lac8D for NBI
---
matlab/TCV_IMAS/tcv_get_ids_core_sources.m | 128 ++++++++++++++++-----
1 file changed, 102 insertions(+), 26 deletions(-)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index ea30cb6a..afeb88d1 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -124,7 +124,7 @@ for ii = 1:bs_n_t
% globals
ids_core_sources.source{last_index+1}.global_quantities{ii}.time = bs_t_grid(ii);
- ids_core_sources.source{last_index+1}.global_quantities{ii}.current_parallel = integrated_jpar_0_tmp(end);
+ ids_core_sources.source{last_index+1}.global_quantities{ii}.current_parallel = integrated_jpar_0_tmp(end);
end
last_index = last_index+1; % add if statement to only increment if bs source has been added
@@ -147,10 +147,10 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
% find times of injected EC power to interpolate power & current densities
% on p_ec_injected tgrid
- ec_powers_tgrid = powers_gdat.ec.t; %nt_ec_powers = numel(ec_powers_tgrid);
- % find times where EC is on to define time grid with extra time slice just
+ ec_powers_tgrid = powers_gdat.ec.t; %nt_ec_powers = numel(ec_powers_tgrid);
+ % find times where EC is on to define time grid with extra time slice just
% before/after EC power and at start/end of shot
- itime_ec = find(powers_gdat.ec.data(:,end)>0);
+ itime_ec = find(powers_gdat.ec.data(:,end)>0);
ec_tgrid_out = [ohm_t_grid(1),ec_powers_tgrid(itime_ec(1)-1:itime_ec(end)+1)',ohm_t_grid(end)];
nt_ec_out = numel(ec_tgrid_out);
p_ec_injected = interpos(ec_powers_tgrid,powers_gdat.ec.data(:,end),ec_tgrid_out);
@@ -246,7 +246,7 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
interp_norm_cd_integrated(irho,i_lau,3:end-2) = interpos(ec_tgrid_toray,trace_cd_integrated,ec_tgrid_out(3:end-2));
end
end
-
+
% interpolate the rho_pol & rho_tor on ec_powers_tgrid
interp_rho_pol = zeros(n_rho,nt_ec_out);
interp_rho_tor = zeros(n_rho,nt_ec_out);
@@ -254,7 +254,7 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
interp_rho_pol(irho,3:end-2) = interpos(ec_tgrid_toray,rho_pol_norm(irho,:),ec_tgrid_out(3:end-2));
interp_rho_tor(irho,3:end-2) = interpos(ec_tgrid_toray,rho_tor_norm(irho,:),ec_tgrid_out(3:end-2));
end
- % fill rho_pol just outside of TORAY times(zero power) with the known profiles one time slice after/before
+ % fill rho_pol just outside of TORAY times(zero power) with the known profiles one time slice after/before
interp_rho_pol(:,2) = interp_rho_pol(:,3); interp_rho_pol(:,end-1) = interp_rho_pol(:,end-2);
% normalised & interpolated profiles * p_ec_injected on interp_tgrid
@@ -278,30 +278,30 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
if ii~=1 || ii~=nt_ec_out
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
- ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.j_parallel = ...
- interp_cd_dens(:,i_lau,ii);
- % integrated current density
- ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.current_parallel_inside = ...
- interp_cd_integrated(:,i_lau,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
+ ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.j_parallel = ...
+ interp_cd_dens(:,i_lau,ii);
+ % integrated current density
+ ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.current_parallel_inside = ...
+ interp_cd_integrated(:,i_lau,ii);
end
-
+
% globals
ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.time = ec_tgrid_out(ii);
ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.power = ...
interp_p_integrated(end,i_lau,ii);
ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.current_parallel = ...
- interp_cd_integrated(end,i_lau,ii);
+ interp_cd_integrated(end,i_lau,ii);
end
end
- %add empty sources for rest of unsused launchers
+ %add empty sources for rest of unused launchers
if numel(ids_core_sources.source)-last_index ~= nb_launchers
ids_core_sources.source{last_index+nb_launchers} = [];
end
@@ -309,13 +309,89 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
last_index = last_index+nb_launchers;
end
-%% nbi
+%% NBI. same tgrid for NBI1 and NBI2
+
+check_nbi = [~isempty(powers_gdat.nbi1.data),~isempty(powers_gdat.nbi2.data)];
+active_nbi = find(check_nbi==1);
+nb_nbi = numel(check_nbi);
+nbi_names = {'nbi1','nbi2'};
+
+if numel(active_nbi)>0
+ % get tgrid (same for NBI1 & 2 if both active) from active_nbi(1)
+ nbi_powers_tgrid = powers_gdat.(nbi_names{active_nbi(1)}).t;
+ % find times where NBI is on to define time grid with extra time slice just
+ % before/after NBI power and at start/end of shot
+ itime_nbi = find((powers_gdat.(nbi_names{active_nbi(1)}).data>0));% & ...
+ % powers_gdat.(nbi_names{active_nbi(1)}).t<(ohm_t_grid(end)));
+ if nbi_powers_tgrid(itime_nbi(end))>=ohm_t_grid(end)
+ i_time_end = iround(nbi_powers_tgrid,ohm_t_grid(end));
+ nbi_tgrid_out = [ohm_t_grid(1),nbi_powers_tgrid(itime_nbi(1)-1:i_time_end)'];
+ else
+ nbi_tgrid_out = [ohm_t_grid(1),nbi_powers_tgrid(itime_nbi(1)-1:itime_nbi(end)+1)',ohm_t_grid(end)];
+ end
+ nt_nbi_out = numel(nbi_tgrid_out);
+
+ % Setup source structs for active nbi from template
+ main_desc = 'Source from Neutral Beam Injection';
+ id_nbi.index = 2; id_nbi.name = 'nbi';
+ for i_nbi = active_nbi
+ id_nbi.description = sprintf('NBI%i %s',i_nbi,main_desc);
+ ids_core_sources.source{last_index+i_nbi} = source_template;
+ ids_core_sources.source{last_index+i_nbi}.identifier = id_nbi;
+ ids_core_sources.source{last_index+i_nbi}.profiles_1d(1:nt_nbi_out) = {profiles_template};
+ ids_core_sources.source{last_index+i_nbi}.global_quantities(1:nt_nbi_out) = {globals_template};
+ end
+
+ sprintf('Loading of current & power densities from ASTRA not implemented yet. \n')
+ sprintf('Checking if ASTRA run available on partition /Lac8_D:')
+ [~,hostname] = unix('hostname');
+ if strcmp(hostname,'lac8.epfl.ch')
+ unix(sprintf('ls /Lac8_D/ASTRA/ | grep ''%i'' && echo File for shotnumber exists! || echo File for shotnumber does not exist!',shot));
+ else
+ unix(sprintf('ssh $(whoami)@lac8 "ls /Lac8_D/ASTRA/ | grep ''%i'' && echo File for shotnumber exists! || echo File for shotnumber does not exist! && exit"',shot));
+ end
+
+ for i_nbi = active_nbi
+ p_nbi_injected_tmp = interpos(nbi_powers_tgrid,powers_gdat.(nbi_names{i_nbi}).data,nbi_tgrid_out);
+ for ii = 1:nt_nbi_out
+ % globals
+ ids_core_sources.source{last_index+i_nbi}.global_quantities{ii}.time = nbi_tgrid_out(ii);
+ ids_core_sources.source{last_index+i_nbi}.global_quantities{ii}.power = p_nbi_injected_tmp(ii);
+ end
+ end
+
+ % add empty for unused NBI
+ if numel(ids_core_sources.source)-last_index ~= nb_nbi
+ ids_core_sources.source{last_index+nb_nbi} = [];
+ end
+
+ last_index = last_index+nb_nbi;
+end
+
+%% DNBI has it's own time grid
+% if ~isempty(powers_gdat.nbi1)
+%
% 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};
+% ids_core_sources.source{last_index+1} = source_template;
+% ids_core_sources.source{last_index+1}.identifier = id_nbi;
+% ids_core_sources.source{last_index+1}.profiles_1d(1:n_t) = {profiles_template};
+% ids_core_sources.source{last_index+1}.global_quantities(1:n_t) = {globals_template};
+%
+% % find times of injected EC power to interpolate power & current densities
+% % on p_ec_injected tgrid
+% ec_powers_tgrid = powers_gdat.ec.t; %nt_ec_powers = numel(ec_powers_tgrid);
+% % find times where EC is on to define time grid with extra time slice just
+% % before/after EC power and at start/end of shot
+% itime_ec = find(powers_gdat.ec.data(:,end)>0);
+% ec_tgrid_out = [ohm_t_grid(1),ec_powers_tgrid(itime_ec(1)-1:itime_ec(end)+1)',ohm_t_grid(end)];
+% nt_ec_out = numel(ec_tgrid_out);
+% p_ec_injected = interpos(ec_powers_tgrid,powers_gdat.ec.data(:,end),ec_tgrid_out);
+%
+% last_index = last_index+nb_launchers;
+% end
+
+
%% total
% id_total.description = 'Total source; combines all sources';
--
GitLab