From 193c7fe6fc63b3bbaffe83bf261f1b8062bfa0ed Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Mon, 13 May 2019 08:21:28 +0000
Subject: [PATCH] finished implementing and testing
torbeam_prepare_inputs_and_run and related functions, will test reading
TORBEAM angles from shotfile from Giovanni
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@11874 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/AUG/torbeam_prepare_inputs_and_run.m | 71 +++++++++++---------
1 file changed, 41 insertions(+), 30 deletions(-)
diff --git a/crpptbx/AUG/torbeam_prepare_inputs_and_run.m b/crpptbx/AUG/torbeam_prepare_inputs_and_run.m
index 9925badf..34cad038 100644
--- a/crpptbx/AUG/torbeam_prepare_inputs_and_run.m
+++ b/crpptbx/AUG/torbeam_prepare_inputs_and_run.m
@@ -3,6 +3,11 @@ function [torbeam_out,torbeam_in,varargout] = torbeam_prepare_inputs_and_run(sho
% [torbeam_out,torbeam_in,varargout] = aaa(shot,times_in,varargin);
%
% varargin{1}: torbeam_in (from a previous run, so avoids to reload all the data)
+% leave then times_in empty to rerun at all times in torbeam_in.inbeam.t_some_beam_on or select times points from that array
+%
+% Examples:
+% [torbeam_out,torbeam_in] = torbeam_prepare_inputs_and_run(30594,linpace(0,10,101));
+% [torbeam_out2] = torbeam_prepare_inputs_and_run(30594,torbeam_in.inbeam.t_some_beam_on(1:30:end),torbeam_in); % to rerun a few cases
%
if nargin==0 || isempty(shot)
@@ -32,8 +37,8 @@ else
end
% base all times on Te data (most important? although ne as well), note take Thomson by default to avoid cut-off
-if ~exist('times_for_torbeam_eff') || isempty(times_for_torbeam_eff),
- times_for_torbeam_eff = torbeam_in.nete.fit.t;
+if doload && (~exist('times_for_torbeam_eff') || isempty(times_for_torbeam_eff)),
+ times_for_torbeam_eff = torbeam_in.nete.fit.t(torbeam_in.nete.fit.t>0);
end
% load beam releated data including theta,phi in torbeam reference
if doload,
@@ -47,12 +52,16 @@ end
% if torbeam_in provided as input (with its own time input sequence) and times_in provided as well: need to adapt times_in to match available inputs
% (otherwise would take closest eqdsk)
if doload==0
- times_for_torbeam_eff = intersect(torbeam_in.inbeam.t_some_beam_on,times_for_torbeam_eff);
-else
- times_for_torbeam_eff = torbeam_in.inbeam.t_some_beam_on;
+ if ~isempty(times_for_torbeam_eff)
+ times_for_torbeam_eff = intersect(torbeam_in.inbeam.t_some_beam_on,times_for_torbeam_eff);
+ else
+ times_for_torbeam_eff = torbeam_in.inbeam.t_some_beam_on;
+ end
end
if doload,
torbeam_in.eqd_all=gdat(shot,'eqdsk','time',torbeam_in.inbeam.t_some_beam_on,'write',0);
+ torbeam_in.dir_torbeam = ['/tmp/' getenv('USER') '/' num2str(shot) '_torbeam'];
+ unix(['mkdir ' torbeam_in.dir_torbeam]);
else
if ~isfield(torbeam_in,'eqd_all') || isempty(torbeam_in.eqd_all)
error('expects eqd_all as field of torbeam_in (varargin{1}) from gdat(shot,''eqdsk'',''time'',torbeam_in.inbeam.t_some_beam_on,''write'',0);')
@@ -65,8 +74,6 @@ if doload,
nb_points_prof = 46;
rhofit = linspace(0,1,nb_points_prof);
- dir_torbeam = ['/tmp/' getenv('USER') '/' num2str(shot) '_torbeam'];
- unix(['mkdir ' dir_torbeam]);
tb_par_fields=fieldnames(torbeam_in.inbeam.tb_par);
itime_te = iround_os(torbeam_in.nete.fit.t,times_for_torbeam_eff);
@@ -77,7 +84,7 @@ if doload,
for it=1:length(torbeam_in.inbeam.t_some_beam_on)
% Te
it_te = itime_te(it);
- torbeam_in.fname{it}.te = fullfile(dir_torbeam,['Te' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it)) '.dat']);
+ torbeam_in.fname{it}.te = fullfile(torbeam_in.dir_torbeam,['Te' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it)) '.dat']);
fid=fopen(torbeam_in.fname{it}.te,'w'); % te.dat file etc per time
fprintf(fid,'%d\n',nb_points_prof);
% tefit = interpos(torbeam_in.nete.fit.rhopolnorm(:,it_te) , torbeam_in.nete.fit.te.data(:,it_te),rhofit,-1,[1 2],[0 0]);
@@ -87,7 +94,7 @@ if doload,
fclose(fid);
% ne
it_ne = itime_ne(it);
- torbeam_in.fname{it}.ne = fullfile(dir_torbeam,['ne' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it)) '.dat']);
+ torbeam_in.fname{it}.ne = fullfile(torbeam_in.dir_torbeam,['ne' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it)) '.dat']);
fid=fopen(torbeam_in.fname{it}.ne,'w'); % ne.dat file etc per time
fprintf(fid,'%d\n',nb_points_prof);
nefit = interpos(torbeam_in.nete.fit.rhopolnorm(:,it_ne) , torbeam_in.nete.fit.ne.data(:,it_ne),rhofit,-0.1,[1 0],[0 0]);
@@ -99,7 +106,7 @@ if doload,
eqd_with_b2=eqdsk_cocos_transform(torbeam_in.eqd_all.eqdsk{it},[17 cocos_out]);
eqd_with_b1=eqdsk_transform(eqd_with_b2,65,129,-0.1);
eqd_with_b = read_eqdsk(eqd_with_b1,cocos_out,1); % to compute BR, etc
- torbeam_in.fname{it}.topfile = fullfile(dir_torbeam,['topfile' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it))]);
+ torbeam_in.fname{it}.topfile = fullfile(torbeam_in.dir_torbeam,['topfile' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it))]);
fid=fopen(torbeam_in.fname{it}.topfile,'w');
fprintf(fid,'%s %d %.6f\n',' Number of radial and vertical grid points in ',shot,torbeam_in.inbeam.t_some_beam_on(it));
fprintf(fid,'%d %d\n',length(eqd_with_b.rmesh),length(eqd_with_b.zmesh));
@@ -133,7 +140,7 @@ if doload,
% extract beam info at given time
for ibeam=find(torbeam_in.inbeam.beam_on(:,it_inb(it)))',
- torbeam_in.fname{it}.inbeam{ibeam} = fullfile(dir_torbeam,['inbeam' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it)) 'g' num2str(ibeam) '.dat']);
+ torbeam_in.fname{it}.inbeam{ibeam} = fullfile(torbeam_in.dir_torbeam,['inbeam' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it)) 'g' num2str(ibeam) '.dat']);
[aa,bb]=unix(['cp ~/TORBEAM/inbeam_default.dat ' torbeam_in.fname{it}.inbeam{ibeam}]);
fid=fopen(torbeam_in.fname{it}.inbeam{ibeam},'a');
fprintf(fid,' xf = %f, ! wave frequency om=2*pi*xf\n',torbeam_in.inbeam.freq(ibeam));
@@ -154,48 +161,52 @@ end
% run TORBEAM
current_dir = pwd;
-eval(['cd ' dir_torbeam])
+eval(['cd ' torbeam_in.dir_torbeam])
% namelist, ol proj 3 rays(R,z)*100, not used now, tor proj 3 rays(R,Y), vol() dP/dV(r/a) j(r/a)
%tb_out_files_to_cp = {'echo.dat', 't1_LIB.dat', 't2_LIB.dat', 't1tor_LIB.dat', 'volumes.dat', 't2_new_LIB.dat', 'cntpr.dat'};
tb_out_files_to_cp = {'t1_LIB.dat', 't1tor_LIB.dat', 'volumes.dat', 't2_new_LIB.dat'};
tb_out_files_to_cp_varname = strrep(tb_out_files_to_cp,'.dat','');
torbeam_out.run_torbeam = false .* ones(size(torbeam_in.inbeam.power_t,1),length(times_for_torbeam_eff));
it_inb_beam_on = iround_os(torbeam_in.inbeam.t_some_beam_on,times_for_torbeam_eff);
+it_inb = iround_os(torbeam_in.inbeam.t,times_for_torbeam_eff);
torbeam_out.t = torbeam_in.inbeam.t_some_beam_on(it_inb_beam_on);
tic
+icount = 0;
for it_eff=1:length(times_for_torbeam_eff)
it = it_inb_beam_on(it_eff);
unix(['cp ' torbeam_in.fname{it}.te ' Te.dat']);
unix(['cp ' torbeam_in.fname{it}.ne ' ne.dat']);
unix(['cp ' torbeam_in.fname{it}.topfile ' topfile']);
- disp([' t=' num2str(torbeam_in.inbeam.t(it_inb(it))) ' beam=' num2str(find(torbeam_in.inbeam.beam_on(:,it_inb(it)))','%d')]);
- for ibeam=find(torbeam_in.inbeam.beam_on(:,it_inb(it)))',
+ disp([' t=' num2str(torbeam_in.inbeam.t(it_inb(it_eff))) ' beam=' num2str(find(torbeam_in.inbeam.beam_on(:,it_inb(it_eff)))','%d')]);
+ for ibeam=find(torbeam_in.inbeam.beam_on(:,it_inb(it_eff)))',
unix(['cp ' torbeam_in.fname{it}.inbeam{ibeam} ' inbeam.dat']);
[a,b]=unix('/afs/ipp-garching.mpg.de/home/o/osauter/TORBEAM/lib-OUT/a.out');
+ icount = icount + 1;
imw = findstr('(MW)',b);
if ~isempty(imw)
- torbeam_out.run_torbeam(ibeam,it) = true;
- torbeam_out.pabs(ibeam,it) = 1e6 * sscanf(b(imw+4:end),'%f',1);
+ torbeam_out.run_torbeam(ibeam,it_eff) = true;
+ torbeam_out.pabs(ibeam,it_eff) = 1e6 * sscanf(b(imw+4:end),'%f',1);
ima = findstr('(MA)',b);
- torbeam_out.icd(ibeam,it) = 1e6 * sscanf(b(ima+4:end),'%f',1);
+ torbeam_out.icd(ibeam,it_eff) = 1e6 * sscanf(b(ima+4:end),'%f',1);
imaxat = findstr('max. at',b);
- torbeam_out.pabs_roamax(ibeam,it) = sscanf(b(imaxat(1)+7:end),'%f',1);
- torbeam_out.icd_roamax(ibeam,it) = sscanf(b(imaxat(2)+7:end),'%f',1);
+ torbeam_out.pabs_roamax(ibeam,it_eff) = sscanf(b(imaxat(1)+7:end),'%f',1);
+ torbeam_out.icd_roamax(ibeam,it_eff) = sscanf(b(imaxat(2)+7:end),'%f',1);
for j=1:length(tb_out_files_to_cp)
- torbeam_out.(tb_out_files_to_cp_varname{j}){ibeam,it} = load(fullfile(dir_torbeam,[tb_out_files_to_cp{j}]));
+ torbeam_out.(tb_out_files_to_cp_varname{j}){ibeam,it_eff} = load(fullfile(torbeam_in.dir_torbeam,[tb_out_files_to_cp{j}]));
end
nnn = load('cntpr.dat');
- torbeam_out.cntpr.tbegabs(ibeam,it) = nnn(1);
- torbeam_out.cntpr.dtabs(ibeam,it) = nnn(2);
- torbeam_out.pdens_jcd_rho(:,ibeam,it) = torbeam_out.t2_new_LIB{ibeam,it}(:,1);
- torbeam_out.pdens(:,ibeam,it) = torbeam_out.t2_new_LIB{ibeam,it}(:,2);
- torbeam_out.jcd(:,ibeam,it) = torbeam_out.t2_new_LIB{ibeam,it}(:,3);
+ torbeam_out.cntpr.tbegabs(ibeam,it_eff) = nnn(1);
+ torbeam_out.cntpr.dtabs(ibeam,it_eff) = nnn(2);
+ torbeam_out.pdens_jcd_rho(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,1);
+ torbeam_out.pdens(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,2);
+ torbeam_out.jcd(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,3);
else
disp(['TORBEAM did not run properly, no MW result'])
end
end
end
-disp(['time for running all TORBEAMs: ' num2str(toc)])
+aaa=toc;
+disp(['time for running all TORBEAMs: ' num2str(aaa) ' for ' num2str(icount) ' runs so ' num2str(aaa/icount) 's per TORBEAM run'])
for ibeam=find(any(torbeam_out.run_torbeam,2))'
figure;
subplot(2,1,1)
@@ -211,16 +222,16 @@ end
figure
subplot(3,1,1)
-plot(torbeam_out.t,torbeam_out.pabs/1e6)
+plot(torbeam_out.t,torbeam_out.pabs/1e6,'*-')
ylabel('Pabs [MW]')
subplot(3,1,2)
-plot(torbeam_out.t,torbeam_out.icd/1e6)
+plot(torbeam_out.t,torbeam_out.icd/1e6,'*-')
ylabel('Icd [MA]')
subplot(3,1,3)
-plot(torbeam_out.t,torbeam_out.pabs_roamax/1e6)
-ylabel('r/a (max Pdens)')
+plot(torbeam_out.t,torbeam_out.pabs_roamax,'*-')
+ylabel('r(max Pdens) / a')
xlabel('t [s]')
--
GitLab