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

icdbseval version added to test; with plotting to compare

parent 0f8108c4
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
+ 83
− 35
View file @ e569476b
......@@ -208,7 +208,7 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
% load LIUQE data to convert
% \tilde{V,TORAY}_// = dI_\phi/dV = (1/2\pi)<j_\phi/R> to j_//0 = <jdotB>/B0;
% conversion using j_//,cd = <j_cd/B> B, thus j_// = 2\pi <B^2> / (B_0 <B_\phi/R>)
% conversion using j_//,cd = <j_cd/B> B, thus j_// = 2\pi <B^2> / (B_0 <B_\phi/R>)
mu0 = 4*pi*10^-7;
% interpolate liuqe outputs on ohm_tgrid
T = interp1(liuqe_time,T_liu.data.', ec_tgrid_toray)';
......@@ -216,14 +216,17 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
vol = interp1(liuqe_time,vol_liu.data.', ec_tgrid_toray)';
Ip = interp1(liuqe_time,Ip_liu.data, ec_tgrid_toray);
q = interp1(liuqe_time,q_liu.data.', ec_tgrid_toray)';
mVprime = -2*pi*q./T./Rm2_fs;
B2_fs = mu0*Ip./mVprime + T.^2.*Rm2_fs;
Vprime = -2*pi*q./T./Rm2_fs;
B2_fs = -mu0*Ip./Vprime + T.^2.*Rm2_fs;
% Iphi? needs to be a profile?
% R is also a profile in icdbseval? should not be! It's definitely R0 see
% eq.(30) in dok
% get vacuum field data from ids
R0 = ids_core_sources.vacuum_toroidal_field.r0;
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;
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));
......@@ -233,6 +236,10 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
interp1(rho_pol_norm_liu.data,jtoray_to_jpar0(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
vol_mapped(:,ii) = ...
interp1(rho_pol_norm_liu.data,vol(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
% for testing
T_mapped(:,ii) = interp1(rho_pol_norm_liu.data,T(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
[~,dTdV_mapped(:,ii)] = interpos(vol_mapped(:,ii),T_mapped(:,ii));
Rm2_fs_mapped(:,ii) = interp1(rho_pol_norm_liu.data,Rm2_fs(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
end
% load 1d_profiles from ec_data
......@@ -242,17 +249,58 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
p_integrated(isnan(p_integrated)) = 0;
cd_dens = ec_data.cd_dens_doublewidth.data; % use double width to have realistic deposition broadening
cd_dens(isnan(cd_dens)) = 0;
% initialize variables
jpar0 = zeros(size(cd_dens));
jpar0_integrated = zeros(size(cd_dens));
cd_integrated = ec_data.cd_integrated.data;
cd_integrated = zeros(size(ec_data.cd_integrated.data));
for i_lau = active_launchers
% convert to j// with conversion factor
cd_dens(:,i_lau,:) = squeeze(cd_dens(:,i_lau,:)).*jtoray_to_jpar0_mapped;
jpar0(:,i_lau,:) = squeeze(cd_dens(:,i_lau,:)).*jtoray_to_jpar0_mapped;
jpar_tilde(:,i_lau,:) = 2*pi/R0./Rm2_fs_mapped.*(squeeze(cd_dens(:,i_lau,:))-dTdV_mapped./T_mapped.*squeeze(cd_integrated(:,i_lau,:)));
jpar02(:,i_lau,:) = squeeze(jpar_tilde(:,i_lau,:)).*R0.*T_mapped.*Rm2_fs_mapped./B0;
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);
[~,~,~,jpar0_integrated(:,i_lau,ii)] = interpos(vol_mapped(:,ii),jpar0(:,i_lau,ii)/(2*pi*R0));
[~,~,~,jpar02_integrated(:,i_lau,ii)] = interpos(vol_mapped(:,ii),jpar02(:,i_lau,ii)/(2*pi*R0));
end
end
figure; hold on; plot(squeeze(jpar0(:,1,50)));
plot(squeeze(jpar02(:,1,50)),'x');
legend('ids','icdbseval'); ylabel('j||0')
figure; plot(squeeze(cd_integrated(end,1,:))); hold on;
plot(squeeze(jpar0_integrated(end,1,:)));
plot(squeeze(jpar02_integrated(end,1,:)));
legend('I_\phi','I||,ids','I||,icdbseval');
% test to integrate but still is not the same
ilau = 4;
jtilde_test = B0./R0./T_mapped./Rm2_fs_mapped.*squeeze(jpar0(:,ilau,:));
jtilde_test2 = squeeze(jpar_tilde(:,ilau,:));
for ii = 1:nt_ec_toray
[~,~,~,icd_test(:,ii)] = interpos(vol_mapped(:,ii),Rm2_fs_mapped(:,ii)./T_mapped(:,ii).*jtilde_test(:,ii));
out(:,ii) = R0/2/pi.*T_mapped(:,ii).*icd_test(:,ii);
% icdbseval
[~,~,~,icd_test2(:,ii)] = interpos(vol_mapped(:,ii),Rm2_fs_mapped(:,ii)./T_mapped(:,ii).*jtilde_test2(:,ii));
out2(:,ii) = R0/2/pi.*T_mapped(:,ii).*icd_test2(:,ii);
end
figure;
plot(ec_tgrid_toray,out(end,:)); hold on;
plot(ec_tgrid_toray,squeeze(cd_integrated(end,ilau,:)));
plot(ec_tgrid_toray,out2(end,:),'x');
ylabel('I_\phi'); legend('from ids','from icdbseval','TORAY')
figure;
plot(squeeze(out(end,:))-squeeze(cd_integrated(end,1,:))'); hold on; plot(squeeze(out2(end,:))-squeeze(cd_integrated(end,1,:))');
ylabel('diff (I_\phi,I_TORARY)'); legend('from ids','from icdbseval')
%% 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
......@@ -262,35 +310,35 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
rho_tor_norm = ec_data.cd_dens.grids.rho_tor_norm;
% 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);
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_jpar0 = zeros(n_rho,nb_launchers+1,nt_ec_toray);
norm_jpar0_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);
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_jpar0(:,:,it) = jpar0(:,:,it) ./sparse_p_ec_injected(it);
norm_jpar0_integrated(:,:,it) = jpar0_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_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);
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_jpar0 = zeros(n_rho,nb_launchers+1,nt_ec_out);
interp_norm_jpar0_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
trace_p_dens = squeeze(norm_p_dens(irho,i_lau,:));
trace_p_integrated = squeeze(norm_p_integrated(irho,i_lau,:));
trace_cd_dens = squeeze(norm_cd_dens(irho,i_lau,:));
trace_cd_integrated = squeeze(norm_cd_integrated(irho,i_lau,:));
trace_jpar0 = squeeze(norm_jpar0(irho,i_lau,:));
trace_jpar0_integrated = squeeze(norm_jpar0_integrated(irho,i_lau,:));
% interpolate on gdat powers tgrid, interpos 21 to extrapolate,
% taking constant value y_edge; needed for time window between start of EC and first TS time
interp_norm_p_dens( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_p_dens, ec_tgrid_out(3:end-2));
interp_norm_p_integrated( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_p_integrated, ec_tgrid_out(3:end-2));
interp_norm_cd_dens( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_cd_dens, ec_tgrid_out(3:end-2));
interp_norm_cd_integrated(irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_cd_integrated,ec_tgrid_out(3:end-2));
interp_norm_p_dens( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_p_dens, ec_tgrid_out(3:end-2));
interp_norm_p_integrated( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_p_integrated, ec_tgrid_out(3:end-2));
interp_norm_jpar0( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_jpar0, ec_tgrid_out(3:end-2));
interp_norm_jpar0_integrated(irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_jpar0_integrated,ec_tgrid_out(3:end-2));
end
end
......@@ -307,13 +355,13 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
% normalised & interpolated profiles * p_ec_injected on interp_tgrid
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);
interp_jpar0 = zeros(n_rho,nb_launchers+1,nt_ec_out);
interp_jpar0_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);
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_jpar0(:,:,it) = interp_norm_jpar0(:,:,it) .*p_ec_injected(it);
interp_jpar0_integrated(:,:,it) = interp_norm_jpar0_integrated(:,:,it).*p_ec_injected(it);
end
% fill the ids entries
......@@ -333,10 +381,10 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
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);
interp_jpar0(:,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);
interp_jpar0_integrated(:,i_lau,ii);
end
% globals
......@@ -344,7 +392,7 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
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_jpar0_integrated(end,i_lau,ii);
end
end
......
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