## Contents

- RAPTOR tutorial -- Troubleshooting and debugging
- Numerical check of state equation jacobians
- Convergence check of a single module derivative
- Other options
- Detailed convergence check of all module derivatives
- "Keyboard" at pre-set time.
- Stop at each Newton iteration, plot and display salient features
- Numerical checking of Jacobians of global solver.
- Plotting Jacobian and state equation residual structure

## 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 = 1; % only 1 newton iterations to avoid overcrowding this tutorial with results
    % heating & current drive parameters

    config.echcd = echcd('none');

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

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

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

    % 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));

## Numerical check of state equation jacobians

    params.debug.checkstateequation = true;
    params.debug.checkprofilejacobians = false;
    simres = RAPTOR_predictive(x0,g,v,U,model,params);

    psi te
    3.947e-09  8.626e-08
    3.346e-07  9.159e-09  Warning: Exceeded maximum newton iterations
    Warning: exceeded maximum newton iterations
         it  telaps newt     res   t[ms]  dt[ms]  Ip[kA] Icd[kA] Ibs[kA] Ioh[kA]      qe    qmin      q0   Vl[V] Te0[keV] Ti0[keV] ne0[e19]   f_ss
          1    0.24    1 3.2e-01       0       1      80       0    1.41    78.6    15.4    5.41    5.41 4.1e+00    0.22    0.22    1.00 1.8e+01

