# RAPTOR tutorial: ITER scenario with alpha particle heating

In this tutorial we simulate the alpha particle heating for a standard ITER scenario.

## Contents

- ITER scenario parameters
- Actuator inputs
- Kinetic profile settings
- Run RAPTOR
- Plot results

## ITER scenario parameters

    clear; close all hidden;
    run(fullfile(pwd,'..','RAPTOR_path.m')); % add RAPTOR path
    % ramp function
    rampfun = @(t,tmin,ymin,tmax,ymax) max(ymin,min((ymax-ymin)/(tmax-tmin)*(t-tmin),ymax-ymin)+ymin); % anonymous function for ramps
    [config] = RAPTOR_config('ITER'); % load default ITER params

    % RAPTOR run options
    config.numerics.restol = 1e-8;
    config.debug.iterplot = 0;
    config.debug.iterdisp = 10;

    % Physics modules for standard ITER scenario
    % chi_e
    config.chi_e = chi_e('BgB');
    config.chi_e.params.cneo = 0.3;
    config.chi_e.params.aegb = 1.8 * config.chi_e.params.aegb;
    config.chi_e.params.aeb  = 0.3 * config.chi_e.params.aeb;

    config.pbrem.active = 1;
    config.pei.active = 1;
    config.prad.active = 0;
    % ECHCD
    config.echcd = echcd('echcd_gaussian');
    config.echcd.params.active = true;
    config.echcd.params.rdep = [0.05 0.2 0.3]; %
    config.echcd.params.wdep = [.1 .15 .15];  %
    config.echcd.params.cd_eff = [5 5 15];  % current drive efficiency factor
    config.echcd.params.uindices = [3 4 5]; % index of power for each actuator in input vector

    % NBI
    % add physics-based NBI module
    config.nbhcd = nbhcd('nbhcd_gaussian');
    config.nbhcd.params.uindices=[2];
    config.nbhcd.params.rdep = [0.23]; %
    config.nbhcd.params.wdep = [0.15]; % broad heating
    config.nbhcd.params.wdep_out = [1.1]; % broad heating
    config.nbhcd.params.cd_eff = [4]; % current drive
    config.nbhcd.params.active = true;
    config.nbhcd.params.frac_to_electrons = 0.72;

    % switch on alpha heating
    config.palpha.active = 1;
    config.palpha.calpha = 1.0; % factor scaling alpha power (see P_alpha.m for details)
    config.palpha.check = 0; % optional check of alpha particle module derivatives

    % Timing
    t_0    = 0.0; t_end  = 140;

    t_flattop = 100; % flattop current time
    t_nbistart = 85; t_nbimax = 110; % NBI timing
    % H mode
    t_LH = 100; % LH transition time
    d_LH = 10; % LH transition duration

    % define time grid
    Ts1 = 0.2; % fastest sample time
    Ts2 = 1;    % slow sample time

    % time grid
    config.grid.tgrid = unique([t_0  :Ts1 :t_0+1,...
                                t_0+1:Ts2 :t_LH,...
                                t_LH:Ts1:t_LH+2,...
                                t_LH+2:Ts2:t_end]);

    % spatial grid
    nrho=21;
    config.grid.rhogrid = linspace(0,1,nrho);
    config.grid.rhogrid = [0:0.05:0.9,0.925:0.025:1];
    config.te.BCval = 300;
    config.ni.method = 'nescal';
    config.ti.method = 'tescal';

    % Hmode
    config.hmode.active = true;
    config.hmode.KTeped = 1; % feedback gain on pedestal height
    config.hmode.rhoped = .90;
    config.hmode.Teped  = 4.3e3; % pedestal height
    config.hmode.chiefact = 1; % core transport decrease during Hmode
    config.hmode.check  = false;

    config.saw.active=false;

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

## Actuator inputs

plasma current

    Ip0 = 4e6;
    Ipflattop = 12e6;

    U(1,:) = rampfun(params.tgrid,t_0,Ip0,t_flattop,Ipflattop);

    % NBI
    Pnb_max = 33e6;
    U(2,:) = rampfun(params.tgrid,t_nbistart,0,t_nbimax,Pnb_max); % NBI

    % ECCD
    % rdep = [0.05 0.2 0.3]; %
    U(3,:) = interp1([40 60 95 105],[0e6 2e6 0 0],params.tgrid,'linear',0);
    U(4,:) = interp1([40 60 95 105],[5e6 5e6  14e6 0],params.tgrid,'linear',0);
    U(5,:) = interp1([95 105 140],[0 18e6 37e6],params.tgrid,'linear',0);

## Kinetic profile settings

