# RAPTOR tutorial 3: Customizing heating & current drive

In this tutorial we will add and customize auxiliary heating and current drive systems. We will have to specify the deposition location, as well as the power traces for each actuator.

This tutorial uses simple EC and NBI modules where the deposition is modeled by gaussian distributions. The EC current drive efficiency is scaled by Te/ne, while the NBI current drive efficiency is constant.

## Contents

- Example 1: Add 2 EC sources at fixed location
- Example 2: Time-varying rdep and add NBI
- Example 3: Add gaussian central ICRH
- Example 4: Completely manual specification of heating & current drive

## Example 1: Add 2 EC sources at fixed location

Preparation Start by repeating the first two steps of the previous tutorial

    clear; close all hidden;
    run ../RAPTOR_path.m

    [config] = RAPTOR_config; % load default config

heating & current drive parameters To specify a heating and current drive actuator, we must choose a module

    config.echcd = echcd('echcd_gaussian');
    disp(config.echcd);
    disp(config.echcd.params);

        modeltype: 'echcd_gaussian'
           params: [1x1 struct]

           active: 0
             rdep: [0 0.4000]
             wdep: [0.3000 0.3500]
           cd_eff: [0 1]
        eta_eccd0: 7.0000e+14
          cd_trap: 0
         uindices: [2 3]
            check: 0

Inspecting the result, we see that the config.echcd contains * the name of the module * the default parameter structure associated with the module

Now we can change the default module parameters

    config.echcd.params.active = true; % activate the module, otherwise it does nothing.
    config.echcd.params.rdep = [0 0.4]; % two actuators, set rdep=0 for first, rdep=0.4 for second
    config.echcd.params.wdep = [0.35, 0.35];  % wdep =0.35 for each
    config.echcd.params.cd_eff = [0 1];  % current drive efficiency factor: pure ECH for first, co-ECCD for second
    config.echcd.params.uindices = [2 3]; % index of power for each actuator in input vector

Other modifications make some further modifications to the configuration structure, look in RAPTOR_config.m for exact parameter definitions

    config.grid.tgrid = [0:0.001:0.1]; % time grid
    config.grid.rhogrid = linspace(0,1,21); % spatial grid
    config.debug.iterplot = 0; %10; % plot profiles every N iterations
    config.debug.iterdisp = 10; % display progress every N iterations

    % Create RAPTOR model,params and init structure with new parameters
    [model,params,init,g,v,U] = build_RAPTOR_model(config);

Define inputs We must now define more inputs than the plasma current alone

    % 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));
    % first EC actautor: on-axis heating 1MW starting at 60ms
    U(2,:) = zeros(size(params.tgrid)); U(2,params.tgrid>0.06) = 1e6;
    % second EC actuator power: off-axis current drive, 1MW start at 40ms
    U(3,:) = zeros(size(params.tgrid)); U(3,params.tgrid>0.04) = 1e6;

