%% RAPTOR tutorial - H mode

close all hidden;
clear; run('../RAPTOR_path.m');

verbosity = 1;

[config] = RAPTOR_config('TCV'); % load default config
config.grid.tgrid = [0:2e-3:0.2];
config.grid.rhogrid = linspace(0,1,21); % spatial grid

% set ECH/ECCD
config.echcd.Configuration.n_units = 2;

% config.ne.method = 'state';
% config.ne.BC.rhoSource = 'v';
% config.ni.method = 'nescal';

config.chi_e = RAPTORmodule('chi_FF');

% rerun the config file with the new parameter set
[model,params,init,g,v,U] = build_RAPTOR_model(config);

% set ECH/ECCD parameters
params.echcd.active = true;
params.echcd.rdep = [0.4,0.0]; %
params.echcd.wdep = [0.35,0.35]; % deposition width
params.echcd.cd_eff = [1,0]; % CD efficiency

% plasma current
rampfun = @(t,tmin,ymin,tmax,ymax) max(ymin,min((ymax-ymin)/(tmax-tmin)*(t-tmin),ymax-ymin)+ymin); % anonymous function for ramps

U(1,:)= 200e3*ones(size(params.tgrid)); % input Ip trace: ramp from 80kA to 200kA

% EC actuator power
U(2,:) = rampfun(params.tgrid,0,0.05, 0.1,0.5e6);
U(3,:) = rampfun(params.tgrid,0,0.1 , 0.15,0.5e6);

% initial condition
init.Ip0 = U(1,1);
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); 

% tweak transport
params.chi_e.cano = 10;

% run RAPTOR
[simres_noenhance] = RAPTOR_predictive(x0,g,v,U,model,params,'verbosity',verbosity);
out = RAPTOR_out(simres_noenhance,model,params);

%% Activate H-mode module
% Activate H-mode
Hmode_activate = (params.tgrid>0.08 & params.tgrid<0.18);
plot(params.tgrid,Hmode_activate);

v = v*ones(size(params.tgrid)); % make time-varying V
v(model.hmode.vind.activation,:) = Hmode_activate; % activation flag
v(model.hmode.vind.rhoped,:) = 0.9; % set desired pedestal rho value
v(model.hmode.vind.gteped,:) = -5e3; % set desired dTe/drho value
params.hmode.active = true;

% re-run simulation
[simres_Hmode] = RAPTOR_predictive(x0,g,v,U,model,params,'verbosity',verbosity);
out_Hmode = RAPTOR_out(simres_Hmode,model,params);
RAPTOR_plot_GUI({out,out_Hmode},'overview')

%% more details and debugging options:
params.hmode.doplot = true;
[simres_Hmode] = RAPTOR_predictive(x0,g,v,U,model,params,'verbosity',verbosity);

%% H-mode module: impose a linear pedestal for Te
% During H-mode, the computational domain for Te is reduced to rho=[0,rhoped].
% This is done by compressing the original spline basis in the reduced
% domain, imposing the BC for the Te equation at rho=rhoped. As the
% alternative spline matrices are pre-calculated, rhoped is not allowed to
% vary during the H-mode.
% The psi equation is still solved on the full radial grid rho=[0,1].
% Optionally, one can also solve for Ti and/or ne. These state variables
% are also solved for on the reduced domain, the user imposes a linear
% pedestal.

[config] = RAPTOR_config; 

t_LH = .01;
t_HL = .04;

hmode_activate = (config.grid.tgrid>t_LH & config.grid.tgrid<t_HL);
% flag to identify L-mode and H-mode time windows

te_BC_lmode = 200; 
ti_BC_lmode = 200;
ne_BC_lmode = 1e18;
% boundary condition during L-mode window, both rho_BC and Te(rho_BC) are allowed to vary

config.grid.rhogrid = linspace(0,1,31);

config.ne.method = 'state';
config.ti.method = 'state';

config.hmode = hmode('imposed');
config.hmode.params.active = true; 
config.hmode.params.te_rhoped = 1e3;   
config.hmode.params.te_rhoedge = 200;
config.hmode.params.ti_rhoped = 0.8e3;   
config.hmode.params.ti_rhoedge = 200;
config.hmode.params.ne_rhoped = 5e18;   
config.hmode.params.ne_rhoedge = 1e18;
config.hmode.params.rhoped = .90;
% These values impose the linear pedestal:
% * Te(rho=rhoped) = Te_rhoped
% * Te(rho=1) = Te_rhoedge
% (same for Ti, ne)

config.echcd = RAPTORmodule('echcd_gaussian');
config.echcd.Configuration.n_units = 1;

config.nbhcd = RAPTORmodule('hcd_gaussian');
config.nbhcd.Configuration.n_units = 1;

[model,params,init,g,v,U] = build_RAPTOR_model(config);

params.pei.active = true; % equipartition heating

params.echcd.rdep = 0;
params.echcd.wdep = 0.3;


params.nbhcd.rdep = 0.2;
params.nbhcd.frac_to_electrons = 0.2;
params.nbhcd.cd_eff = 0;

init.te0 = 1e3;
init.ti0 = .1*1e3;
init.ni0 = 1e19;
init.Ip0 = 250e3;

% extend the v vector to a V matrix (allowing to vary v in time)
V = repmat(v,1,numel(params.tgrid));
V(model.hmode.vind.activation, ~hmode_activate) = 0;
% H-mode time window
V(model.te.BC.vind_value, ~hmode_activate) = te_BC_lmode;
V(model.ti.BC.vind_value, ~hmode_activate) = ti_BC_lmode;
V(model.ne.BC.vind_value, ~hmode_activate) = ne_BC_lmode;
% Te BC during L-mode 
% By default, Te(rho=rhoped) is copied in V(model.te.BC.vind_value, :),
% will be overwritten for L-mode time window (same for Ti, ne)

U(1,:) = init.Ip0;
U(2,:) = .5e6;
U(3,:) = 1e6;

x0 = RAPTOR_initial_conditions(model,init,g(:,1),V(:,1)); 

simres = RAPTOR_predictive(x0,g,V,U,model,params,'verbosity',2);

out = RAPTOR_out(simres,model,params);

%%
tic
% plot of the time evolution of Te(rho)
for i = 1:numel(params.tgrid)
    figure(1)
    subplot(3,1,1); plot(out.rho, out.te(:,i))
    title('Te');
    hold on
    subplot(3,1,2); plot(out.rho, out.ti(:,i))
    title('Ti')
    hold on
    subplot(3,1,3); plot(out.rho, out.ne(:,i))
    title('ne')
    hold on
    pause(.1)
end
toc


