# Tutorial explaining how to simulate density evolution

The default setting in RAPTOR is to simulate only flux diffusion and electron temperature evolution. Here we explain how to add also an equation for evolving electron density.

## Contents

- Solve for ne with manual transport coefficients, fix ni/ne fraction
- Solve for ne with manual transport coefficients, compute ni from Zeff and quasineutrality
- ITPA-IOS benchmark coordinated by Y.S. Na [paper in preparation]

    % First, add the relevant paths
    close all hidden;
    run(fullfile(pwd,'..','RAPTOR_path.m')); % add RAPTOR path

## Solve for ne with manual transport coefficients, fix ni/ne fraction

Load the default configuration structure.

    [config] = RAPTOR_config; % load default params
    config.equi.mode = 'free'; % required to solve for ne

    % To solve ne equation, set these options:
    config.ne.method = 'state';
    % Set boundary condition type (Dirichlet, specify ne directly) and value
    config.ne.BCtype = 'ne';
    config.ne.BCval =  0.1e19;

    config.ni.method = 'nescal'; % also specify ion density method: scale w.r.t. ne
    init.ni0scal = 0.9; init.niescal = 1; % set core ni to be 90% of core ne and edge ni=ne;

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

    % inspect transport model used for ne;
    disp(model.vpdn_e); % manual model
    disp(params.vpdn_e);

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

    % Run RAPTOR
    simres = RAPTOR_predictive(x0,g,v,U,model,params);
    % To compute some physical outputs, we run the RAPTOR_out function.
    out = RAPTOR_out(simres,model,params);
    % Plots
    hf=figure(1); set(hf,'position',[50 50 700 300]); clf(hf);
    subplot(131); plot(out.time,out.ne/1e19); xlabel('time [s]'); title('n_e at various rho [10^{19}m^{-3}]');
    subplot(132); plot(out.rho,out.ne(:,1:5:end)/1e19); xlabel('\rho_{tor}'); title('n_e profiles [10^{19}m^{-3}]');
    subplot(133); plot(out.time,out.ne(1,:)/1e19,'b',out.time,out.ni(1,:)/1e19,'r--'); xlabel('time [s]'); title('ne(-),ni(--)');

    manual
        dne_pr: 1
        vpe_pr: -5

         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.063    4 3.0e-09       0       1      80       0     1.5    78.5    15.3    5.68    5.68 3.1e+00    0.21    0.21    1.16 8.8e+00
         11    0.38    2 3.9e-12      10       1      80       0    3.33    76.6    15.3    3.27    3.27 2.4e+00    0.35    0.35    1.05 4.4e+00
         21    0.67    2 3.5e-13      20       1      80       0    3.32    76.6    15.3    2.35    2.35 2.1e+00    0.46    0.46    0.80 3.2e+00
         31    0.95    2 1.2e-13      30       1      80       0    3.16    76.8    15.3    1.88    1.88 1.8e+00    0.57    0.57    0.64 2.6e+00
         41     1.2    2 4.2e-14      40       1      80       0    2.98      77    15.3     1.6     1.6 1.7e+00    0.67    0.67    0.54 2.2e+00
         51     1.5    2 1.8e-14      50       1      80       0    2.84    77.1    15.3    1.41    1.41 1.5e+00    0.75    0.75    0.48 1.8e+00

## Solve for ne with manual transport coefficients, compute ni from Zeff and quasineutrality

To solve for ni via Zeff and quasineutrality, we need to introduce one impurity species (since we have two constraints).

    % Change the method for ne to 'qnze' (quasi-neutrality and Z-effective).
    config.ni.method = 'qnze';
    config.n1.method = 'qnze'; % solve for one impurity species as wel.
    config.atom.Z1 = 6; config.atom.A1 = 12; % choose carbon impurity

    [model,params,init,g,v,U] = build_RAPTOR_model(config);
    init.ze0 = 2.0; init.zee = 3.0; % set Zeff at axis and edge.
    v = build_kinetic_profiles(model,init); % update Zeff in v
    x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1));

    % Run RAPTOR
    simres = RAPTOR_predictive(x0,g,v,U,model,params);
    % To compute some physical outputs, we run the RAPTOR_out function.
    out = RAPTOR_out(simres,model,params);

    hf=figure(1); set(hf,'position',[50 50 700 300]); clf(hf);
    subplot(131); plot(out.rho,[out.ne(:,end),out.ni(:,end)]/1e19); xlabel('\rho_{tor}'); title('n_e,n_i [10^{19}m^{-3}]');
    subplot(132); plot(out.rho,out.ze,'r',out.rho,out.ne(:,end)./out.ni(:,end)); xlabel('\rho_{tor}'); title('Z_{eff}, n_i/n_e');
    subplot(133); plot(out.time,out.ne(1,:)/1e19,'b',out.time,out.ni(1,:)/1e19,'r--'); xlabel('time [s]'); title('ne(-),ni(--)');

         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.081    4 2.0e-10       0       1      80       0    1.64    78.3    15.3    5.56    5.56 2.7e+00    0.20    0.20    1.16 9.2e+00
         11    0.43    2 3.2e-12      10       1      80       0    3.06    76.9    15.3    3.07    3.07 2.2e+00    0.31    0.31    1.05 4.7e+00
         21    0.77    2 3.1e-13      20       1      80       0    3.02    76.9    15.3    2.15    2.15 1.9e+00    0.41    0.41    0.80 3.4e+00
         31     1.1    2 1.3e-13      30       1      80       0    2.84    77.1    15.3    1.69    1.69 1.6e+00    0.51    0.51    0.64 2.7e+00
         41     1.4    2 4.5e-14      40       1      80       0    2.66    77.3    15.3    1.43    1.43 1.5e+00    0.60    0.60    0.54 2.3e+00
         51     1.8    2 1.7e-14      50       1      80       0    2.51    77.5    15.3    1.26    1.26 1.4e+00    0.68    0.68    0.48 1.9e+00