## Convergence check of a single module derivative

    params.debug.checkstateequation = false;
    params.debug.checkprofilejacobians = true;
    params.debug.profilestocheck = 'jbsB';
    simres = RAPTOR_predictive(x0,g,v,U,model,params);

     **** Jacobian residuals ******

       dpsi/x(psi) = 4.64e-09      dpsi/x(te) = 0
    dpsi/xdot(psi) = 0          dpsi/xdot(te) = 0
      dpsip/x(psi) = 6.2e-09      dpsip/x(te) = 0
    dpsip/xdot(psi) = 0         dpsip/xdot(te) = 0
     dpsipp/x(psi) = 5.56e-09    dpsipp/x(te) = 0
    dpsipp/xdot(psi) = 0         dpsipp/xdot(te) = 0
       dupl/x(psi) = 0             dupl/x(te) = 0
    dupl/xdot(psi) = 9.28e-11   dupl/xdot(te) = 0
        dne/x(psi) = 0              dne/x(te) = 0
     dne/xdot(psi) = 0           dne/xdot(te) = 0
     dnedot/x(psi) = 0           dnedot/x(te) = 0
    dnedot/xdot(psi) = 0         dnedot/xdot(te) = 0
        dni/x(psi) = 0              dni/x(te) = 0
     dni/xdot(psi) = 0           dni/xdot(te) = 0
     dnidot/x(psi) = 0           dnidot/x(te) = 0
    dnidot/xdot(psi) = 0         dnidot/xdot(te) = 0
       dnep/x(psi) = 0             dnep/x(te) = 0
    dnep/xdot(psi) = 0          dnep/xdot(te) = 0
       dnip/x(psi) = 0             dnip/x(te) = 0
    dnip/xdot(psi) = 0          dnip/xdot(te) = 0
        dn1/x(psi) = 0              dn1/x(te) = 0
     dn1/xdot(psi) = 0           dn1/xdot(te) = 0
        dn2/x(psi) = 0              dn2/x(te) = 0
     dn2/xdot(psi) = 0           dn2/xdot(te) = 0
        dn3/x(psi) = 0              dn3/x(te) = 0
     dn3/xdot(psi) = 0           dn3/xdot(te) = 0
        dte/x(psi) = 0              dte/x(te) = 4.91e-09
     dte/xdot(psi) = 0           dte/xdot(te) = 0
        dti/x(psi) = 0              dti/x(te) = 4.91e-09
     dti/xdot(psi) = 0           dti/xdot(te) = 0
        dze/x(psi) = 0              dze/x(te) = 0
     dze/xdot(psi) = 0           dze/xdot(te) = 0
       dtep/x(psi) = 0             dtep/x(te) = 6.1e-09
    dtep/xdot(psi) = 0          dtep/xdot(te) = 0
       dtip/x(psi) = 0             dtip/x(te) = 6.1e-09
    dtip/xdot(psi) = 0          dtip/xdot(te) = 0
     dtedot/x(psi) = 0           dtedot/x(te) = 0
    dtedot/xdot(psi) = 0         dtedot/xdot(te) = 0
     dtidot/x(psi) = 0           dtidot/x(te) = 0
    dtidot/xdot(psi) = 0         dtidot/xdot(te) = 9.7e-17
         dp/x(psi) = 0               dp/x(te) = 0
      dp/xdot(psi) = 0            dp/xdot(te) = 0
         dq/x(psi) = 1.08e-05        dq/x(te) = 0
      dq/xdot(psi) = 0            dq/xdot(te) = 0
      diota/x(psi) = 6.71e-09     diota/x(te) = 0
    diota/xdot(psi) = 0         diota/xdot(te) = 0
     dshear/x(psi) = 9.66e-07    dshear/x(te) = 0
    dshear/xdot(psi) = 0         dshear/xdot(te) = 0
      dlamH/x(psi) = 0            dlamH/x(te) = 0
    dlamH/xdot(psi) = 0         dlamH/xdot(te) = 0
     dwntms/x(psi) = 0           dwntms/x(te) = 0
    dwntms/xdot(psi) = 0         dwntms/xdot(te) = 0
    dwdotntms/x(psi) = 0         dwdotntms/x(te) = 0
    dwdotntms/xdot(psi) = 0         dwdotntms/xdot(te) = 0
      dchie/x(psi) = 1.86e-06     dchie/x(te) = 0
    dchie/xdot(psi) = 0         dchie/xdot(te) = 0
      dchie/u(psi) = 0            dchie/u(te) = 0
      dchii/x(psi) = 1.86e-06     dchii/x(te) = 0
    dchii/xdot(psi) = 0         dchii/xdot(te) = 0
      dchii/u(psi) = 0            dchii/u(te) = 0
       dvpe/x(psi) = 0             dvpe/x(te) = 0
    dvpe/xdot(psi) = 0          dvpe/xdot(te) = 0
       dvpe/u(psi) = 0             dvpe/u(te) = 0
       ddne/x(psi) = 0             ddne/x(te) = 0
    ddne/xdot(psi) = 0          ddne/xdot(te) = 0
       ddne/u(psi) = 0             ddne/u(te) = 0
       dsne/x(psi) = 0             dsne/x(te) = 0
    dsne/xdot(psi) = 0          dsne/xdot(te) = 0
       dsne/u(psi) = 0             dsne/u(te) = 0
      djbsB/x(psi) = 4.58e-07     djbsB/x(te) = 2.13e-08
    djbsB/xdot(psi) = 0         djbsB/xdot(te) = 0
      djbsB/u(psi) = 0            djbsB/u(te) = 0
    dsigneo/x(psi) = 3.09e-07   dsigneo/x(te) = 9.64e-09
    dsigneo/xdot(psi) = 0         dsigneo/xdot(te) = 0
    dsigneo/u(psi) = 0          dsigneo/u(te) = 0
    djphioR/x(psi) = 3.43e-09   djphioR/x(te) = 0
    djphioR/xdot(psi) = 0         djphioR/xdot(te) = 0
    djphioR/u(psi) = 0          djphioR/u(te) = 0
     dpetot/x(psi) = 3.43e-09    dpetot/x(te) = 0
    dpetot/xdot(psi) = 8.74e-11  dpetot/xdot(te) = 0
     dpetot/u(psi) = 1.99e-09    dpetot/u(te) = 1.55e-09
     dpitot/x(psi) = 0           dpitot/x(te) = 0
    dpitot/xdot(psi) = 0         dpitot/xdot(te) = 0
     dpitot/u(psi) = 0           dpitot/u(te) = 0
     djauxB/x(psi) = 0           djauxB/x(te) = 0
    djauxB/xdot(psi) = 0         djauxB/xdot(te) = 0
     djauxB/u(psi) = 0           djauxB/u(te) = 0
    Warning: Exceeded maximum newton iterations
    Warning: exceeded maximum newton iterations
         it  telaps newt     res   t[ms]  dt[ms]  Ip[kA] Icd[kA] Ibs[kA] Ioh[kA]      qe    qmin      q0   Vl[V] Te0[keV] Ti0[keV] ne0[e19]   f_ss
          1   0.073    1 3.2e-01       0       1      80       0    1.41    78.6    15.4    5.41    5.41 4.1e+00    0.22    0.22    1.00 1.8e+01