Initial conditions and run model initial conditions

    init.Ip0 = U(1,1);
    % Define the initial condition for state
    x0 = RAPTOR_initial_conditions(model,init,g,v);
    % run RAPTOR
    simres = RAPTOR_predictive(x0,g,v,U,model,params);

         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.06    4 5.1e-10       0       1      80       0    1.41    78.6    15.4    5.59    5.59 3.4e+00    0.21    0.21    1.00 1.2e+01
         11    0.31    2 2.6e-10      10       1     143       0    4.06     139    8.61    3.19    3.19 5.3e+00    0.40    0.40    1.00 2.0e+01
         21    0.55    3 9.0e-13      20       1     200       0    6.24     194    6.16    2.32    2.32 4.0e+00    0.64    0.64    1.00 1.6e+01
         31    0.77    2 1.8e-11      30       1     200       0    6.75     193    6.16    1.79    1.79 2.7e+00    0.83    0.83    1.00 1.1e+01
         41       1    2 9.2e-13      40       1     200       0    6.67     193    6.16    1.38    1.38 2.4e+00    0.95    0.95    1.00 8.4e+00
         51     1.3    2 2.5e-10      50       1     200    36.1    25.1     139    6.16    1.29    1.29 5.3e-01    3.17    3.17    1.00 2.5e+00
         61     1.5    4 4.3e-14      60       1     200      39    27.2     134    6.16    1.28    1.28 3.8e-01    6.54    6.54    1.00 2.2e+00
         71     1.8    2 2.1e-11      70       1     200    86.8    81.9    31.2    6.16    1.33    1.33 1.2e-01   13.49   13.49    1.00 1.1e+00
         81       2    2 4.9e-14      80       1     200    94.2    89.5    16.3    6.16    1.37    1.37 1.0e-01   14.61   14.61    1.00 1.0e+00
         91     2.3    2 3.8e-14      90       1     200    96.9    93.4    9.66    6.16     1.4     1.4 9.1e-02   15.32   15.32    1.00 9.6e-01
        101     2.5    2 4.4e-14     100       1     200    99.2    97.5    3.28    6.16    1.43    1.44 8.3e-02   16.11   16.11    1.00 9.4e-01

plot outputs

    out=RAPTOR_out(simres,model,params);
    clf;
    subplot(211); surf(out.time,out.rho,out.pauxe);
    xlabel('t [s]');ylabel('\rho'); title('P_{aux,electrons} [W/m^3]'); shading flat;
    subplot(212); surf(out.time,out.rho,out.jaux);
    xlabel('t [s]');ylabel('\rho'); title('j_{aux} [A/m^2]'); shading flat;

The attentive reader will notice that the driven current increases then EC heating is added. This is due to the the ECCD current drive efficiency which increases with the higher temperature.

## Example 2: Time-varying rdep and add NBI

We will now let the deposition locations be specified through the input vector. We also add NBI.

    clear; close all;
    addpath(fullfile(pwd,'..','code'),'-end'); % add RAPTOR path
    [config] = RAPTOR_config; % load default config
    config.grid.tgrid = [0:0.001:0.1]; % time grid
    config.grid.rhogrid = linspace(0,1,21); % spatial grid
    config.debug.iterplot = 0; %10; % plot profiles every N iterations
    config.debug.iterdisp = 10; % display progress every N iterations

    % heating & current drive parameters

EC with time-varying deposition location To define the deposition location as time-varying, set rdep to -1. Then the actuator index assigning the deposition locations are appended to uindex.

    % set ECH/ECCD parameters
    config.echcd.params.active = true;
    config.echcd.params.rdep = [-1 -1]; % -1 to indicate it is commanded extrernally
    config.echcd.params.wdep = [0.3, 0.35]; % deposition width
    config.echcd.params.cd_eff = [0 1]; % CD efficiency
    config.echcd.params.uindices = [2 3 4 5]; % index in input vector
    % 2,3 are for the two powers, 4,5 is for the deposition location

NBI Configure the NBI module like the EC module

    config.nbhcd = nbhcd('nbhcd_gaussian');

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

Create model rerun the config file with the new parameter set

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

Define inputs: Must define additional inputs for rdep and for the NBI power.

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

    % first EC actautor:  1MW starting at 60ms
    U(2,:) = zeros(size(params.tgrid)); U(2,params.tgrid>0.06) = 1e6;
    % second EC actuator power:  1MW start at 40ms
    U(3,:) = zeros(size(params.tgrid)); U(3,params.tgrid>0.04) = 1e6;

    % deposition location for first EC actuator: on-axis with small movement
    U(4,:) = zeros(size(params.tgrid));
    U(4,params.tgrid>0.08) = linspace(0,0.2,sum(params.tgrid>0.08));

    % deposition location for first EC actuator: off-axis with sweep inward
    U(5,:) = 0.4*ones(size(params.tgrid));
    U(5,params.tgrid>0.05) = linspace(0.4,0.2,sum(params.tgrid>0.05));

    % NBI power: ramp slowly
    U(6,:) = 1e6*ones(size(params.tgrid));
    U(6,params.tgrid<0.06) = linspace(0,500e3,sum(params.tgrid<0.06));

    % Second NBI source: off
    U(7,:) = 0e6*ones(size(params.tgrid));