## ITPA-IOS benchmark coordinated by Y.S. Na [paper in preparation]

    [config] = RAPTOR_config('ITER'); % load default params
    config.grid.tgrid = [0:.1:10];
    config.grid.rhogrid = linspace(0,1,41);
    config.equi.mode = 'free';

    config.ne.method = 'state';
    config.ne.scal = 1e19;
    config.ne.BCtype = 'ne';
    config.ne.BCval =  4.6e19;

    config.ni.method = 'nescal';

    config.init.ne0 = 6e19;
    config.init.ni0 = 6e19;

    % Transport and sources
    config.vpdn_e = vpdn_e('benchmark');
    config.sn_e = sn_e('benchmark');

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

    x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1));
    simres = RAPTOR_predictive(x0,g,v,U,model,params);
    out = RAPTOR_out(simres,model,params);

    % plots
    hf=figure(1); set(hf,'position',[50 50 700 300]); clf(hf);
    subplot(131); plot(out.rho,out.ne(:,end)/1e19,'b',out.rho,out.ni(:,end)/1e19,'r--');
    legend({'ne','ni'},'location','southwest'); xlabel('\rho_{tor}'); title('Final n_e,n_i profiles'); ylabel('m^{-3}');
    subplot(132); plot(out.sne(:,end)/1e19,'b');
    title('Source'); ylabel('m^{-3}/s'); xlabel('\rho_{tor}'); set(gca,'ylim',[0,5]);
    subplot(133); plot(out.rho,out.dne(:,end),'r',out.rho,out.vne(:,end),'b');
    legend({'D_{ne} [m^2/s]','V_{ne} [m/s]'},'location','northwest'); title('diffusion, pinch')

         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.14    5 6.3e-12       0     0.1       4       0   0.724    3.28    14.9    2.63    2.63 5.1e-02    3.60    3.60    6.01 3.3e+00
         11     1.1    3 3.7e-12       1     0.1       4       0   0.163    3.84    14.9    2.54    2.54 2.9e-01    1.45    1.45    6.54 4.3e+00
         21     1.9    3 4.0e-12       2     0.1       4       0   0.106    3.89    14.9    2.46    2.46 4.1e-01    1.08    1.08    7.50 4.3e+00
         31     2.7    2 4.5e-09       3     0.1       4       0  0.0951     3.9    14.9    2.38    2.38 4.9e-01    1.00    1.00    8.24 3.9e+00
         41     3.4    2 1.4e-09       4     0.1       4       0  0.0905    3.91    14.9    2.31    2.31 5.4e-01    0.97    0.97    8.74 3.5e+00
         51       4    2 4.8e-10       5     0.1       4       0  0.0871    3.91    14.9    2.26    2.26 5.9e-01    0.94    0.94    9.06 3.2e+00
         61     4.6    2 1.7e-10       6     0.1       4       0  0.0842    3.92    14.9    2.24    2.24 6.2e-01    0.92    0.92    9.26 3.1e+00
         71     5.3    2 6.6e-11       7     0.1       4       0  0.0815    3.92    14.9    2.25    2.25 6.5e-01    0.89    0.89    9.40 3.0e+00
         81       6    2 2.5e-11       8     0.1       4       0   0.079    3.92    14.9     2.3     2.3 6.7e-01    0.87    0.87    9.48 3.0e+00
         91     6.6    2 1.1e-11       9     0.1       4       0  0.0767    3.92    14.9    2.37    2.37 7.0e-01    0.84    0.84    9.54 2.9e+00
        101     7.3    2 5.1e-12      10     0.1       4       0  0.0745    3.93    14.9    2.47    2.47 7.1e-01    0.81    0.81    9.57 2.9e+00