From d3d59d95abac21dad681e12c478f530e367a5e79 Mon Sep 17 00:00:00 2001
From: Antonia Frank <antonia.frank@epfl.ch>
Date: Fri, 9 Aug 2024 18:40:21 +0200
Subject: [PATCH] Update ec_tgrid_out for interpolation saveing to ids
---
matlab/TCV_IMAS/tcv_get_ids_core_sources.m | 134 +++++++++++----------
1 file changed, 69 insertions(+), 65 deletions(-)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index a904007c..ea30cb6a 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -58,7 +58,7 @@ 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;
ids_core_sources.source{last_index+1}.profiles_1d(1:ohm_n_t) = {profiles_template};
-ids_core_sources.source{last_index+1}.global_quantities(1:ohm_n_t) = globals_template;
+ids_core_sources.source{last_index+1}.global_quantities(1:ohm_n_t) = {globals_template};
ohm_data = ohm_gdat.ohm.ohm_data;
% load LIUQE data to convert
@@ -107,7 +107,7 @@ 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;
ids_core_sources.source{last_index+1}.profiles_1d(1:bs_n_t) = {profiles_template};
-ids_core_sources.source{last_index+1}.global_quantities(1:bs_n_t) = globals_template;
+ids_core_sources.source{last_index+1}.global_quantities(1:bs_n_t) = {globals_template};
for ii = 1:bs_n_t
% profiles_1d
@@ -139,7 +139,7 @@ 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_tgrid = ec_data.p_dens.t; nt_ec_data = numel(ec_tgrid);
+ ec_tgrid_toray = ec_data.p_dens.t; nt_ec_toray = numel(ec_tgrid_toray);
% retrieve active launcher information from ec_inputs
nb_launchers = numel(ec_inputs.launchers_active.data);
@@ -147,10 +147,13 @@ 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);
- 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);
+ 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);
% Setup source structs for active launchers from template
main_desc = 'Source from electron cyclotron heating and current drive';
@@ -160,13 +163,13 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
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_interp) = {profiles_template};
- ids_core_sources.source{last_index+i_lau}.global_quantities(1:nt_interp) = globals_template;
+ ids_core_sources.source{last_index+i_lau}.profiles_1d(1:nt_ec_out) = {profiles_template};
+ ids_core_sources.source{last_index+i_lau}.global_quantities(1:nt_ec_out) = {globals_template};
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_tgrid,'iterq',50,'ilim',3,'icsint',true);
+ [L,~,LY] = liuqe(shot,ec_tgrid_toray,'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;
@@ -176,13 +179,13 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
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);
+ B0 = interp1(ids_core_sources.time,ids_core_sources.vacuum_toroidal_field.b0,ec_tgrid_toray);
jtoray_to_jpar0 = 2*pi./B0./(T.*Rm2_fs).*B2_fs;
% interpolate on ec_data rho_pol grid
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
+ for ii = 1:nt_ec_toray
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) = ...
@@ -198,41 +201,37 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
for i_lau = active_launchers
% 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_toray
% integrate j//
cd_integrated(:,i_lau,ii) = cumtrapz(vol_mapped(:,ii),cd_dens(:,i_lau,ii))/(2*pi*R0);
end
end
-%% interpolating 1d_profiles from 'ec_tgrid' grid (toray tgrid) to 'interp_tgrid' (powers tgrid when powers.data>0)
- ij = iround_os(interp_tgrid,ec_tgrid);
- sparse_p_ec_injected = interp_p_ec_injected(ij,:); % injected ec power vals corresponding to ec_tgrid
+%% interpolating 1d_profiles from 'ec_tgrid_toray' grid (toray tgrid) to 'interp_tgrid' (powers tgrid when powers.data>0)
+ ij = iround_os(ec_tgrid_out,ec_tgrid_toray);
+ sparse_p_ec_injected = p_ec_injected(ij); % injected ec power vals corresponding to ec_tgrid_toray
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_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);
- norm_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_data);
- for it = 1:nt_ec_data
- norm_p_dens(:,:,it) = p_dens(:,:,it)./sparse_p_ec_injected(it,:);
- norm_p_integrated(:,:,it) = p_integrated(:,:,it)./sparse_p_ec_injected(it,:);
- norm_cd_dens(:,:,it) = cd_dens(:,:,it)./sparse_p_ec_injected(it,:);
- norm_cd_integrated(:,:,it) = cd_integrated(:,:,it)./sparse_p_ec_injected(it,:);
+ % calculate normalised profiles on ec_tgrid_toray grid
+ norm_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_toray);
+ norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_toray);
+ norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_toray);
+ norm_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_toray);
+ for it = 1:nt_ec_toray
+ norm_p_dens(:,:,it) = p_dens(:,:,it) ./sparse_p_ec_injected(it);
+ norm_p_integrated(:,:,it) = p_integrated(:,:,it) ./sparse_p_ec_injected(it);
+ norm_cd_dens(:,:,it) = cd_dens(:,:,it) ./sparse_p_ec_injected(it);
+ norm_cd_integrated(:,:,it) = cd_integrated(:,:,it)./sparse_p_ec_injected(it);
end
% interpolate normalised p_dens profiles on ec_powers_tgrid
- 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;
+ interp_norm_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ interp_norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ interp_norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ interp_norm_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
for i_lau = active_launchers
for irho = 1:n_rho
% get power and current density at each rho
@@ -241,59 +240,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,:) = 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
+ interp_norm_p_dens( irho,i_lau,3:end-2) = interpos(ec_tgrid_toray,trace_p_dens, ec_tgrid_out(3:end-2));
+ interp_norm_p_integrated( irho,i_lau,3:end-2) = interpos(ec_tgrid_toray,trace_p_integrated, ec_tgrid_out(3:end-2));
+ interp_norm_cd_dens( irho,i_lau,3:end-2) = interpos(ec_tgrid_toray,trace_cd_dens, ec_tgrid_out(3:end-2));
+ interp_norm_cd_integrated(irho,i_lau,3:end-2) = interpos(ec_tgrid_toray,trace_cd_integrated,ec_tgrid_out(3:end-2));
end
- interp_rho_flag = 0;
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);
+ for irho = 1:n_rho
+ 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
+ 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
- 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,:);
+ interp_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ interp_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ interp_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ interp_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ for it = 1:nt_ec_out % fill profiles only for time slices when EC is nonzero, thus it+2 assignment
+ 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);
end
% fill the ids entries
for i_lau = active_launchers
- for ii = 1:nt_interp
+ for ii = 1:nt_ec_out
% 1d_profiles
- 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);
+ ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.time = ec_tgrid_out(ii);
+ % grids (do not fill for first and last time slice)
+ 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);
- power_coupled_tmp = 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
ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.current_parallel_inside = ...
interp_cd_integrated(:,i_lau,ii);
- current_parallel_tmp = interp_cd_integrated(end,i_lau,ii);
-
+ end
+
% globals
- ids_core_sources.source{last_index+i_lau}.globa_quantities{ii}.time = interp_tgrid(ii);
- ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.power = power_coupled_tmp;
- ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.current_parallel = current_parallel_tmp;
+ 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);
end
end
--
GitLab