rundke

Script for running the DKE solver (can be modified by the user for specific simulations)
by Y.Peysson CEA-DRFC <yves.peysson@cea.fr> and Joan Decker MIT-RLE (jodecker@mit.edu)

0001 %Script for running the DKE solver (can be modified by the user for specific simulations)
0002 %by Y.Peysson CEA-DRFC <yves.peysson@cea.fr> and Joan Decker MIT-RLE (jodecker@mit.edu)
0003 %
0004 clear all
0005 clear mex
0006 clear functions
0007 close all
0008 warning off
0009 global nfig
0010 %
0011 permission = test_permissions_yp;
0012 %
0013 if ~permission 
0014     disp('Please move the script to a local folder where you have write permission before to run it')
0015     return;
0016 end
0017 %
0018 % ***********************This part must be specified by the user, run make files if necessary) *****************************
0019 %
0020 id_simul = 'LH_karney_transp_bounce_JETh77_PETSc';%Simulation ID
0021 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0022 %
0023 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0024 rho_S = [0.05:0.02:0.2,0.22:0.03:0.6,0.62:0.05:0.95];%Normalized radial flux grid where calculations are performed (0 < rho_S < 1) (If one value: local calculation only, not used if empty)
0025 %
0026 id_path = '';%For all paths used by DKE solver
0027 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0028 %
0029 id_equil = 'JETh77';%For plasma equilibrium
0030 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0031 %
0032 id_dkeparam = 'BOUNCE3D_MUMPS';%For DKE code parameters
0033 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0034 %
0035 id_display = 'NO_DISPLAY';%For output code display
0036 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0037 %
0038 id_ohm = '';%For Ohmic electric contribution
0039 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0040 %
0041 ids_wave = {'LH_karney_offaxis_2'};%For RF waves contribution (put all the type of waves needed)
0042 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0043 %
0044 id_transpfaste = 'magneticturbulence';%For fast electron radial transport
0045 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0046 %
0047 id_ripple = '';%For fast electron magnetic ripple losses
0048 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0049 %
0050 %************************************************************************************************************************************
0051 %************************************************************************************************************************************
0052 %************************************************************************************************************************************
0053 %
0054 [dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple] = load_structures_yp('dkepath',id_path,path_path,'equil',id_equil,path_equil,'dkeparam',id_dkeparam,path_dkeparam,'dkedisplay',id_display,path_display,'ohm',id_ohm,path_ohm,'waves',ids_wave,paths_wave,'transpfaste',id_transpfaste,path_transpfaste,'ripple',id_ripple,path_ripple);
0055 %
0056 %************************************************************************************************************************************
0057 %
0058 dkeparam.psin_S = psin_S;
0059 dkeparam.rho_S = rho_S;
0060 %
0061 if isfield(dkepath,'PETSC'),
0062 %
0063 %dkeparam.invproc = -30;%Number of processors
0064 %dkeparam.PETScparam.nproc = 1;%number of processor used by PETSc (depends of the computer framework). Default is 1
0065 %dkeparam.PETScparam.matrixtype = '-matload_type superlu_dist';%Matrix type used by PETSc. For SUPERLU use '-matload_type superlu_dist' otherwise ' '. See PETSc documentation
0066 %dkeparam.PETScparam.kspmethod = '-ksp_type preonly';%KSP method used by PETSc. '-ksp_type preonly': no iteration, direct matrix solve. Otherwise '-ksp_type bicg' or '-ksp_type bcgs'  for example. See PETSc documentation
0067 %dkeparam.PETScparam.pcmethod = '-pc_type lu';%Preconditioning method used by PETSc. '-pc_type lu' for SUPERLU or '-pc-type hypre' for HYPRE.  See PETSc documentation
0068 %
0069     dkeparam.invproc = -30;%Number of processors
0070     dkeparam.PETScparam.nproc = 1;%number of processor used by PETSc (depends of the computer framework). Default is 1
0071     dkeparam.PETScparam.matrixtype = ' ';%Matrix type used by PETSc. For SUPERLU use '-matload_type superlu_dist' otherwise ' '. See PETSc documentation
0072     dkeparam.PETScparam.kspmethod = '-ksp_type bcgs';%KSP method used by PETSc. '-ksp_type preonly': no iteration, direct matrix solve. Otherwise '-ksp_type bicg' or '-ksp_type bcgs'  for example. See PETSc documentation
0073     dkeparam.PETScparam.pcmethod = ' ';%Preconditioning method used by PETSc. '-pc_type lu' for SUPERLU or '-pc-type hypre' for HYPRE.  See PETSc documentation
0074     dkeparam.PETScparam.cloop = 1;%First order collision selfconsistency is done internaly in C (1) or in MatLab (0)
0075 %
0076     equil.pTe = equil.pTe(1) + equil.pTe*0;%Uniform profile enforced for comparison with TScirc case
0077     equil.pne = equil.pne(1) + equil.pne*0;%Uniform profile enforced for comparison with TScirc case
0078     equil.pzni(1,:) = equil.pzni(1,1) + equil.pzni(1,:)*0;%Uniform profile enforced for comparison with TScirc case
0079     equil.pzni(2,:) = equil.pzni(2,1) + equil.pzni(2,:)*0;%Uniform profile enforced for comparison with TScirc case
0080     equil.pzni(3,:) = equil.pzni(3,1) + equil.pzni(3,:)*0;%Uniform profile enforced for comparison with TScirc case
0081     equil.pzni(4,:) = equil.pzni(4,1) + equil.pzni(4,:)*0;%Uniform profile enforced for comparison with TScirc case
0082     equil.pzTi(1,:) = equil.pzTi(1,1) + equil.pzTi(1,:)*0;%Uniform profile enforced for comparison with TScirc case
0083     equil.pzTi(2,:) = equil.pzTi(2,1) + equil.pzTi(2,:)*0;%Uniform profile enforced for comparison with TScirc case
0084     equil.pzTi(3,:) = equil.pzTi(3,1) + equil.pzTi(3,:)*0;%Uniform profile enforced for comparison with TScirc case
0085     equil.pzTi(4,:) = equil.pzTi(4,1) + equil.pzTi(4,:)*0;%Uniform profile enforced for comparison with TScirc case
0086 %
0087     transpfaste.Dr0 = 0;
0088     [Znorm0,Zcurr0,ZP0_0,dke_out0,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0089     disp('D=0 m2/s done');
0090     transpfaste.Dr0 = 0.1;
0091     [Znorm01,Zcurr01,ZP0_01,dke_out01,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0092     disp('D=0.1 m2/s done');
0093     transpfaste.Dr0 = 1;
0094     [Znorm1,Zcurr1,ZP0_1,dke_out1,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0095     disp('D=1 m2/s done');
0096     transpfaste.Dr0 = 10;
0097     [Znorm10,Zcurr10,ZP0_10,dke_out10,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0098     disp('D=10 m2/s done');
0099 %
0100     vparmin = round(1/(waves{1}.rays{1}.sNpar(1)+imag(waves{1}.rays{1}.sdNpar(1))/2)/mksa.betath_ref);
0101     vparmax = round(1/(waves{1}.rays{1}.sNpar(1)-imag(waves{1}.rays{1}.sdNpar(1))/2)/mksa.betath_ref);
0102 %
0103     figure(1),clf,graph1D_jd(gridDKE.xrho,Zcurr0.x_tot_fsav,0,0,'','','',NaN,'','','-','none','r',2,20,gca,1,0.6,0.6);hold on
0104     graph1D_jd(equilDKE.xrho,Zcurr01.x_tot_fsav,0,0,'','','',NaN,'','','-','none','m',1,20,gca,1,0.6,0.6);hold on
0105     graph1D_jd(equilDKE.xrho,Zcurr1.x_tot_fsav,0,0,'','','',NaN,'','','-','none','m',1,20,gca,1,0.6,0.6);hold on
0106     graph1D_jd(equilDKE.xrho,Zcurr10.x_tot_fsav,0,0,'\rho','J_{LH}',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2)],NaN,'','','--','none','b',1,20,gca,1,0.6,0.6);hold on
0107     legend('D_{\psi\psi} = 0 (m^{2}/s)','D_{\psi\psi} = 0.1 (m^{2}/s)','D_{\psi\psi} = 1 (m^{2}/s)','D_{\psi\psi} = 10 (m^{2}/s)')
0108     savefig_dke_yp('profcurr_Dpsipsi')
0109 %
0110 %************************************************************************************************************************************
0111     eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0112     info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0113 %
0114     proc_rundke_yp(1,'DKE_RESULTS_JETh77_LH_karney_transp_bounce_JETh77.mat',1,1)
0115     proc_rundke_yp(1,'DKE_RESULTS_JETh77_LH_karney_transp_bounce_JETh77.mat',14,1)
0116     proc_rundke_yp(1,'DKE_RESULTS_JETh77_LH_karney_transp_bounce_JETh77.mat',24,1)
0117 else
0118     info_dke_yp(2,'PETSc matrix inversion package is not installed on this computer !')
0119 end