From 6a1fb1b5d0f1704eb04d21a4d2fe12347f9c5e50 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Wed, 27 Mar 2019 14:47:02 +0000
Subject: [PATCH] further debug aug equil
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@11654 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/AUG_IMAS/aug_get_ids_equilibrium.m | 111 ++++++++++++---------
1 file changed, 63 insertions(+), 48 deletions(-)
diff --git a/crpptbx/AUG_IMAS/aug_get_ids_equilibrium.m b/crpptbx/AUG_IMAS/aug_get_ids_equilibrium.m
index 077554d3..5fb8c79d 100644
--- a/crpptbx/AUG_IMAS/aug_get_ids_equilibrium.m
+++ b/crpptbx/AUG_IMAS/aug_get_ids_equilibrium.m
@@ -15,18 +15,29 @@ machine = 'aug';
% vacuum_toroidal_field and time, using homogeneous and q_rho time base
%
% time may be "too long" in shotfile compared to effective equilibrium reconstructed times
-temp.q_rho = gdat(params_equilibrium.shot,'q_rho','machine',machine);
-temp_desc.q_rho = 'q_rho';
-ids_equilibrium.time = temp.q_rho.t;
-ids_equilibrium_description.time = '.t subfield from: [''q_rho'']';
+
+% use "equil" for AUG as reference at this stage, it contains most info, so no need for individual calls
+temp_2d.equil = gdat(params_equilibrium.shot,'equil','machine',machine);
+R0exp = 1.65; % to decide B0 vacuum from Jpol
+
+% temp.q_rho = gdat(params_equilibrium.shot,'q_rho','machine',machine);
+ids_equilibrium.time = temp_2d.equil.t;
+ids_equilibrium_description.time = '.t subfield from: [''equil'']';
%
tens_time = -1;
-vacuum_toroidal_field.b0=gdat(params_equilibrium.shot,'b0','machine',machine);
-vacuum_toroidal_field_desc.b0 = ['''b0'' with interpos using tens_time = ' num2str(tens_time)];
-vacuum_toroidal_field_desc.r0 = '.r0 subfield from: [''b0'']';
-ids_equilibrium.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
-ids_equilibrium.vacuum_toroidal_field.b0 = interpos(vacuum_toroidal_field.b0.t,vacuum_toroidal_field.b0.data, ids_equilibrium.time,tens_time);
+% $$$ vacuum_toroidal_field.b0=gdat(params_equilibrium.shot,'b0','machine',machine);
+% $$$ vacuum_toroidal_field_desc.b0 = ['''b0'' with interpos using tens_time = ' num2str(tens_time)];
+% $$$ vacuum_toroidal_field_desc.r0 = '.r0 subfield from: [''b0'']';
+% $$$ ids_equilibrium.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
+% $$$ ids_equilibrium.vacuum_toroidal_field.b0 = interpos(vacuum_toroidal_field.b0.t,vacuum_toroidal_field.b0.data, ids_equilibrium.time,tens_time);
+% $$$ ids_equilibrium_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;
+vacuum_toroidal_field.b0.data = temp_2d.equil.jpol(end,:).*2e-7 / R0exp;
+vacuum_toroidal_field.b0.t = temp_2d.equil.t;
+vacuum_toroidal_field_desc.b0 = ['''equil.jpol'' *2e-7 divided by R0exp=' num2str(R0exp)];
+vacuum_toroidal_field_desc.r0 = 'set by hand in aug_get_ids_equilibrium';
+ids_equilibrium.vacuum_toroidal_field.r0 = R0exp;
+ids_equilibrium.vacuum_toroidal_field.b0 = vacuum_toroidal_field.b0.data;
ids_equilibrium_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;
ids_equilibrium.time_slice(1:length(ids_equilibrium.time)) = ids_equilibrium.time_slice(1);
@@ -42,7 +53,8 @@ ids_equilibrium.time_slice(1:length(ids_equilibrium.time)) = ids_equilibrium.tim
% $$$ end
% $$$ temp_desc.eqdsks{1} = '''eqdsk'',''time'',ids_equilibrium.time(it)';
-global_quantities.area = gdat(params_equilibrium.shot,'area_edge','machine',machine);
+global_quantities.area.data = temp_2d.equil.area(end,:);
+global_quantities.area.t = temp_2d.equil.t;
global_quantities_desc.area = 'area_edge';
global_quantities.beta_normal = gdat(params_equilibrium.shot,'betan','machine',machine);
global_quantities_desc.beta_normal = 'betan';
@@ -61,8 +73,8 @@ temp.r_magnetic_axis = gdat(params_equilibrium.shot,'r_axis','machine',machine);
temp_desc.r_magnetic_axis = 'r_axis';
temp.z_magnetic_axis = gdat(params_equilibrium.shot,'z_axis','machine',machine);
temp_desc.z_magnetic_axis = 'z_axis';
-temp.psi_axis = gdat(params_equilibrium.shot,'psi_axis','machine',machine); % needs to add psi_edge since psi_axis ??? assuming 0 dege value
-temp_desc.psi_axis = 'psi_axis';
+global_quantities.psi_axis = gdat(params_equilibrium.shot,'psi_axis','machine',machine);
+global_quantities_desc.psi_axis = 'psi_axis';
global_quantities.psi_boundary = gdat(params_equilibrium.shot,'psi_edge','machine',machine);
global_quantities_desc.psi_boundary = 'psi_edge';
global_quantities.q_95 = gdat(params_equilibrium.shot,'q95','machine',machine);
@@ -76,14 +88,14 @@ global_quantities.w_mhd = gdat(params_equilibrium.shot,'w_mhd','machine',machine
global_quantities_desc.w_mhd = 'w_mhd';
global_quantities_fieldnames = fieldnames(global_quantities);
-special_fields = {'magnetic_axis', 'psi_axis', 'q_min'}; % fields needing non-automatic treatments
+special_fields = {'magnetic_axis', 'q_min'}; % fields needing non-automatic treatments
for i=1:length(global_quantities_fieldnames)
try
if ~any(strcmp(global_quantities_fieldnames{i},special_fields))
if ~isstruct(ids_equilibrium.time_slice{1}.global_quantities.(global_quantities_fieldnames{i}))
if length(global_quantities.(global_quantities_fieldnames{i}).data) > length(ids_equilibrium_description.time)
% assume need interpolation with constant extrapolation
- aa = interpos(global_quantities.(global_quantities_fieldnames{i}).t, ...
+ aa = interpos(63,global_quantities.(global_quantities_fieldnames{i}).t, ...
global_quantities.(global_quantities_fieldnames{i}).data,ids_equilibrium.time,tens_time);
for it=1:length(ids_equilibrium.time)
ids_equilibrium.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}) = aa(it);
@@ -109,20 +121,19 @@ for i=1:length(global_quantities_fieldnames)
special_fields{end+1} = global_quantities_fieldnames{i};
end
catch ME_global
- throw(ME_global)
+ rethrow(ME_global)
keyboard
end
end
-
% special case
+aar = interpos(63,temp.r_magnetic_axis.t,temp.r_magnetic_axis.data,ids_equilibrium.time,tens_time);
+aaz = interpos(63,temp.z_magnetic_axis.t,temp.z_magnetic_axis.data,ids_equilibrium.time,tens_time);
for it=1:length(ids_equilibrium.time)
- ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.r = temp.r_magnetic_axis.data(it);
- ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.z = temp.z_magnetic_axis.data(it);
- ids_equilibrium.time_slice{it}.global_quantities.psi_axis = temp.psi_axis.data(it) + ...
- ids_equilibrium.time_slice{it}.global_quantities.psi_boundary;
- [zz,izz]=min(temp.q_rho.data(:,it));
+ ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.r = aar(it);
+ ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.z = aaz(it);
+ [zz,izz]=min(temp_2d.equil.qvalue(:,it));
ids_equilibrium.time_slice{it}.global_quantities.q_min.value = zz;
- ids_equilibrium.time_slice{it}.global_quantities.q_min.rho_tor_norm = temp.q_rho.grids_1d.rhotornorm(izz);
+ ids_equilibrium.time_slice{it}.global_quantities.q_min.rho_tor_norm = temp_2d.equil.rhotornorm(izz,it);
end
% for boundary in addition to lcfs
@@ -162,7 +173,6 @@ temp_desc.zgeom = 'zgeom';
% $$$ temp_desc.z_lcfs = 'z_contour_edge';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-temp_2d.equil = gdat(params_equilibrium.shot,'equil','machine',machine);
% at this stage use geteqdsk to get R,Z boundary since already included
for it=1:length(temp_2d.equil.t)
[eqdskAUG{it}, dummy, equil_t_index]=geteqdskAUG(temp_2d.equil,temp_2d.equil.t(it),[],[],[],[],[],'fignb',999);
@@ -170,15 +180,16 @@ end
boundary_fieldnames = fieldnames(boundary);
special_fields = {'lcfs', 'geometric_axis', 'x_point','outline'}; % fields needing non-automatic treatments
-for it=1:length(ids_equilibrium.time)
- for i=1:length(boundary_fieldnames)
- if ~any(strcmp(boundary_fieldnames{i},special_fields))
- if ~isstruct(ids_equilibrium.time_slice{it}.boundary.(boundary_fieldnames{i}))
- ids_equilibrium.time_slice{it}.boundary.(boundary_fieldnames{i}) = ...
- boundary.(boundary_fieldnames{i}).data(it);
- else
- special_fields{end+1} = boundary_fieldnames{i};
+for i=1:length(boundary_fieldnames)
+ if ~any(strcmp(boundary_fieldnames{i},special_fields))
+ if ~isstruct(ids_equilibrium.time_slice{1}.boundary.(boundary_fieldnames{i}))
+ aa = interpos(63,boundary.(boundary_fieldnames{i}).t,boundary.(boundary_fieldnames{i}).data, ...
+ ids_equilibrium.time,tens_time);
+ for it=1:length(ids_equilibrium.time)
+ ids_equilibrium.time_slice{it}.boundary.(boundary_fieldnames{i}) = aa(it);
end
+ else
+ special_fields{end+1} = boundary_fieldnames{i};
end
end
end
@@ -191,13 +202,20 @@ for it=1:length(ids_equilibrium.time)
ids_equilibrium.time_slice{it}.boundary.lcfs.r = ids_equilibrium.time_slice{it}.boundary.outline.r;
ids_equilibrium.time_slice{it}.boundary.lcfs.z = ids_equilibrium.time_slice{it}.boundary.outline.z;
end
- if ~isempty(temp.rgeom.data)
- ids_equilibrium.time_slice{it}.boundary.geometric_axis.r = temp.rgeom.data(it);
+end
+if ~isempty(temp.rgeom.data)
+ aa = interpos(63,temp.rgeom.t,temp.rgeom.data,ids_equilibrium.time,tens_time);
+ for it=1:length(ids_equilibrium.time)
+ ids_equilibrium.time_slice{it}.boundary.geometric_axis.r = aa(it);
end
- if ~isempty(temp.zgeom.data)
- ids_equilibrium.time_slice{it}.boundary.geometric_axis.z = temp.zgeom.data(it);
+end
+if ~isempty(temp.zgeom.data)
+ aa = interpos(63,temp.zgeom.t,temp.zgeom.data,ids_equilibrium.time,tens_time);
+ for it=1:length(ids_equilibrium.time)
+ ids_equilibrium.time_slice{it}.boundary.geometric_axis.z = aa(it);
end
- % could use temp_2d.equil.Xpoints but not sure which points to take
+end
+% could use temp_2d.equil.Xpoints but not sure which points to take
% $$$ if ~isempty(temp.n_x_point.data) && temp.n_x_point.data(it) > 0
% $$$ ids_equilibrium.time_slice{it}.boundary.x_point(1:temp.n_x_point.data(it)) = ids_equilibrium.time_slice{it}.boundary.x_point(1);
% $$$ for i=1:length(temp.n_x_point.data(it))
@@ -206,8 +224,6 @@ for it=1:length(ids_equilibrium.time)
% $$$ end
% $$$ end
-end
-
%
%% profiles_1d (cannot use eqdsk since not same radial mesh), but can interpolate and use base from equil
%
@@ -229,7 +245,7 @@ dvolume_dpsi = temp_2d.equil.dvoldpsi;
% dvolume_drho_tor = gdat(params_equilibrium.shot,'dvolume_drho_tor','machine',machine);
% elongation = gdat(params_equilibrium.shot,'elongation','machine',machine);
profiles_1d.f_df_dpsi.data = temp_2d.equil.ffprime;
-profiles_1d.f = gdat(params_equilibrium.shot,'rbphi_rho','machine',machine);
+profiles_1d.f.data = temp_2d.equil.jpol*2e-7;
% geometric_axis = gdat(params_equilibrium.shot,'geometric_axis','machine',machine);
% gm1 = gdat(params_equilibrium.shot,'gm1','machine',machine);
% gm2 = gdat(params_equilibrium.shot,'gm2','machine',machine);
@@ -244,14 +260,16 @@ profiles_1d.f = gdat(params_equilibrium.shot,'rbphi_rho','machine',machine);
% j_tor = gdat(params_equilibrium.shot,'j_tor','machine',machine);
% magnetic_shear = gdat(params_equilibrium.shot,'magnetic_shear','machine',machine);
% mass_density = gdat(params_equilibrium.shot,'mass_density','machine',machine);
-profiles_1d.phi = gdat(params_equilibrium.shot,'phi_tor','machine',machine);
-profiles_1d.pressure = gdat(params_equilibrium.shot,'pressure','machine',machine);
+profiles_1d.phi.data = temp_2d.equil.pressure;
+profiles_1d.pressure.data = temp_2d.equil.pressure;
% psi = gdat(params_equilibrium.shot,'psi_rho','machine',machine); % (could take from .x of any like rhotor and psi_axis, psi_edge from global_quantities)
profiles_1d.psi.data = temp_2d.equil.psi;
-profiles_1d.q = gdat(params_equilibrium.shot,'q_rho','machine',machine);
-profiles_1d.rho_tor = gdat(params_equilibrium.shot,'rhotor','machine',machine);
+profiles_1d.q.data = temp_2d.equil.qvalue;
+for it=1:length(temp_2d.equil.t)
+ profiles_1d.rho_tor.data(:,it) = sqrt(temp_2d.equil.phi(:,it)./vacuum_toroidal_field.b0.data(it)/pi);
+end
%rho_tor_norm = gdat(params_equilibrium.shot,'rhotor_norm','machine',machine); % from rho_tor
-profiles_1d.rho_volume_norm = gdat(params_equilibrium.shot,'rhovol','machine',machine);
+profiles_1d.rho_volume_norm.data = temp_2d.equil.rhovolnorm;
% r_inboard = gdat(params_equilibrium.shot,'r_inboard','machine',machine);
% r_outboard = gdat(params_equilibrium.shot,'r_outboard','machine',machine);
% squareness_lower_inner = gdat(params_equilibrium.shot,'squareness_lower_inner','machine',machine);
@@ -262,7 +280,7 @@ profiles_1d.rho_volume_norm = gdat(params_equilibrium.shot,'rhovol','machine',ma
% trapped_fraction = gdat(params_equilibrium.shot,'trapped_fraction','machine',machine);
% triangularity_lower = gdat(params_equilibrium.shot,'triangularity_lower','machine',machine);
% triangularity_upper = gdat(params_equilibrium.shot,'triangularity_upper','machine',machine);
-profiles_1d.volume = gdat(params_equilibrium.shot,'volume_rho','machine',machine);
+profiles_1d.volume.data = temp_2d.equil.vol;
profiles_1d_fieldnames = fieldnames(profiles_1d);
special_fields = {'geometric_axis', 'rho_tor_norm', 'psi'}; % fields needing non-automatic treatments
@@ -283,15 +301,12 @@ for it=1:length(ids_equilibrium.time)
end
end
% special cases
-nrho = length(profiles_1d.rho_tor.x);
-ntime = length(temp.psi_axis.data);
for it=1:length(ids_equilibrium.time)
ids_equilibrium.time_slice{it}.profiles_1d.rho_tor_norm = ids_equilibrium.time_slice{it}.profiles_1d.rho_tor./ ...
ids_equilibrium.time_slice{it}.profiles_1d.rho_tor(end);
% ids_equilibrium.time_slice{it}.profiles_1d.psi = temp.psi_axis.data(it) + ...
% profiles_1d.rho_tor.x(:,it).^2.*(global_quantities.psi_boundary.data(it)-temp.psi_axis.data(it));
end
-
%
%% profiles_2d{1} ala eqdsk, only this one thus grid_type=1
%
--
GitLab