initial condition

    init.Ip0 = U(1,1);
    x0 = RAPTOR_initial_conditions(model,init,g,v);
    % Run RAPTOR
    simres = RAPTOR_predictive(x0,g,v,U,model,params);

         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.058    4 5.1e-10       0       1      80       0    1.41    78.6    15.4    5.59    5.59 3.4e+00    0.21    0.21    1.00 1.2e+01
         11    0.43    2 3.3e-09      10       1     143       0    5.38     138    8.61    3.12    3.12 5.0e+00    0.54    0.54    1.00 2.1e+01
         21     0.7    3 8.5e-13      20       1     200       0    9.94     190    6.16    2.37    2.37 3.6e+00    1.09    1.09    1.00 1.7e+01
         31    0.97    2 4.2e-11      30       1     200       0    12.9     187    6.16    1.91    2.02 2.0e+00    1.69    1.69    1.00 1.0e+01
         41     1.2    2 2.9e-12      40       1     200       0    14.6     185    6.16    1.58    1.78 1.5e+00    2.20    2.20    1.00 7.3e+00
         51     1.5    2 1.2e-10      50       1     200    51.1      40     109    6.16    1.57    1.74 3.2e-01    5.38    5.38    1.00 1.8e+00
         61     1.8    4 4.0e-14      60       1     200    63.1    48.4    88.4    6.16    1.54    1.76 2.0e-01   10.92   10.92    1.00 1.7e+00
         71     2.1    2 2.9e-11      70       1     200     163     140    -103    6.16     1.5     1.9 4.6e-02   22.50   22.50    1.00 1.2e+00
         81     2.3    2 1.2e-12      80       1     200     207     162    -168    6.16    1.46    2.09 4.3e-02   26.01   26.01    1.00 1.2e+00
         91     2.6    2 9.7e-13      90       1     200     250     176    -226    6.16    1.41    2.36 4.3e-02   28.10   28.10    1.00 1.3e+00
        101     2.8    2 1.0e-11     100       1     200     289     184    -273    6.16    1.37     2.9 4.0e-02   29.07   29.07    1.00 1.3e+00

plot results

    out = RAPTOR_out(simres,model,params);
    subplot(321); plot(out.time,out.U(1,:)); title('Ip');
    subplot(322); plot(out.time,out.U([2:3,6],:)); title('P_{aux}')
    subplot(323); plot(out.time,out.U([4:5],:)); title('\rho_{dep}'); set(gca,'Ylim',[0,1]);
    subplot(324); pcolor(out.time,out.rho,bsxfun(@times,out.Vp,out.pec)); shading flat; title('p_{EC}*dV/d\rho')
    subplot(325); [cs,h] = contour(out.time,out.rho,out.q,[1 2 3 4],'color','k'); clabel(cs,h,'labelspacing',72); title('rational q locations');
    xlabel('t [s]'); ylabel('\rho');
    subplot(326); [cs,h] = contour(out.time,out.rho,out.te/1e3,[0.5 1 2 4 8],'color','k'); clabel(cs,h,'labelspacing',72); title('T_e contours [keV]')
    xlabel('t [s]'); ylabel('\rho');

## Example 3: Add gaussian central ICRH

    clear; close all;
    [config] = RAPTOR_config; % load default config
    config.grid.tgrid = [0:0.001:0.1]; % time grid
    config.grid.rhogrid = linspace(0,1,21); % spatial grid
    config.debug.iterplot = 0; %10; % plot profiles every N iterations
    config.debug.iterdisp = 10; % display progress every N iterations

