paramconsistency_dke_yp

 Cross-check input parameter consistency

by Y.Peysson (CEA-DRFC) <yves.peysson@cea.fr> and J. Decker (MIT/RLE) <joan.decker@cea.fr>

0001 function [dkeparam] = paramconsistency_dke_yp(dkeparam,equil)
0002 %
0003 % Cross-check input parameter consistency
0004 %
0005 %by Y.Peysson (CEA-DRFC) <yves.peysson@cea.fr> and J. Decker (MIT/RLE) <joan.decker@cea.fr>
0006 %
0007 if isinf(equil.Rp),
0008     dkeparam.bounce_mode = 0;
0009 end
0010 if dkeparam.bounce_mode == 0,
0011     dkeparam.dke_mode = 0;
0012 end
0013 if dkeparam.dke_mode == 0, 
0014     dkeparam.syn_mode =0;
0015 end
0016 if dkeparam.rt_mode == 0,
0017     dkeparam.nit_rf = 1;
0018 end
0019 %
0020 if dkeparam.invproc == -1, 
0021     dkeparam.luproc = -1;
0022 end%MUMPS performs direct inversion only
0023 %
0024 if isempty(dkeparam.psin_S) && isempty(dkeparam.rho_S),%No iteration if one radial point
0025     error('No input radial flux grid is specified !');
0026 end
0027 if length(dkeparam.rho_S) + length(dkeparam.psin_S) == 1,%No iteration if one radial point
0028     dkeparam.rt_mode = 0;
0029     dkeparam.nit_rf = 1;
0030 else
0031     if (dkeparam.ref_mode == 1),dkeparam.ref_mode = 0;end
0032 end
0033 %
0034 if strcmp(equil.id,'LOCAL'),%local calculations, assume flux surface is circular
0035     dkeparam.equil_mode = 0;
0036 end
0037 %
0038 if ~isfield(dkeparam,'clustermode'),
0039     dkeparam.clustermode = '';
0040 end
0041 %
0042 if ~isfield(dkeparam,'mfactor'),
0043     dkeparam.mfactor = 10;
0044 end
0045 %
0046 if ~isfield(dkeparam,'avalanche_mode'),
0047     dkeparam.avalanche_mode = 0;
0048 end
0049 %
0050 if ~isfield(dkeparam,'synchro_mode'),
0051     dkeparam.synchro_mode = 0;
0052 end
0053 %
0054 if ~isfield(dkeparam,'nit_ohm'),
0055     dkeparam.nit_ohm = 1;
0056 elseif dkeparam.nit_ohm > 1,
0057     %
0058     if ~isfield(dkeparam,'tol_ohm'),
0059         dkeparam.tol_ohm = 1e-3;
0060     end
0061     if ~isfield(dkeparam,'efield_ohm'),
0062         dkeparam.efield_ohm = 1e-2;
0063     end
0064     if ~isfield(dkeparam,'Itot_ohm'),
0065         dkeparam.Itot_ohm = NaN;
0066     end
0067     if ~isfield(dkeparam,'Jtot_ohm'),
0068         dkeparam.Jtot_ohm = NaN;
0069         dkeparam.rho_ohm = NaN;
0070     elseif ~isfield(dkeparam,'rho_ohm') || any(size(dkeparam.Jtot_ohm) ~= size(dkeparam.rho_ohm)),
0071         error('error on the current prescription (efield adjusment mode)')
0072     end
0073 end
0074 %
0075 if dkeparam.pnmax_S < 0,
0076     dkeparam.pnmax_S = -dkeparam.pnmax_S;
0077 end
0078 %
0079 if ~isfield(dkeparam,'method'),
0080     dkeparam.method = 'spline';%For the pitch-angle interpolation in the prjectioon of the electron velocity distribution on the Legendre polynomial basis
0081     dkeparam.n_gauleg = 35;%Number of Legendre polynomials for projections of the electron velocity distribution
0082     dkeparam.epsi_gauleg = 1e-14;%Accuracy of the Legendre projection
0083 end
0084 %
0085 if ~isfield(dkeparam,'finitabs'),
0086     dkeparam.finitabs = 0;%The initial distribution is used to first estimate the ray absorption
0087 end
0088 %
0089 if ~isfield(dkeparam,'postest_rf'),
0090     dkeparam.postest_rf = true;%Test that the power absorbed is positive everywhere, or continue convergence process
0091 end
0092 %
0093 if ~isfield(dkeparam,'finsert'),
0094     dkeparam.finsert = 0;%The initial distribution is used to first estimate the ray absorption
0095 end
0096 %
0097 % FP time management
0098 %
0099 if ~isfield(dkeparam,'tn'),
0100     dkeparam.tn = NaN;
0101 end
0102 if ~isfield(dkeparam,'dtn'),
0103     dkeparam.dtn = NaN;
0104 end
0105 %
0106 if isscalar(dkeparam.tn),
0107     if isnan(dkeparam.tn),
0108         if isscalar(dkeparam.dtn) && isnan(dkeparam.dtn),
0109             dkeparam.tn = 100000;%asymptotic calculation - quasi steady-state
0110         elseif isscalar(dkeparam.dtn) && imag(dkeparam.dtn) > 0,%number of ierations fixed by imag(dtn)
0111             nit = imag(dkeparam.dtn);
0112             dkeparam.tn = 100000.*(1:nit)/nit;%
0113         else
0114             dkeparam.tn = cumsum(dkeparam.dtn);%tn fixed by scalar or vector dtn
0115         end
0116     else
0117         if ~isscalar(dkeparam.dtn),
0118             info_dke_yp(2,'Warning : dkeparam.dtn prescription overridden by dkeparam.tn');
0119             dkeparam.dtn = NaN;
0120         end
0121         if isnan(dkeparam.dtn),% automatic time step
0122             nit = min([100,max([1,ceil(60./dkeparam.tn.^(3/8))])]);%number of internal iterations, fitted result from TS/32299 analysis
0123             dkeparam.tn = dkeparam.tn.*(1:nit)/nit;%
0124         elseif imag(dkeparam.dtn) > 0,%number of ierations fixed by imag(dtn)
0125             nit = imag(dkeparam.dtn);
0126             dkeparam.tn = dkeparam.tn.*(1:nit)/nit;%
0127         else% time step determined by dtn
0128             if dkeparam.dtn > dkeparam.tn,
0129                 info_dke_yp(2,['Warning : the prescribed time step (',num2str(dkeparam.dtn),...
0130                 ') is reduced to tn = ',num2str(dkeparam.tn)]);
0131             end
0132             dkeparam.tn = unique([dkeparam.dtn:dkeparam.dtn:dkeparam.tn,dkeparam.tn]);
0133             if length(dkeparam.tn) > 1 && dkeparam.tn(end) - dkeparam.tn(end-1) < dkeparam.dtn,
0134                 info_dke_yp(2,['Warning : the last time step (',num2str(dkeparam.tn(end) - dkeparam.tn(end-1)),...
0135                 ') is shorter than dtn = ',num2str(dkeparam.dtn)]);
0136             end
0137         end
0138     end
0139 else%time steps determined by tn vector
0140     if ~isscalar(dkeparam.dtn) || ~isnan(dkeparam.dtn),
0141         info_dke_yp(2,'Warning : dkeparam.dtn prescription overridden by dkeparam.tn');
0142     end
0143 end
0144 %
0145 dkeparam = rmfield(dkeparam,'dtn');%dtn determined from tn in main_dke_yp
0146 %
0147 if dkeparam.coll_mode ~= 2,
0148     dkeparam.nit_f = 1;%no Legendre iterations without RHS
0149 end
0150 %
0151 if dkeparam.nit_f == 1,
0152     dkeparam.prec0_f = -1;%no residue comparison with only one iteration
0153 end
0154 %
0155 if ~isfield(dkeparam,'timevol'),
0156     dkeparam.timevol = 0;
0157 end
0158 %
0159 if dkeparam.avalanche_mode == 1 || dkeparam.nit_ohm > 1 || imag(dkeparam.norm_mode_f) > 0
0160     dkeparam.timevol = max([dkeparam.timevol,1]);
0161 end
0162 %
0163 if ~isfield(dkeparam,'pnmin0_KO'),
0164     dkeparam.pnmin0_KO = NaN;% equal to pnmax_S
0165 end
0166 if ~isfield(dkeparam,'pnmax1_KO'),
0167     dkeparam.pnmax1_KO = 0;% equal to pnmin2_KO
0168 end
0169 if ~isfield(dkeparam,'pnmin2_KO'),
0170     dkeparam.pnmin2_KO = 1i;% suprathermal domain
0171 end
0172 if ~isfield(dkeparam,'pnmax2_KO'),
0173     dkeparam.pnmax2_KO = NaN;% equal to pnmin0_KO
0174 end
0175 %
0176 if real(dkeparam.norm_mode_f) == 3,%compensate for boundary losses with p=0 source term
0177     dkeparam.timevol = max([dkeparam.timevol,1]);
0178     if ~isfield(dkeparam,'nit_norm'),
0179         dkeparam.nit_norm = 10;
0180     end
0181     if ~isfield(dkeparam,'tol_norm'),
0182         dkeparam.tol_norm = 1e-3;
0183     end
0184     if ~isfield(dkeparam,'thres_norm'),
0185         dkeparam.thres_norm = eps;
0186     end
0187 else
0188     dkeparam.nit_norm = 1;
0189     dkeparam.tol_norm = Inf;    
0190     dkeparam.thres_norm = Inf;
0191 end
0192 %
0193 if ~isfield(dkeparam,'redfac'),
0194     dkeparam.redfac = 0;
0195 end
0196 %
0197 if dkeparam.bounce_mode == 1 && isfield(dkeparam,'mhu_S') && ~isempty(dkeparam.mhu_S),
0198      info_dke_yp(2,'Warning : the pitch angle grid cannot be externally enforced with bounce_mode = 1.');
0199      dkeparam = rmfield(dkeparam,'mhu_S');
0200 end
0201 %
0202 if ~isfield(dkeparam,'opt_load'),
0203     dkeparam.opt_load = 0;
0204 end
0205 %
0206 if ~isfield(dkeparam,'opt_struct'),
0207     dkeparam.opt_struct = 0;
0208 end
0209 %
0210 if ~isfield(dkeparam,'opt_flegendre'),
0211     dkeparam.opt_flegendre = 0;
0212 end
0213 %