## 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.

## Detailed convergence check of all module derivatives

    params.debug.checkstateequation = false;
    params.debug.checkprofilejacobians = true;
    simres = RAPTOR_predictive(x0,g,v,U,model,params);

     **** Jacobian residuals ******

       dpsi/x(psi) = 4.17e-09      dpsi/x(te) = 0
    dpsi/xdot(psi) = 0          dpsi/xdot(te) = 0
      dpsip/x(psi) = 6.31e-09     dpsip/x(te) = 0
    dpsip/xdot(psi) = 0         dpsip/xdot(te) = 0
     dpsipp/x(psi) = 8.93e-09    dpsipp/x(te) = 0
    dpsipp/xdot(psi) = 0         dpsipp/xdot(te) = 0
       dupl/x(psi) = 0             dupl/x(te) = 0
    dupl/xdot(psi) = 9.11e-11   dupl/xdot(te) = 0
        dne/x(psi) = 0              dne/x(te) = 0
     dne/xdot(psi) = 0           dne/xdot(te) = 0
     dnedot/x(psi) = 0           dnedot/x(te) = 0
    dnedot/xdot(psi) = 0         dnedot/xdot(te) = 0
        dni/x(psi) = 0              dni/x(te) = 0
     dni/xdot(psi) = 0           dni/xdot(te) = 0
     dnidot/x(psi) = 0           dnidot/x(te) = 0
    dnidot/xdot(psi) = 0         dnidot/xdot(te) = 0
       dnep/x(psi) = 0             dnep/x(te) = 0
    dnep/xdot(psi) = 0          dnep/xdot(te) = 0
       dnip/x(psi) = 0             dnip/x(te) = 0
    dnip/xdot(psi) = 0          dnip/xdot(te) = 0
        dn1/x(psi) = 0              dn1/x(te) = 0
     dn1/xdot(psi) = 0           dn1/xdot(te) = 0
        dn2/x(psi) = 0              dn2/x(te) = 0
     dn2/xdot(psi) = 0           dn2/xdot(te) = 0
        dn3/x(psi) = 0              dn3/x(te) = 0
     dn3/xdot(psi) = 0           dn3/xdot(te) = 0
        dte/x(psi) = 0              dte/x(te) = 4.64e-09
     dte/xdot(psi) = 0           dte/xdot(te) = 0
        dti/x(psi) = 0              dti/x(te) = 4.64e-09
     dti/xdot(psi) = 0           dti/xdot(te) = 0
        dze/x(psi) = 0              dze/x(te) = 0
     dze/xdot(psi) = 0           dze/xdot(te) = 0
       dtep/x(psi) = 0             dtep/x(te) = 4.99e-09
    dtep/xdot(psi) = 0          dtep/xdot(te) = 0
       dtip/x(psi) = 0             dtip/x(te) = 4.99e-09
    dtip/xdot(psi) = 0          dtip/xdot(te) = 0
     dtedot/x(psi) = 0           dtedot/x(te) = 0
    dtedot/xdot(psi) = 0         dtedot/xdot(te) = 0
     dtidot/x(psi) = 0           dtidot/x(te) = 0
    dtidot/xdot(psi) = 0         dtidot/xdot(te) = 1.02e-16
         dp/x(psi) = 0               dp/x(te) = 0
      dp/xdot(psi) = 0            dp/xdot(te) = 0
         dq/x(psi) = 5.39e-06        dq/x(te) = 0
      dq/xdot(psi) = 0            dq/xdot(te) = 0
      diota/x(psi) = 9.93e-09     diota/x(te) = 0
    diota/xdot(psi) = 0         diota/xdot(te) = 0
     dshear/x(psi) = 5.25e-07    dshear/x(te) = 0
    dshear/xdot(psi) = 0         dshear/xdot(te) = 0
      dlamH/x(psi) = 0            dlamH/x(te) = 0
    dlamH/xdot(psi) = 0         dlamH/xdot(te) = 0
     dwntms/x(psi) = 0           dwntms/x(te) = 0
    dwntms/xdot(psi) = 0         dwntms/xdot(te) = 0
    dwdotntms/x(psi) = 0         dwdotntms/x(te) = 0
    dwdotntms/xdot(psi) = 0         dwdotntms/xdot(te) = 0
      dchie/x(psi) = 5.74e-06     dchie/x(te) = 0
    dchie/xdot(psi) = 0         dchie/xdot(te) = 0
      dchie/u(psi) = 0            dchie/u(te) = 0
      dchii/x(psi) = 5.74e-06     dchii/x(te) = 0
    dchii/xdot(psi) = 0         dchii/xdot(te) = 0
      dchii/u(psi) = 0            dchii/u(te) = 0
       dvpe/x(psi) = 0             dvpe/x(te) = 0
    dvpe/xdot(psi) = 0          dvpe/xdot(te) = 0
       dvpe/u(psi) = 0             dvpe/u(te) = 0
       ddne/x(psi) = 0             ddne/x(te) = 0
    ddne/xdot(psi) = 0          ddne/xdot(te) = 0
       ddne/u(psi) = 0             ddne/u(te) = 0
       dsne/x(psi) = 0             dsne/x(te) = 0
    dsne/xdot(psi) = 0          dsne/xdot(te) = 0
       dsne/u(psi) = 0             dsne/u(te) = 0
      djbsB/x(psi) = 2.76e-07     djbsB/x(te) = 2.4e-08
    djbsB/xdot(psi) = 0         djbsB/xdot(te) = 0
      djbsB/u(psi) = 0            djbsB/u(te) = 0
    dsigneo/x(psi) = 2.07e-07   dsigneo/x(te) = 9.34e-09
    dsigneo/xdot(psi) = 0         dsigneo/xdot(te) = 0
    dsigneo/u(psi) = 0          dsigneo/u(te) = 0
    djphioR/x(psi) = 8.53e-09   djphioR/x(te) = 0
    djphioR/xdot(psi) = 0         djphioR/xdot(te) = 0
    djphioR/u(psi) = 0          djphioR/u(te) = 0
     dpetot/x(psi) = 8.53e-09    dpetot/x(te) = 0
    dpetot/xdot(psi) = 1.02e-10  dpetot/xdot(te) = 0
     dpetot/u(psi) = 9.32e-10    dpetot/u(te) = 2.88e-09
     dpitot/x(psi) = 0           dpitot/x(te) = 0
    dpitot/xdot(psi) = 0         dpitot/xdot(te) = 0
     dpitot/u(psi) = 0           dpitot/u(te) = 0
     djauxB/x(psi) = 0           djauxB/x(te) = 0
    djauxB/xdot(psi) = 0         djauxB/xdot(te) = 0
     djauxB/u(psi) = 0           djauxB/u(te) = 0
    Warning: Exceeded maximum newton iterations
    Warning: exceeded maximum newton iterations
         it  telaps newt     res   t[ms]  dt[ms]  Ip[kA] Icd[kA] Ibs[kA] Ioh[kA]      qe    qmin      q0   Vl[V] Te0[keV] Ti0[keV] ne0[e19]   f_ss
          1   0.073    1 3.2e-01       0       1      80       0    1.41    78.6    15.4    5.41    5.41 4.1e+00    0.22    0.22    1.00 1.8e+01

## "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.