Heating & current drive parameters *keep default EC and NBI, add IC Configure the IC module similarly to EC module

    config.ichcd = ichcd('ichcd_gaussian');

    % set NBI parameters
    config.ichcd.params.active = true;
    config.ichcd.params.rdep = [0];
    config.ichcd.params.wdep = [0.4]; % broad heating
    config.ichcd.params.wdep_out = [0.4]; % broad heating
    config.ichcd.params.cd_eff = [0]; % no current drive
    config.ichcd.params.uindices = [5]; % index in input vector (after NBH = 4 by default)

Create model rerun the config file with the new parameter set

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

Define inputs: Must define additional inputs for rdep and for the NBI power.

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

    % IC power
    U(5,:) = 1e6*ones(size(params.tgrid));
    U(5,params.tgrid<0.06) = linspace(0,1e6,sum(params.tgrid<0.06));

initial condition

    x0 = RAPTOR_initial_conditions(model,init,g,v);
    % Run RAPTOR
    simres = RAPTOR_predictive(x0,g,v,U,model,params);

         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.05    4 5.1e-10       0       1      80       0    1.41    78.6    15.4    5.59    5.59 3.4e+00    0.21    0.21    1.00 1.2e+01
         11    0.41    3 6.0e-15      10       1     143       0    6.98     136    8.61    3.17    3.17 4.7e+00    0.70    0.70    1.00 2.1e+01
         21    0.65    3 6.8e-13      20       1     200       0    14.8     185    6.16    2.53    2.61 3.1e+00    1.62    1.62    1.00 1.6e+01
         31     0.9    2 7.1e-11      30       1     200       0    22.4     178    6.16    1.99    2.46 1.5e+00    2.84    2.84    1.00 8.5e+00
         41     1.2    2 3.4e-12      40       1     200       0    29.2     171    6.16    1.71    2.43 9.8e-01    4.11    4.11    1.00 5.5e+00
         51     1.4    2 3.1e-13      50       1     200       0    35.2     165    6.16    1.56    2.43 6.7e-01    5.30    5.30    1.00 3.8e+00
         61     1.7    2 1.5e-10      60       1     200       0    40.8     159    6.16    1.45    2.44 4.8e-01    6.37    6.37    1.00 2.8e+00
         71     1.9    2 6.3e-14      70       1     200       0    41.6     158    6.16    1.38    2.44 4.2e-01    6.65    6.65    1.00 2.4e+00
         81     2.2    2 3.0e-14      80       1     200       0    40.2     160    6.16    1.32     2.4 3.9e-01    6.62    6.62    1.00 2.2e+00
         91     2.4    2 3.0e-14      90       1     200       0    38.6     161    6.16    1.27    2.33 3.7e-01    6.55    6.55    1.00 2.1e+00
        101     2.7    2 6.2e-14     100       1     200       0    37.3     163    6.16    1.23    2.24 3.5e-01    6.47    6.47    1.00 1.9e+00

plot results

    out = RAPTOR_out(simres,model,params);
    subplot(311); plot(out.time,out.U(1,:)); title('Ip');
    subplot(312); plot(out.time,out.U([5],:)); title('P_{ic}')
    subplot(325); [cs,h] = contour(out.time,out.rho,out.te/1e3,[0.5 1 2 4 8],'color','k'); clabel(cs,h,'labelspacing',72); title('T_e contours [keV]')
    subplot(326); [cs,h] = contour(out.time,out.rho,out.pice/1e3); clabel(cs,h,'labelspacing',72); title('p_{ic} contours [W/m^3]')
    xlabel('t [s]'); ylabel('\rho');

## Example 4: Completely manual specification of heating & current drive

    clear; close all;
    [config] = RAPTOR_config; % load default config
    config.grid.tgrid = [0:0.001:0.1]; % time grid
    config.grid.rhogrid = linspace(0,1,21); % spatial grid
    config.debug.iterplot = 0; %10; % plot profiles every N iterations
    config.debug.iterdisp = 10; % display progress every N iterations

    % Manual heating and current drive specification
    config.echcd = echcd('none');
    config.nbhcd = nbhcd('nbhcd_manual'); % define the manual module
    rho = linspace(0,1,21); % arbitrary rho grid

    %
    disp(nbhcd_manual)

          active: 1
             rho: [21x1 double]
            pnbe: [21x1 double]
            pnbi: [21x1 double]
             jnb: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
        uindices: 2
           check: 0