initial conditions

    init.newidth = 2;
    init.ni0scal = 0.9;
    init.niescal = 0.9;
    init.ti0scal = 1;
    init.tiescal = 1;

    init.te0 = 2e3;
    init.tewidth = 0.6;
    init.ze0 = 1.6;
    init.zee = 1.6;
    init.jpow = 1;
    init.Ip0 = U(1,1);

    % Construct kinetic profiles according to the new settings.
    v = build_kinetic_profiles(model,init);
    % Define the initial condition
    x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1));

    % Assign V
    vmod = v*ones(1,numel(params.tgrid));

    % density, assign by hand to V
    neparam = neHmodeEvo; % LH transition effect on ne evolution
    neparam.LHtime = t_LH;
    neparam.LHduration = d_LH; % duration
    neparam.HLtime = t_end+10; % not in simulation

    neparam.peaking = 1;
    neparam.neL_t = rampfun(params.tgrid,t_LH,3e19,t_LH+d_LH,6.5e19);
    neparam.neH_t = 6.5e19;

    ne = neHmodeEvo(model,params,neparam);
    vmod(model.ne.vind,:)    =  model.ne.Lamgauss\ne;

    % hmode timing
    vmod(model.hmode.vind,:) = rampfun(params.tgrid,t_LH,0,t_LH+d_LH,1);

## Run RAPTOR

H-mode simulation

    simres = RAPTOR_predictive(x0,g,vmod,U,model,params);
    % output
    out = RAPTOR_out(simres,model,params);

         it  telaps newt     res   t[ s]  dt[ s]  Ip[MA] Icd[MA] Ibs[MA] Ioh[MA]      qe    qmin      q0   Vl[V] Te0[keV] Ti0[keV] ne0[e19]   f_ss
          1   0.066    4 2.8e-12       0     0.2       4       0   0.154    3.85      15    2.67    2.67 2.3e-01    1.92    1.92    3.00 3.6e+00
         11    0.46    2 6.6e-09       6       1    4.48       0   0.134    4.35    13.4    2.17    2.17 7.1e-01    1.80    1.80    3.00 1.4e+00
         21    0.73    2 2.5e-11      16       1    5.28       0    0.15    5.13    11.4    1.87    1.87 7.8e-01    2.04    2.04    3.00 1.9e+00
         31       1    2 1.3e-11      26       1    6.08       0   0.166    5.91     9.9    1.65    1.65 8.0e-01    2.35    2.35    3.00 2.1e+00
         41     1.3    2 8.8e-12      36       1    6.88       0   0.178     6.7    8.75    1.45    1.45 8.0e-01    2.69    2.69    3.00 2.2e+00
         51     1.7    3 1.6e-11      46       1    7.68   0.111   0.445    7.12    7.84    1.37    1.37 5.6e-01    8.34    8.34    3.00 2.5e+00
         61       2    2 4.0e-11      56       1    8.48   0.162   0.492    7.83     7.1    1.35    1.35 5.1e-01   10.47   10.47    3.00 2.3e+00
         71     2.3    2 3.7e-11      66       1    9.28   0.235   0.536    8.51    6.49    1.21    1.24 4.7e-01   12.19   12.19    3.00 2.2e+00
         81     2.6    2 3.5e-11      76       1    10.1   0.337   0.585    9.16    5.97    1.12    1.12 4.3e-01   13.67   13.67    3.00 2.0e+00
         91     2.9    3 3.1e-11      86       1    10.9   0.457   0.627     9.8    5.53    1.06    1.06 3.9e-01   15.11   15.11    3.00 1.9e+00
        101     3.2    3 5.8e-11      96       1    11.7    2.01   0.773    8.89    5.15    1.03    1.03 2.8e-01   17.96   17.96    3.00 1.4e+00
        111     3.6    3 3.1e-11     101     0.2      12    3.54    1.08    7.38    5.01    1.04    1.04 3.9e-02   17.09   17.09    3.79 1.0e+00
        121       4    3 3.5e-11     108       1      12    4.22    2.05    5.73    5.01    1.05    1.05 -8.2e-03   18.30   18.30    6.36 7.8e-01
        131     4.3    2 4.9e-11     118       1      12    5.06    2.41    4.54    5.01    1.06    1.06 1.3e-02   20.15   20.15    6.50 5.3e-01
        141     4.6    2 8.7e-10     128       1      12    5.52    2.49    3.99    5.01    1.05    1.05 2.0e-02   21.02   21.02    6.50 4.9e-01
        151     4.9    2 2.9e-11     138       1      12    6.02    2.57    3.42    5.01    1.05    1.05 2.2e-02   21.84   21.84    6.50 5.0e-01

## Plot results

    subplot(321); plot(out.time,out.U(1,:)/1e6); title('Ip [MA]');
    subplot(322); plot(out.time,out.U(2,:)/1e6); title('P_{aux} [MW]')
    subplot(323); plot(out.time,out.Palpha(end,:)/1e6); title('P_{alpha} [MW]')
    subplot(324); plot(out.time,out.f_ss); title('f_{ss} (steady stateness)')
    subplot(325); [cs,h] = contour(out.time,out.rho,out.q,[1 3/2 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 16 32],'color','k'); clabel(cs,h,'labelspacing',72); title('T_e contours [keV]')
    xlabel('t [s]'); ylabel('\rho');