%% RAPTOR tutorial -- Troubleshooting and debugging
% Things do not always go perfectly, so debugging options are useful
% When things go wrong numerically, it is usually because of some
% 'unreasonable' physical setting. For example a too rapid change of plasma
% current, violent addition of too much power at the same time, or other sudden jumps.
% In this case, try changing the grid sizes (time and rho) or go back to
% settings that you know worked, and evolve slowly from there. The tools
% presented here are for advanced users who want to check things in
% detail.

% First, configure a standard run with some central NBI power
close all hidden;
clear; run ../RAPTOR_path.m;

[config] = RAPTOR_config; % load default params
config.grid.tgrid = 0:0.001:0.002; % choose a smaller time grid just for demo
config.grid.rhogrid = linspace(0,1,21); % spatial grid
config.debug.iterplot = 0; %10; % plot profiles every N iterations
config.debug.iterdisp = 20; % display progress every N iterations
config.numerics.nmax = uint32(1); % only 1 newton iterations to avoid overcrowding this tutorial with results
% heating & current drive parameters
warning('off','RAPTOR:ExceedMaxNewtonIterations')

config.echcd = RAPTORmodule('none'); % switch off echcd

config.nbhcd.Configuration.n_units = 1;

% create model, params, init, geometry g, kinetic profiles v
[model,params,init,g,v,U] = build_RAPTOR_model(config);

% set NBI parameters
params.nbhcd.active = true;
params.nbhcd.rdep = [0]; % -1 to indicate commanded extrernally
params.nbhcd.wdep = [0.8]; % broad heating
params.nbhcd.cd_eff = [0]; % no current drive

% plasma current
U(1,:)= 200e3*ones(size(params.tgrid)); % input Ip trace: ramp from 40kA to 200kA
U(1,params.tgrid<0.02) = linspace(80e3,200e3,sum(params.tgrid<0.02));

% NBI
U(2,:) = zeros(size(params.tgrid)); U(2,params.tgrid>0.04) = 2e6;

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

%% Verbosity
% Set verbosity parameter to get more or less information during run
disp(params)
simres = RAPTOR_predictive(x0,g,v,U,model,params);
% To use user-defined iterdisp options, set verbosity to -1, details in
% "set_verbosity.m"
params.debug.iterdisp=1;
simres = RAPTOR_predictive(x0,g,v,U,model,params,'verbosity',-1);
%% Numerical check of state equation jacobians
params.debug.checkstateequation = true;
params.debug.checkprofilejacobians = false;
simres = RAPTOR_predictive(x0,g,v,U,model,params);

%% Convergence check of a single module derivative
params.debug.checkstateequation = false;
params.debug.checkprofilejacobians = true;
params.debug.trapProfilesToCheck = 'jbsB'; 
params.debug.stapProfilesToCheck = 'none'; 

simres = RAPTOR_predictive(x0,g,v,U,model,params);

%% Other options
% Many other debugging options are listed below but are not executed here.
% file but not published as html. Please inspect the .m to learn more.
return
%% Detailed convergence check of all module derivatives
params.debug.checkstateequation = false;
params.debug.checkprofilejacobians = true;
simres = RAPTOR_predictive(x0,g,v,U,model,params);

%% "Keyboard" at pre-set time.
% You may want to pause the simulation at a given time step, to observe in
% detail what is happening. For example, we will pause for 10 iterations after the 
% at the time of NBI switch on and advance slowly (by 'continuing' the simulation each time)

params.debug.keyboard_at = find(U(2,:)>0,10,'first');
% We do not run the simulation for the sake of tutorial "publishing" but 
% you can try yourself by uncommenting.
%[simres] = RAPTOR_predictive(x0,g0,v0,U,model,params); % Uncomment this

%% Stop at each Newton iteration, plot and display salient features
params.debug.keyboard_at = []; % remove this check
params.debug.plotEachNewtonStep = 1;
% [simres] = RAPTOR_predictive(x0,g0,v0,U,model,params); We do not run it
% for the tutorial "publishing" but you can try yourself by uncommenting

%% Numerical checking of Jacobians of global solver.
% This advanced option is to check that the state equation jacobians are
% computed properly, to check if any bugs have been introduced by new additions to the code.
params.debug.plotEachNewtonStep = false; % remove previous check
params.debug.checkstateequation = true;
%[simres] = RAPTOR_predictive(x0,g0,v0,U,model,params); % We do not run it
% for the tutorial "publishing" but you can try yourself by uncommenting.

%% Plotting Jacobian and state equation residual structure
params.debug.plotJacobians = true; % remove this check
params.debug.checkstateequation = false;
%[simres] = RAPTOR_predictive(x0,g0,v0,U,model,params); % We do not run it
% for the tutorial "publishing" but you can try yourself by uncommenting.

%% Stopping in case of a given warning:
dbstop if warning RAPTOR:ExceedMaxNewtonIterations