We now set the nbh deposition parameters manually

    config.nbhcd.params.active = 1;
    config.nbhcd.params.rho  = rho;
    config.nbhcd.params.pnbe  = 1*(1-linspace(0,1,21).^2)'; % heating profile to electrons
    config.nbhcd.params.pnbi  = 0.8*(1-linspace(0,1,21).^2)'; % heating profile to ions (ignored anyway since we don't solve for ions now)
    config.nbhcd.params.jnb = 0.1*config.nbhcd.params.pnbe; % current drive profile
    config.nbhcd.params.uindices = 2; % index of U to be used for scaling the prescribed power deposition profiles
    % *Create model*
    % rerun the config file with the new parameter set
    [model,params,init,g,v] = build_RAPTOR_model(config);
    init.te0 = 300;

    % define inputs
    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));
    % first EC actautor: on-axis heating 1MW starting at 60ms
    U(2,:) = zeros(size(params.tgrid)); U(2,params.tgrid>0.06) = 1e6;

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

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

    % *plot outputs*
    out = RAPTOR_out(simres,model,params);
    figure;
    subplot(221); surf(out.time,out.rho,out.pauxe);
    xlabel('t [s]');ylabel('\rho'); title('P_{aux,electrons} from NBI [W/m^3]'); shading flat;
    subplot(222); surf(out.time,out.rho,out.jaux);
    xlabel('t [s]');ylabel('\rho'); title('j_{aux} NBI [A/m^2]'); shading flat;
    subplot(223); plot(out.time,[out.Pnbe;out.Pnbi]/1e6);
    xlabel('t [s]');ylabel('[MW]'); title('Total power to electrons (-), ions (--)');

    return

         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.075    4 8.7e-15       0       1      80       0     2.9    77.1    15.4    5.61    5.62 3.2e+00    0.26    0.26    1.00 1.3e+01
         11     0.4    2 2.7e-10      10       1     143       0    4.06     139    8.61    3.16    3.16 5.3e+00    0.40    0.40    1.00 2.0e+01
         21    0.71    3 9.1e-13      20       1     200       0    6.24     194    6.16     2.3     2.3 4.0e+00    0.64    0.64    1.00 1.6e+01
         31    0.99    2 1.8e-11      30       1     200       0    6.75     193    6.16    1.78    1.78 2.7e+00    0.83    0.83    1.00 1.1e+01
         41     1.3    2 9.2e-13      40       1     200       0    6.66     193    6.16    1.37    1.37 2.3e+00    0.95    0.95    1.00 8.3e+00
         51     1.6    2 1.5e-13      50       1     200       0    6.56     193    6.16    1.09    1.09 2.1e+00    1.05    1.05    1.00 6.7e+00
         61     1.9    4 1.5e-11      60       1     200       0    6.49     193    6.16   0.917   0.917 1.1e+00    1.46    1.46    1.00 2.6e+00
         71     2.2    2 2.7e-11      70       1     200    13.2    15.7     171    6.16   0.862   0.862 6.0e-01    2.46    2.46    1.00 1.4e+00
         81     2.5    2 4.6e-14      80       1     200    13.2      16     171    6.16   0.828   0.828 5.7e-01    2.52    2.52    1.00 1.2e+00
         91     2.8    2 2.5e-14      90       1     200    13.2    15.9     171    6.16   0.799   0.799 5.5e-01    2.54    2.54    1.00 1.1e+00
        101       3    1 6.9e-09     100       1     200    13.2    15.8     171    6.16   0.775   0.775 5.4e-01    2.56    2.56    1.00 1.0e+00