diff --git a/.gitignore b/.gitignore index c1a8011..cc72666 100644 --- a/.gitignore +++ b/.gitignore @@ -1,52 +1,57 @@ *.h5 *.h5_old *_genmod.f90 *.mod *.o *.optrpt *.d *.i *.i90 *.out *.m~ dataanalysis/**/* *.optrpt wk/**/*.pdf wk/**/*.m wk/**/*.fig wk/**/*.log wk/**/*.err wk/**/*.out wk/advi/ wk/hotspots*/ wk/threading*/ wk/ahotspots*/ espic2d *.fig *.trace *.log dataanalysis/**/* src/release src/debug wk/advi f2matlab .VSCodeCounter/** matlab/C2xyz_v2 matlab/export_fig # Files for intel profiler wk/espic2d wk/hotspots* wk/ahotspots* wk/threading* *.th *.th.aux *.options *.options1.1 *.cfg *.cfg.1 wk/** + +# CMake build directory +build +build_profile +build_debug diff --git a/.vscode/launch.json b/.vscode/launch.json index 998d867..d8e1fd2 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -1,140 +1,140 @@ { // Use IntelliSense to learn about possible attributes. // Hover to view descriptions of existing attributes. // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 "version": "0.0.1", "configurations": [ { "name": "Python: Current File", "type": "python", "request": "launch", "program": "${file}", "args": [ "/misc/espic2dsimus/openmp/trajectories.h5" ], "console": "integratedTerminal" }, { "name": "Fortran Launch (GDB)", "type": "cppdbg", "request": "launch", "program": "${workspaceFolder}/src/fennecs", "args": ["${workspaceFolder}/wk/no_ela tst/job.in"], "stopAtEntry": false, "cwd": "${workspaceFolder}/wk", //"externalConsole": true, "preLaunchTask": "make debug", "miDebuggerPath": "gdb", }, { "name": "Fortran test.in (GDB)", "type": "cppdbg", "request": "launch", "program": "${workspaceFolder}/src/fennecs", "args": ["${workspaceFolder}/wk/test.in"], "stopAtEntry": false, "cwd": "${workspaceFolder}/wk", //"externalConsole": true, "preLaunchTask": "", "miDebuggerPath": "gdb", "environment": [{"name":"OMP_NUM_THREADS","Value":"1"}] }, { "name": "Fortran gt170.in (GDB)", "type": "cppdbg", "request": "launch", "program": "${workspaceFolder}/src/fennecs", "args": ["${workspaceFolder}/wk/gt170.in"], "stopAtEntry": false, "cwd": "${workspaceFolder}/wk", //"externalConsole": true, "preLaunchTask": "make debug", "miDebuggerPath": "gdb", "environment": [{"name":"OMP_NUM_THREADS","Value":"6"}] }, { "name": "Fortran gt170_coll.in (GDB)", "type": "cppdbg", "request": "launch", "program": "${workspaceFolder}/src/fennecs", "args": ["${workspaceFolder}/wk/gt170_coll.in"], "stopAtEntry": false, "cwd": "${workspaceFolder}/wk", //"externalConsole": true, "preLaunchTask": "make debug", "miDebuggerPath": "gdb", "environment": [{"name":"OMP_NUM_THREADS","Value":"6"}] }, { "name": "Fortran Launch stable.in", "type": "cppdbg", "request": "launch", "program": "${workspaceFolder}/src/fennecs", "args": ["${workspaceFolder}/wk/stable.in"], "stopAtEntry": false, "cwd": "${workspaceFolder}/wk", //"externalConsole": true, "preLaunchTask": "make debug", "miDebuggerPath": "gdb", "environment": [{"name":"OMP_NUM_THREADS","Value":"1"}] }, { "name": "Fortran Launch splinetest.in", "type": "cppdbg", "request": "launch", "program": "${workspaceFolder}/src/fennecs", - "args": ["${workspaceFolder}/wk/splinetest.in"], + "args": ["${workspaceFolder}/wk/Article_test_case/job.in"], "stopAtEntry": false, - "cwd": "${workspaceFolder}/wk", + "cwd": "${workspaceFolder}/wk/Article_test_case", //"externalConsole": true, "preLaunchTask": "make debug", "miDebuggerPath": "gdb", "environment": [{"name":"OMP_NUM_THREADS","Value":"6"}] }, { "name": "Fortran Launch Arc_splines", "type": "cppdbg", "request": "launch", "program": "${workspaceFolder}/src/fennecs", - "args": ["${workspaceFolder}/wk/arc_splines/job.in"], + "args": ["${workspaceFolder}/wk/T-REX/job.in"], "stopAtEntry": false, - "cwd": "${workspaceFolder}/wk/arc_splines", + "cwd": "${workspaceFolder}/wk/T-REX", //"externalConsole": true, "preLaunchTask": "make debug", "miDebuggerPath": "gdb", - "environment": [{"name":"OMP_NUM_THREADS","Value":"6"}] + "environment": [{"name":"OMP_NUM_THREADS","Value":"2"}] }, { "name": "Fortran Launch limup", "type": "cppdbg", "request": "launch", "program": "${workspaceFolder}/src/fennecs", "args": ["${workspaceFolder}/wk/limup/job.in"], "stopAtEntry": false, "cwd": "${workspaceFolder}/wk/limup", //"externalConsole": true, "preLaunchTask": "make debug", "miDebuggerPath": "gdb", "environment": [{"name":"OMP_NUM_THREADS","Value":"6"}] }, { "name": "Fortran Launch ellipellip.in", "type": "cppdbg", "request": "launch", "program": "${workspaceFolder}/src/fennecs", "args": ["${workspaceFolder}/wk/ellipellip.in"], "stopAtEntry": false, "cwd": "${workspaceFolder}/wk", //"externalConsole": true, "preLaunchTask": "make debug", "miDebuggerPath": "gdb", }, { "name": "Intel Debug Attach", "type": "cppdbg", "request": "attach", "processId": "${command:pickProcess}", "program": "${workspaceFolder}/src/fennecs" } ] } diff --git a/.vscode/tasks.json b/.vscode/tasks.json index eaa071f..7ba3d03 100644 --- a/.vscode/tasks.json +++ b/.vscode/tasks.json @@ -1,51 +1,51 @@ { "version": "2.0.0", "tasks": [ { "label": "make debug", "command": "bash", "type": "shell", "options": { "cwd": "${workspaceFolder}/src", "env": { - "PLATFORM": "intel", + "PLATFORM": "intel_spcpc607", } }, "args": [ "-c", "make debug" ], "group": { "kind": "build", "isDefault": true } }, { "label": "make clean", "command": "bash", "type": "shell", "args": [ "-c", "cd ../src && make clean" ], "group": "build" }, { "label": "make release", "command": "bash", "type": "shell", "options": { "cwd": "${workspaceFolder}/src", "env": { "PLATFORM": "intel", } }, "args": [ "-c", "make" ], "group": "build", "problemMatcher": [] } ] } diff --git a/geometries/experiment_arc/exp_arc_withvessel.m b/geometries/experiment_arc/exp_arc_withvessel.m new file mode 100644 index 0000000..9946bb5 --- /dev/null +++ b/geometries/experiment_arc/exp_arc_withvessel.m @@ -0,0 +1,146 @@ +%% Create the geometry 1 for the T-REX experiment +% uses the upper ellipse and constant radius center column. + +C_column_length=400;%mm +C_column_radius=10; %mm +O_electrode_length=400;%mm + +% Outer vacuum vessel +splid=1; +geomcells(splid).Dval=0; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +geomcells(splid).name='vacuum vessel'; +geomcells(splid).type=0; +geomcells(splid).periodic=0; +rmax=50; %mm +lw=linspace(40,rmax,8); +rw=linspace(250,480,20); +uw=flip(linspace(0,rmax,30)); +geomcells(splid).points=flip([250*ones(size(lw')) lw'; + rw' rmax*ones(size(rw')); + 480*ones(size(uw')) uw']/1e3); + + +% Central electrode +splid=2; +geomcells(splid).Dval=-20000; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +%splines(1).epsge=1e-6; +%splines(1).epsce=1e-6; +geomcells(splid).name='vacuum vessel'; +geomcells(splid).type=0; +geomcells(splid).periodic=0; +lw=linspace(-.5,0,10)+C_column_radius; +rw=linspace(270.5,C_column_length-0.5,60); +uw=flip(linspace(0,1,10))*C_column_radius; +geomcells(splid).points=[ 270*ones(size(lw')) lw'; + rw' C_column_radius*ones(size(rw')); + C_column_length*ones(size(uw')) uw']/1e3; + + +% Dielectric ring +splid=4; +geomcells(splid).Dval=0; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +%splines(1).epsge=1e-6; +%splines(1).epsce=1e-6; +geomcells(splid).name='vacuum vessel'; +geomcells(splid).type=2; +geomcells(splid).periodic=0; +lw=flip(linspace(C_column_radius+.1,17.5,10)); +rw=linspace(250,C_column_length-0.5,60); +uw=flip(linspace(0,C_column_radius,10)); +geomcells(splid).points=[ + 250*ones(size(lw')) lw'; + rw' C_column_radius*ones(size(rw')); + C_column_length*ones(size(uw')) uw']/1e3; + +% Bottom plate +splid=3; +geomcells(splid).Dval=0; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +%splines(1).epsge=1e-6; +%splines(1).epsce=1e-6; +geomcells(splid).name='vacuum vessel'; +geomcells(splid).type=0; +geomcells(splid).periodic=0; +rw=linspace(240,250,10); +lw=linspace(C_column_radius-.5,17.5,40); +uw=flip(linspace(240,250,10)); +geomcells(splid).points=flip([ + 250*ones(size(lw')) lw';]/1e3); + + + + + + +% Outer electrode +splid=5; +geomcells(splid).Dval=0; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +%splines(1).epsge=1e-6; +%splines(1).epsce=1e-6; +geomcells(splid).name='vacuum vessel'; +geomcells(splid).type=0; +geomcells(splid).periodic=0; +base=linspace(250,276,80); + +%arc +radius=700; % mm +r_c=716.87399; +z_c=290.26217; + +arcz=linspace(277,O_electrode_length); +arcr=r_c-sqrt(radius^2-(arcz-z_c).^2); + +%tiltedwall + +axialend=linspace(arcr(end)+0.5,39.5,15); +upperend=flip(linspace(250,arcz(end),60)); + +geomcells(splid).points=flip([base' 17*ones(size(base))'; + arcz' arcr'; + arcz(end)*ones(size(axialend')) axialend'; + upperend' 40*ones(size(upperend'))]/1e3); + + + + %% Plots +f=figure; +for k=1:length(geomcells) + plothandle=plot(geomcells(k).points(:,1), geomcells(k).points(:,2),'k-x','linewidth',1.5); + hold on + %geomcells(k).points=[geomcells(k).Z; geomcells(k).R]; + order=geomcells(k).order; + knots=linspace(0,1,size(geomcells(k).points,1)-(order-2)); + knots=augknt(knots, order); + sizec=size(geomcells(k).points,1); + order=length(knots)-sizec(end); + coeffs=geomcells(k).points'; + pp=spmak(knots,coeffs); + s=linspace(0,1,1000); + fittedpos=fnval(pp,s); + plot(fittedpos(1,:),fittedpos(2,:),'x-') +end + + +legend(plothandle,{'Gun geometry'},'location','southwest') +f.PaperUnits='centimeters'; +f.PaperSize=[12,8]; +xlabel('z [m]') +ylabel('r [m]') + +axis equal + +savegeomtoh5('exp_arc_with_vessel_geom.h5',geomcells,1e-2,true); +% print(f,name,'-dpdf','-fillpage') +% savefig(f,name) +% set(f, 'Color', 'w'); +% export_fig(f,name,'-eps') +hold off \ No newline at end of file diff --git a/geometries/experiment_extrude/exp_extrude_withvessel.m b/geometries/experiment_extrude/exp_extrude_withvessel.m new file mode 100644 index 0000000..a993f62 --- /dev/null +++ b/geometries/experiment_extrude/exp_extrude_withvessel.m @@ -0,0 +1,175 @@ +%% Create the geometry 1 for the T-REX experiment +% uses the upper ellipse and constant radius center column. + +% bottom axial position +zbottom=305; % mm + +radialindent=8; %mm + +% Outer vacuum vessel +splid=1; +geomcells(splid).Dval=0; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +geomcells(splid).name='vacuum vessel'; +geomcells(splid).type=0; +geomcells(splid).periodic=0; +rmax=50; %mm +lw=linspace(34,rmax,8); +rwellipse=linspace(zbottom,480,20); +uw=flip(linspace(0,rmax,30)); +geomcells(splid).points=flip([zbottom*ones(size(lw')) lw'; + rwellipse' rmax*ones(size(rwellipse')); + 480*ones(size(uw')) uw']/1e3); + + +% Central electrode +splid=2; +geomcells(splid).Dval=-20000; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +%splines(1).epsge=1e-6; +%splines(1).epsce=1e-6; +geomcells(splid).name='central electrode'; +geomcells(splid).type=0; +geomcells(splid).periodic=0; +lw=linspace(10-radialindent,10,10); +zwellipse=linspace(zbottom+30.5,414.5,150); +rwellipse=10*ones(size(zwellipse)); +uw=flip(linspace(0,10,30)); + + +%tilted ellipse +r_c=10; +z_c=250+125; +cosa=1; +sina=0; +Lz=25; +Lr=2.5; +deltar=(rwellipse-r_c); +deltaz=(zwellipse-z_c); +D=((deltaz)/Lz).^2+((deltar)/Lr).^2; +w=1-D; +ellipsez=zwellipse(w>=0); +deltaz=ellipsez-z_c; +a=(Lr^2*sina^2+Lz^2*cosa^2); +b=2*deltaz*sina*cosa*(Lr^2-Lz^2); +c=-Lr^2*Lz^2 +deltaz.^2*(Lr^2*cosa^2+Lz^2*sina^2); +rwellipse(w>=0)=r_c+(-b-sqrt(b.^2-4*a.*c))/2/a; + +%ellipse= +lowerring=linspace(zbottom,zbottom+29.5,40); +geomcells(splid).points=[ lowerring' (10-radialindent)*ones(size(lowerring')); + (zbottom+30)*ones(size(lw')) lw'; + zwellipse' rwellipse'; + (zwellipse(end)+0.5)*ones(size(uw')) uw']/1e3; + + +% Bottom plate +splid=3; +geomcells(splid).Dval=0; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +%splines(1).epsge=1e-6; +%splines(1).epsce=1e-6; +geomcells(splid).name='vacuum vessel'; +geomcells(splid).type=0; +geomcells(splid).periodic=0; +lw=linspace(10-radialindent,24+radialindent,50); +geomcells(splid).points=flip([ + zbottom*ones(size(lw')) lw';]/1e3); + + + + + + +% Outer electrode +splid=4; +geomcells(splid).Dval=0; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +%splines(1).epsge=1e-6; +%splines(1).epsce=1e-6; +geomcells(splid).name='outer electrode'; +geomcells(splid).type=0; +geomcells(splid).periodic=0; +lbaser=linspace(24+radialindent,24,20); +base=linspace(zbottom+30.5,349,80); +lowerbase=linspace(zbottom,zbottom+29.5,50); + +%tiltedwall +tiltedz=linspace(350,180+250,100); +tiltedr=24*ones(size(tiltedz));%(tiltedz-tiltedz(1))*tan(alpha)+24; + + + +axialend=linspace(tiltedr(end)+0.5,33.5,15); +upperend=flip(linspace(zbottom,180+250,60)); + +geomcells(splid).points=flip([lowerbase' (24+radialindent)*ones(size(lowerbase')); + zbottom+30*ones(size(lbaser')) lbaser'; + base' 24*ones(size(base))'; + tiltedz' tiltedr'; + (180+250)*ones(size(axialend')) axialend'; + upperend' 34*ones(size(upperend'))]/1e3); + + +% Dielectric ring +splid=5; +geomcells(splid).Dval=0; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +%splines(1).epsge=1e-6; +%splines(1).epsce=1e-6; +geomcells(splid).name='dielectric rings'; +geomcells(splid).type=2; +geomcells(splid).periodic=0; +ow=flip(linspace(zbottom,zbottom+30,130)); +lw=flip(linspace(10.1,23.5,10)); +zwellipse=[linspace(zbottom,zwellipse(1)-0.5,30) zwellipse]; +rwellipse=[rwellipse(1)*ones(1,30) rwellipse]; +uw=flip(linspace(0,10,30)); +geomcells(splid).points=[flip([upperend', 34*ones(size(upperend'))]); + flip([(180+250)*ones(size(axialend')) axialend']); + flip([tiltedz' tiltedr']); + flip([base' 24*ones(size(base))']); + ow' 24*ones(size(ow')); + zbottom*ones(size(lw')) lw'; + zwellipse' rwellipse'; + (zwellipse(end)+0.5)*ones(size(uw')) uw']/1e3; + + + %% Plots +f=figure; +for k=1:length(geomcells) + plothandle=plot(geomcells(k).points(:,1), geomcells(k).points(:,2),'k-x','linewidth',1.5); + hold on + %geomcells(k).points=[geomcells(k).Z; geomcells(k).R]; + order=geomcells(k).order; + knots=linspace(0,1,size(geomcells(k).points,1)-(order-2)); + knots=augknt(knots, order); + sizec=size(geomcells(k).points,1); + order=length(knots)-sizec(end); + coeffs=geomcells(k).points'; + pp=spmak(knots,coeffs); + s=linspace(0,1,100000); + fittedpos=fnval(pp,s); + plot(fittedpos(1,:),fittedpos(2,:),'x-') +end + + +legend(plothandle,{'Gun geometry'},'location','southwest') +f.PaperUnits='centimeters'; +f.PaperSize=[12,8]; +xlabel('z [m]') +ylabel('r [m]') + +axis equal + +savegeomtoh5('exp_extrude_with_vessel_geom.h5',geomcells,1e-2,true); +% print(f,name,'-dpdf','-fillpage') +% savefig(f,name) +% set(f, 'Color', 'w'); +% export_fig(f,name,'-eps') +hold off \ No newline at end of file diff --git a/geometries/experiment_upper/exp_outer_withvessel.m b/geometries/experiment_upper/exp_outer_withvessel.m new file mode 100644 index 0000000..88e2f68 --- /dev/null +++ b/geometries/experiment_upper/exp_outer_withvessel.m @@ -0,0 +1,170 @@ +%% Create the geometry 1 for the T-REX experiment +% uses the upper ellipse and constant radius center column. + +% bottom axial position +zbottom=305; % mm + +radialindent=8; %mm + +% Outer vacuum vessel +splid=1; +geomcells(splid).Dval=0; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +geomcells(splid).name='vacuum vessel'; +geomcells(splid).type=0; +geomcells(splid).periodic=0; +rmax=50; %mm +lw=linspace(40,rmax,8); +rw=linspace(zbottom,480,20); +uw=flip(linspace(0,rmax,30)); +geomcells(splid).points=flip([zbottom*ones(size(lw')) lw'; + rw' rmax*ones(size(rw')); + 480*ones(size(uw')) uw']/1e3); + + +% Central electrode +splid=2; +geomcells(splid).Dval=-20000; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +%splines(1).epsge=1e-6; +%splines(1).epsce=1e-6; +geomcells(splid).name='central electrode'; +geomcells(splid).type=0; +geomcells(splid).periodic=0; +lw=linspace(10-radialindent,10,10); +rw=linspace(zbottom+30.5,429.5,60); +uw=flip(linspace(0,10,10)); +%ellipse= +lowerring=linspace(zbottom,zbottom+29.5,40); +geomcells(splid).points=[ lowerring' (10-radialindent)*ones(size(lowerring')); + (zbottom+30)*ones(size(lw')) lw'; + rw' 10*ones(size(rw')); + 430*ones(size(uw')) uw']/1e3; + + +% Bottom plate +splid=3; +geomcells(splid).Dval=0; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +%splines(1).epsge=1e-6; +%splines(1).epsce=1e-6; +geomcells(splid).name='vacuum vessel'; +geomcells(splid).type=0; +geomcells(splid).periodic=0; +lw=linspace(10-radialindent,24+radialindent,50); +geomcells(splid).points=flip([ + zbottom*ones(size(lw')) lw';]/1e3); + + + + + + +% Outer electrode +splid=4; +geomcells(splid).Dval=0; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +%splines(1).epsge=1e-6; +%splines(1).epsce=1e-6; +geomcells(splid).name='outer electrode'; +geomcells(splid).type=0; +geomcells(splid).periodic=0; +lbaser=linspace(24+radialindent,24,20); +base=linspace(zbottom+30.5,349,80); +lowerbase=linspace(zbottom,zbottom+29.5,50); + +%tiltedwall +alpha=0.1745; +tiltedz=linspace(350,180+250,500); +tiltedr=(tiltedz-tiltedz(1))*tan(alpha)+24; + +%tilted ellipse +cosa=cos(alpha); +sina=sin(alpha); +r_c=28; +z_c=375; +Lz=25; +Lr=5; +deltar=(tiltedr-r_c); +deltaz=(tiltedz-z_c); +D=((deltaz*cosa+deltar*sina)/Lz).^2+((deltaz*sina-deltar*cosa)/Lr).^2; +w=1-D; +ellipsez=tiltedz(w>=0); +deltaz=ellipsez-z_c; +a=(Lr^2*sina^2+Lz^2*cosa^2); +b=2*deltaz*sina*cosa*(Lr^2-Lz^2); +c=-Lr^2*Lz^2 +deltaz.^2*(Lr^2*cosa^2+Lz^2*sina^2); +tiltedr(w>=0)=r_c+(-b-sqrt(b.^2-4*a.*c))/2/a; + +axialend=linspace(tiltedr(end)+0.5,39.5,15); +upperend=flip(linspace(zbottom,180+250,60)); + +geomcells(splid).points=flip([lowerbase' (24+radialindent)*ones(size(lowerbase')); + zbottom+30*ones(size(lbaser')) lbaser'; + base' 24*ones(size(base))'; + tiltedz' tiltedr'; + (180+250)*ones(size(axialend')) axialend'; + upperend' 40*ones(size(upperend'))]/1e3); + + +% Dielectric ring +splid=5; +geomcells(splid).Dval=0; % V +geomcells(splid).order=3; +geomcells(splid).dim=2; +%splines(1).epsge=1e-6; +%splines(1).epsce=1e-6; +geomcells(splid).name='dielectric rings'; +geomcells(splid).type=2; +geomcells(splid).periodic=0; +ow=flip(linspace(zbottom,zbottom+30,130)); +lw=flip(linspace(10.1,23.5,10)); +rw=linspace(zbottom,429.5,60); +uw=flip(linspace(0,10,10)); +geomcells(splid).points=[flip([upperend', 40*ones(size(upperend'))]); + flip([(180+250)*ones(size(axialend')) axialend']); + flip([tiltedz' tiltedr']); + flip([base' 24*ones(size(base))']); + ow' 24*ones(size(ow')); + zbottom*ones(size(lw')) lw'; + rw' 10*ones(size(rw')); + 430*ones(size(uw')) uw']/1e3; + + + %% Plots +f=figure; +for k=1:length(geomcells) + plothandle=plot(geomcells(k).points(:,1), geomcells(k).points(:,2),'k-x','linewidth',1.5); + hold on + %geomcells(k).points=[geomcells(k).Z; geomcells(k).R]; + order=geomcells(k).order; + knots=linspace(0,1,size(geomcells(k).points,1)-(order-2)); + knots=augknt(knots, order); + sizec=size(geomcells(k).points,1); + order=length(knots)-sizec(end); + coeffs=geomcells(k).points'; + pp=spmak(knots,coeffs); + s=linspace(0,1,1000); + fittedpos=fnval(pp,s); + plot(fittedpos(1,:),fittedpos(2,:),'x-') +end + + +legend(plothandle,{'Gun geometry'},'location','southwest') +f.PaperUnits='centimeters'; +f.PaperSize=[12,8]; +xlabel('z [m]') +ylabel('r [m]') + +axis equal + +savegeomtoh5('exp_outer_with_vessel_geom.h5',geomcells,1e-2,true); +% print(f,name,'-dpdf','-fillpage') +% savefig(f,name) +% set(f, 'Color', 'w'); +% export_fig(f,name,'-eps') +hold off \ No newline at end of file diff --git a/matlab/@fennecshdf5/displayDensPsi.m b/matlab/@fennecshdf5/displayDensPsi.m new file mode 100644 index 0000000..8346720 --- /dev/null +++ b/matlab/@fennecshdf5/displayDensPsi.m @@ -0,0 +1,79 @@ +function displayDensPsi(obj,tsteps) +%DISPLAYDENSPSI plot the fluid density averaged in time for times t2d(tsteps) +% Displays also the geometry of the metallic boundaries and the +% magnetic field lines +% Show on top the contour of the Psi function for the Davidson +% loading +f=figure('Name', sprintf('%sfluid_dens_psi',obj.name)); +ax1=gca; +if nargin<2 + tsteps=length(obj.t2d); +end + +geomw=obj.geomweight(:,:,1); +dens=mean(obj.N(:,:,tsteps),3); +dens(geomw<=0)=0; +maxdens=max(dens(:)); +[h,curve]=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,dens,30); +hold on +set(curve,'linecolor','none'); +%contour(ax1,obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0],'r-','linewidth',1.5,'Displayname','Boundaries'); +if(obj.conformgeom) + ylim(ax1,[obj.rgrid(1) obj.rgrid(sum(obj.nnr(1:2)))]*1000) +else + ylim(ax1,[obj.rgrid(1) obj.rgrid(end)]*1000) +end +xlim(ax1,[obj.zgrid(1) obj.zgrid(end)]*1000) +xlabel(ax1,'z [mm]') +ylabel(ax1,'r [mm]') +%title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',obj.t2d(fieldstart),obj.t2d(fieldend),double(maxdens))) +c = colorbar(ax1); +c.Label.String= 'electron density [m^{-3}]'; +view(ax1,2) +colormap(flipud(hot)); + +% Draws the magnetic field lines +Blines=obj.rAthet; +levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),20); +Blines(obj.geomweight(:,:,1)<0)=NaN; +[~,h1]=contour(ax1,obj.zgrid*1000,obj.rgrid*1000,Blines,real(levels),'-.','color','k','linewidth',1.5,'Displayname','Magnetic field lines'); + +% Draw the metallic boundaries and the geometry itself +[c1,hContour]=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,-geomw,[0,0],'linewidth',1.5); +drawnow; +xlim(ax1,[obj.zgrid(1)*1000 obj.zgrid(end)*1000]) +% Change the color of the metallic boundaries to grey +hFills=hContour.FacePrims; +[hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' +try + hFills(end).ColorData = uint8([150;150;150;255]); + for idx = 1 : numel(hFills)-1 + hFills(idx).ColorData(4) = 0; % default=255 + end +catch +end + +if(obj.species(1).R.nt<2) + h0=obj.H0; + p0=obj.P0; +else + h0=mean(H(obj,{1:obj.species(1).VR.nparts,obj.species(1).VR.nt,false})); + p0=mean(P(obj,{1:obj.species(1).VR.nparts,obj.species(1).VR.nt,false})); +end + +Mirrorratio=(obj.Rcurv-1)/(obj.Rcurv+1); +locpot=mean(obj.pot(:,:,tsteps),3); +psi=1+obj.qe*locpot/h0-1/(2*obj.me*h0)*(p0./obj.rgrid+obj.qe*0.5*obj.B0.*(obj.rgrid-obj.width/pi*Mirrorratio*cos(2*pi*obj.zgrid'/obj.width).*besseli(1,2*pi*obj.rgrid/obj.width))).^2; +contour(obj.zgrid*1000,obj.rgrid*1000,psi,[0 0],'k--') + +grid on + +f.PaperOrientation='landscape'; +f.PaperUnits='centimeters'; +papsize=[16 14]; +f.PaperSize=papsize; +set(ax1,'fontsize',14) +%axis equal +obj.savegraph(f,sprintf('%sfluid_dens_psi',obj.name)) +end + diff --git a/matlab/@fennecshdf5/fennecshdf5.m b/matlab/@fennecshdf5/fennecshdf5.m index 353ecd0..d1ff868 100644 --- a/matlab/@fennecshdf5/fennecshdf5.m +++ b/matlab/@fennecshdf5/fennecshdf5.m @@ -1,2982 +1,3097 @@ classdef fennecshdf5 %fennecshdf5 General class used to treat hdf5 result files of fennecs code % A result file is loaded with a call to M=fennecshdf5(filename) where filename is the relative or absolute file path % after loading, several quantities and composite diagnostics such as moments of the distribution function or individual particles % quantities can be accessed. properties filename name folder fullpath timestamp info t0d t1d t2d tpart it0 it1 it2 restartsteps restarttimes %% Physical constants vlight=299792458; qe=1.60217662E-19; me=9.109383E-31; eps_0=8.85418781762E-12; kb=1.38064852E-23; %% Run parameters dt % simulation time step nrun % number of time steps simulated nlres nlsave nlclassical % Was the equation of motion solved in the classical framework nlPhis % Was the self-consistent electric field computed nz % number of intervals in the z direction for the grid nnr % number of intervals in the r direction for the grid for each of the 3 mesh regions lz % physical axial dimension of the simulation space nplasma % Number of initial macro particles potinn % Electric potential at the coaxial insert potout % Electric potential at the cylinder surface B0 % Normalization for the magnetic field Rcurv % Magnetic mirror ratio width % Magnetic mirror length n0 % Initial particle density in case of old particle loading temp % Initial particle temperature in case of old particle loading femorder % finite element method order in z and r direction ngauss % Order of the Gauss integration method for the FEM plasmadim % initial dimensions of the plasma for the old particle loading system radii % Radial limits of the three mesh regions coarse,fine,coarse H0 % Initial particle Energy for Davidsons distribution function P0 % Initial particle Angular momentum for Davidsons distribution function normalized % Are the parts quantities normalized in the h5 file nbspecies % Number of species simulated %% Frequencies omepe % Reference plasma frequency used for normalization omece % Reference cyclotronic frequency for normalization %% Normalizations tnorm % Time normalization rnorm % Dimension normalization bnorm % Magnetic field normalization enorm % Electric field normalization phinorm % Electric potential normalization vnorm % Velocity normalization %% Grid data rgrid % Radial grid position points zgrid % Axial grid position points dz % Axial grid step dr % Radial grid step for the three mesh regions CellVol % Volume of the cell used for density calculation celltype % type of cell -1 outside 1 inside 0 border linked_s % location of linked spline bsplinetype %% Magnetic field Br % Radial magnetic field Bz % Axial magnetic field Athet % Azimuthal component of the Magnetic potential vector rAthet % r*Athet used for the representation of magnetic field lines B % Magnetic field amplitude sinthet % ratio to project quantities along the magnetic field lines costhet % ratio to project quantities along the magnetic field lines %% Energies epot % Time evolution of the particles potential energy ekin % Time evolution of the particles kinetic energy etot % Time evolution of the particles total energy etot0 % Time evolution of the reference particle total energy eerr % Time evolution of the error on the energy conservation npart % Time evolution of the number of simulated %% 2D time data evaluated on grid points N % main specie Density fluidUR % main specie radial fluid velocity fluidUZ % main specie axial fluid velocity fluidUTHET % main specie azimuthal fluid velocity pot % Electric potential evaluated at grid points potxt % External Electric potential evaluated at grid points phi % Electric potential in spline form Er % Radial electric field Ez % Axial electric field Erxt % External Radial electric field Ezxt % External Axial electric field Presstens % Pressure tensor fluidEkin % average kinetic energy in each direction %% Splines knotsr % Spline radial knots knotsz % Spline axial knots %% Particle parameters weight % Macro particle numerical weight of the main specie qsim % Macro particle charge msim % Macro particle mass nbparts % Time evolution of the number of simulated particles partepot % Electric potential at the particles positions R % Particles radial position Z % Particles axial position Rindex % Particles radial grid index Zindex % Particles axial grid index partindex % Particles unique id for tracing trajectories VR % Particles radial velocity VZ % Particles axial velocity VTHET % Particles azimuthal velocity THET % Particles azimuthal position species % Array containing the other simulated species %% Celldiag celldiag % Array containing the cell diagnostic data nbcelldiag % Total number of cell diagnostics %% Curvilinear geometry conformgeom % stores if we use the conforming or nonconforming boundary conditions r_a r_b z_r z_0 r_0 r_r L_r L_z Interior above1 above2 interior walltype geomweight dirichletweight gtilde spl_bound %% Maxwell source parameters maxwellsrce %% Collision with neutral parameters neutcol nudcol % effective momentum collision frequency %% Non ideal power supply psupply end methods function file=file(obj) % returns the h5 file name file=obj.filename; end function obj = fennecshdf5(filename,readparts,old) % Reads the new result file filename and read the parts data if readparts==true % adds the helper_classes folder to the path matlabfuncpath = dir([mfilename('fullpath'),'.m']); addpath(sprintf('%s/../helper_classes',matlabfuncpath.folder)); addpath(sprintf('%s/../extrema',matlabfuncpath.folder)); addpath(sprintf('%s/../export_fig',matlabfuncpath.folder)); rehash path % Try catch are there for compatibility with older simulation files filedata=dir(filename); if (isempty(filedata)) error("File: ""%s"" doesn't exist",filename) end obj.folder=filedata.folder; obj.filename=filename; [~, obj.name, ext] = fileparts(obj.filename); obj.filename=[obj.name,ext]; obj.fullpath=[obj.folder,'/',obj.filename]; obj.timestamp=filedata.date; if nargin==1 readparts=true; end if nargin<3 old=false; end %obj.info=h5info(filename); %% Read the run parameters obj.dt = h5readatt(obj.fullpath,'/data/input.00/','dt'); obj.nrun = h5readatt(obj.fullpath,'/data/input.00/','nrun'); obj.nlres = strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlres'),'y'); obj.nlsave = strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlsave'),'y'); obj.nlclassical =strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlclassical'),'y'); obj.nlPhis =strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlPhis'),'y'); obj.nz = h5readatt(obj.fullpath,'/data/input.00/','nz'); obj.nnr = h5read(obj.fullpath,'/data/input.00/nnr'); obj.lz = h5read(obj.fullpath,'/data/input.00/lz'); obj.qsim = h5readatt(obj.fullpath,'/data/input.00/','qsim'); obj.msim = h5readatt(obj.fullpath,'/data/input.00/','msim'); try obj.r_a=h5readatt(obj.fullpath,'/data/input.00/geometry','r_a'); obj.r_b=h5readatt(obj.fullpath,'/data/input.00/geometry','r_b'); obj.z_r=h5readatt(obj.fullpath,'/data/input.00/geometry','z_r'); obj.r_r=h5readatt(obj.fullpath,'/data/input.00/geometry','r_r'); obj.r_0=h5readatt(obj.fullpath,'/data/input.00/geometry','r_0'); obj.z_0=h5readatt(obj.fullpath,'/data/input.00/geometry','z_0'); obj.above1=h5readatt(obj.fullpath,'/data/input.00/geometry','above1'); obj.above2=h5readatt(obj.fullpath,'/data/input.00/geometry','above2'); obj.interior=h5readatt(obj.fullpath,'/data/input.00/geometry','interior'); obj.walltype=h5readatt(obj.fullpath,'/data/input.00/geometry','walltype'); try obj.L_r=h5readatt(obj.fullpath,'/data/input.00/geometry','L_r'); obj.L_z=h5readatt(obj.fullpath,'/data/input.00/geometry','L_z'); catch end obj.conformgeom=false; catch obj.conformgeom=true; obj.walltype=0; obj.r_a=obj.rgrid(1); obj.r_b=obj.rgrid(end); obj.above1=1; obj.above2=-1; obj.L_r=0; obj.L_z=0; end try obj.weight=h5readatt(obj.fullpath,'/data/part/','weight'); catch obj.weight=obj.msim/obj.me; end filesgrpinfo=h5info(obj.fullpath,'/files'); nbrst=h5readatt(obj.fullpath,'/files','jobnum'); obj.restartsteps(1)=0; obj.restarttimes(1)=0; grp=sprintf('/data/input.%02i/',0); obj.dt(1)=h5readatt(obj.fullpath,grp,'dt'); for i=1:nbrst - grp=sprintf('/data/input.%02i/',i); - obj.restartsteps(i+1)= h5readatt(obj.fullpath,grp,'startstep'); - obj.restarttimes(i+1)= obj.restarttimes(i) + (obj.restartsteps(i+1)-obj.restartsteps(i))*obj.dt(i); - obj.dt(i+1)=h5readatt(obj.fullpath,grp,'dt'); + grp=sprintf('/data/input.%02i/',i); + obj.restartsteps(i+1)= h5readatt(obj.fullpath,grp,'startstep'); + obj.restarttimes(i+1)= obj.restarttimes(i) + (obj.restartsteps(i+1)-obj.restartsteps(i))*obj.dt(i); + obj.dt(i+1)=h5readatt(obj.fullpath,grp,'dt'); end obj.nplasma = h5readatt(obj.fullpath,'/data/input.00/','nplasma'); obj.potinn = h5readatt(obj.fullpath,'/data/input.00/','potinn'); obj.potout = h5readatt(obj.fullpath,'/data/input.00/','potout'); obj.B0 = h5readatt(obj.fullpath,'/data/input.00/','B0'); obj.Rcurv = h5readatt(obj.fullpath,'/data/input.00/','Rcurv'); obj.width = h5readatt(obj.fullpath,'/data/input.00/','width'); obj.n0 = h5readatt(obj.fullpath,'/data/input.00/','n0'); obj.temp = h5readatt(obj.fullpath,'/data/input.00/','temp'); try obj.it0 = h5readatt(obj.fullpath,'/data/input.00/','it0d'); obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it2d'); obj.it2 = h5readatt(obj.fullpath,'/data/input.00/','itparts'); catch obj.it0 = h5readatt(obj.fullpath,'/data/input.00/','it0'); obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it1'); obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it2'); end try try obj.nbspecies=h5readatt(obj.fullpath,'/data/part/','nbspecies'); catch obj.nbspecies=h5readatt(obj.fullpath,'/data/input.00/','nbspecies'); end obj.normalized=strcmp(h5readatt(obj.fullpath,'/data/input.00/','rawparts'),'y'); catch obj.nbspecies=1; obj.normalized=false; end try obj.nbcelldiag=h5readatt(obj.fullpath,'/data/celldiag/','nbcelldiag'); catch obj.nbcelldiag=0; end obj.omepe=sqrt(abs(obj.n0)*obj.qe^2/(obj.me*obj.eps_0)); obj.omece=obj.qe*obj.B0/obj.me; obj.npart= h5read(obj.fullpath, '/data/var0d/nbparts'); try obj.nudcol= h5read(obj.fullpath, '/data/var0d/nudcol'); catch end try obj.H0 = h5read(obj.fullpath,'/data/input.00/H0'); obj.P0 = h5read(obj.fullpath,'/data/input.00/P0'); catch obj.H0=3.2e-14; obj.P0=8.66e-25; end % Normalizations if old obj.tnorm=abs(1/obj.omepe); else obj.tnorm=min(abs(1/obj.omepe),abs(1/obj.omece)); end obj.rnorm=obj.vlight*obj.tnorm; obj.bnorm=obj.B0; obj.enorm=obj.vlight*obj.bnorm; obj.phinorm=obj.enorm*obj.rnorm; obj.vnorm=obj.vlight; % Conversion of the bias to V obj.potinn=obj.potinn*obj.phinorm; obj.potout=obj.potout*obj.phinorm; % Grid data obj.rgrid= h5read(obj.fullpath, '/data/var1d/rgrid')*obj.rnorm; obj.zgrid= h5read(obj.fullpath, '/data/var1d/zgrid')*obj.rnorm; obj.dz=(obj.zgrid(end)-obj.zgrid(1))/double(obj.nz); rid=1; for i=1:length(obj.nnr) obj.dr(i)=(obj.rgrid(sum(obj.nnr(1:i))+1)-obj.rgrid(rid))/double(obj.nnr(i)); rid=rid+obj.nnr(i); end Br = h5read(obj.fullpath,'/data/fields/Br')*obj.bnorm; obj.Br= reshape(Br,length(obj.zgrid),length(obj.rgrid)); Bz = h5read(obj.fullpath,'/data/fields/Bz')*obj.bnorm; obj.Bz= reshape(Bz,length(obj.zgrid),length(obj.rgrid)); try - Atheta = h5read(obj.fullpath,'/data/fields/Athet')*obj.bnorm; + Atheta = h5read(obj.fullpath,'/data/fields/Athet'); obj.Athet= reshape(Atheta,length(obj.zgrid),length(obj.rgrid)); [rmeshgrid,~]=meshgrid(obj.rgrid,obj.zgrid); obj.rAthet=(rmeshgrid.*obj.Athet)'; catch end obj.B=sqrt(obj.Bz.^2+obj.Br.^2); obj.costhet=(obj.Br./obj.B)'; obj.sinthet=(obj.Bz./obj.B)'; clear Br Bz try obj.t0d=h5read(obj.fullpath,'/data/var0d/time'); catch obj.t0d=obj.dt.*double(0:length(obj.epot)-1); end try for i=0:nbrst grp=sprintf('/data/input.%02i/',i); obj.Erxt(:,:,i+1)=reshape(h5read(obj.fullpath,[grp,'Erxt']),length(obj.zgrid),length(obj.rgrid))'*obj.enorm; obj.Ezxt(:,:,i+1)=reshape(h5read(obj.fullpath,[grp,'Ezxt']),length(obj.zgrid),length(obj.rgrid))'*obj.enorm; obj.potxt(:,:,i+1)=reshape(h5read(obj.fullpath,[grp,'potxt']),length(obj.zgrid),length(obj.rgrid))'*obj.phinorm; end catch end obj.femorder = h5read(obj.fullpath,'/data/input.00/femorder'); obj.ngauss = h5read(obj.fullpath,'/data/input.00/ngauss'); obj.plasmadim = h5read(obj.fullpath,'/data/input.00/plasmadim'); obj.radii = h5read(obj.fullpath,'/data/input.00/radii'); obj.epot = h5read(obj.fullpath,'/data/var0d/epot'); obj.ekin = h5read(obj.fullpath,'/data/var0d/ekin'); obj.etot = h5read(obj.fullpath,'/data/var0d/etot'); try obj.etot0 = h5read(obj.fullpath,'/data/var0d/etot0'); obj.eerr = obj.etot-obj.etot0; catch obj.eerr = obj.etot-obj.etot(2); end if(obj.normalized) obj.pot=gridquantity(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,1); obj.Er=gridquantity(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,1); obj.Ez=gridquantity(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,1); else obj.pot=gridquantity(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,obj.phinorm); obj.Er=gridquantity(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,obj.enorm); obj.Ez=gridquantity(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,obj.enorm); end try obj.t2d = h5read(obj.fullpath,'/data/fields/time'); catch info=h5info(obj.fullpath,'/data/fields/partdensity'); obj.t2d=obj.dt*(0:info.objspace.Size(2)-1)*double(obj.it1); end try info=h5info(obj.fullpath,'/data/fields/moments'); obj.femorder = h5read(obj.fullpath,'/data/input.00/femorder'); kr=obj.femorder(2)+1; obj.knotsr=augknt(obj.rgrid,kr); kz=obj.femorder(1)+1; obj.knotsz=augknt(obj.zgrid,kz); try obj.CellVol= reshape(h5read(obj.fullpath,'/data/fields/volume'),length(obj.knotsz)-kz,length(obj.knotsr)-kr); obj.CellVol=permute(obj.CellVol,[2,1,3])*obj.rnorm^3; catch zvol=fnder(spmak(obj.knotsz,ones(1,length(obj.knotsz)-kz)), -1 ); rvol=fnder(spmak(obj.knotsr,2*pi*[obj.rgrid' 2*obj.rgrid(end)-obj.rgrid(end-1)]), -1 ); ZVol=diff(fnval(zvol,obj.knotsz)); RVol=diff(fnval(rvol,obj.knotsr)); obj.CellVol=RVol(3:end-1)*ZVol(3:end-1)'; obj.CellVol=padarray(obj.CellVol,[1,1],'replicate','post'); end try obj.geomweight = h5read(obj.fullpath,'/data/input.00/geometry/geomweight'); obj.geomweight= reshape(obj.geomweight,length(obj.zgrid),length(obj.rgrid),[]); obj.geomweight = permute(obj.geomweight,[2,1,3]); catch obj.geomweight=ones(length(obj.rgrid),length(obj.zgrid),3); end try obj.dirichletweight = h5read(obj.fullpath,'/data/input.00/geometry/dirichletweight'); obj.dirichletweight= reshape(obj.dirichletweight,length(obj.zgrid),length(obj.rgrid),[]); obj.dirichletweight = permute(obj.dirichletweight,[2,1,3]); catch obj.dirichletweight=obj.geomweight; end try obj.gtilde = h5read(obj.fullpath,'/data/input.00/geometry/gtilde'); obj.gtilde= reshape(obj.gtilde,length(obj.zgrid),length(obj.rgrid),[]); obj.gtilde = permute(obj.gtilde,[2,1,3]); catch obj.gtilde=zeros(length(obj.rgrid),length(obj.zgrid),3); end geomweight=ones(length(obj.rgrid),length(obj.zgrid)); if(obj.normalized) obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, 1, geomweight, 1); obj.phi=splinequantity(obj.fullpath,'/data/fields/phi', obj.knotsr, obj.knotsz, obj.femorder, 1, obj.geomweight(:,:,1), -1); else obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, abs(obj.qsim/obj.qe), geomweight, 1); end obj.fluidUR=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 2); obj.fluidUTHET=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 3); obj.fluidUZ=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 4); if(obj.normalized) obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.me, geomweight); obj.fluidEkin=splineenergy(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.me*0.5, geomweight); else obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim, geomweight); obj.fluidEkin=splineenergy(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim*0.5, geomweight); end try obj.celltype=h5read(obj.fullpath,'/data/input.00/geometry/ctype')'; obj.linked_s=h5read(obj.fullpath,'/data/input.00/geometry/linked_s'); obj.bsplinetype=h5read(obj.fullpath,'/data/input.00/geometry/bsplinetype'); obj.bsplinetype=reshape(obj.bsplinetype,length(obj.knotsz)-kz,length(obj.knotsr)-kr); catch obj.celltype=[]; obj.linked_s=[]; end catch obj.CellVol=(obj.zgrid(2:end)-obj.zgrid(1:end-1))*((obj.rgrid(2:end).^2-obj.rgrid(1:end-1).^2)*pi)'; obj.CellVol=obj.CellVol'; obj.N=griddensity(obj.fullpath, '/data/fields/partdensity', sum(obj.nnr)+1, obj.nz+1, obj.CellVol, abs(obj.qsim/obj.qe), true); obj.fluidUR=gridquantity(obj.fullpath, '/data/fields/fluidur', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true); obj.fluidUTHET=gridquantity(obj.fullpath, '/data/fields/fluiduthet', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true); obj.fluidUZ=gridquantity(obj.fullpath, '/data/fields/fluiduz', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true); end % If we have a maxwellian source, read its parameters try obj.maxwellsrce.rlim=h5read(obj.fullpath, '/data/input.00/maxwellsource/rlimits'); obj.maxwellsrce.zlim=h5read(obj.fullpath, '/data/input.00/maxwellsource/zlimits'); obj.maxwellsrce.frequency=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','frequency'); obj.maxwellsrce.radialtype=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','radialtype'); obj.maxwellsrce.temperature=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','temperature'); obj.maxwellsrce.time_end=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','time_end'); obj.maxwellsrce.time_start=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','time_start'); obj.maxwellsrce.vth=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','vth'); obj.maxwellsrce.rate=obj.maxwellsrce.frequency*obj.weight/(pi*(diff(obj.maxwellsrce.rlim.^2))*diff(obj.maxwellsrce.zlim)); obj.maxwellsrce.current=obj.maxwellsrce.frequency*obj.weight*obj.qe; obj.maxwellsrce.present=true; catch obj.maxwellsrce.present=false; end %% load neutcol parameters try obj.neutcol.neutdens=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','neutdens')); obj.neutcol.neutpressure=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','neutpressure')); obj.neutcol.scatter_fac=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','scatter_fac')); obj.neutcol.Eion=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','Eion')); obj.neutcol.E0=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','E0')); obj.neutcol.Escale=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','Escale')); try obj.neutcol.io_cross_sec=double(h5read(obj.fullpath, '/data/input.00/neutcol/io_cross_sec')); obj.neutcol.io_cross_sec(:,2)=obj.neutcol.io_cross_sec(:,2)*obj.rnorm^2; obj.neutcol.io_cross_sec(:,3)=[log(obj.neutcol.io_cross_sec(2:end,2)./obj.neutcol.io_cross_sec(1:end-1,2))... ./log(obj.neutcol.io_cross_sec(2:end,1)./obj.neutcol.io_cross_sec(1:end-1,1)); 0]; obj.neutcol.iom_cross_sec=zeros(500,3); obj.neutcol.iom_cross_sec(:,1)=logspace(log10(obj.neutcol.Eion+0.001),log10(5e4),size(obj.neutcol.iom_cross_sec,1)); obj.neutcol.iom_cross_sec(:,2)=obj.sigmiopre(obj.neutcol.iom_cross_sec(:,1),true); obj.neutcol.iom_cross_sec(:,3)=abs([log(obj.neutcol.iom_cross_sec(2:end,2)./obj.neutcol.iom_cross_sec(1:end-1,2))... ./log(obj.neutcol.iom_cross_sec(2:end,1)./obj.neutcol.iom_cross_sec(1:end-1,1)); 0]); catch obj.neutcol.io_cross_sec=[]; obj.neutcol.iom_cross_sec=[]; end try obj.neutcol.ela_cross_sec=double(h5read(obj.fullpath, '/data/input.00/neutcol/ela_cross_sec')); obj.neutcol.ela_cross_sec(:,2)=obj.neutcol.ela_cross_sec(:,2)*obj.rnorm^2; obj.neutcol.ela_cross_sec(:,3)=[log(obj.neutcol.ela_cross_sec(2:end,2)./obj.neutcol.ela_cross_sec(1:end-1,2))... ./log(obj.neutcol.ela_cross_sec(2:end,1)./obj.neutcol.ela_cross_sec(1:end-1,1)); 0]; catch obj.neutcol.ela_cross_sec=[]; end obj.neutcol.present=true; catch obj.neutcol.present=false; end %% load spline boundaries try obj.spl_bound.nbsplines=h5readatt(obj.fullpath, '/data/input.00/geometry_spl','nbsplines'); for i=1:obj.spl_bound.nbsplines splgroup=sprintf('/data/input.00/geometry_spl/%02d',i); obj.spl_bound.boundary(i).knots=h5read(obj.fullpath,sprintf('%s/knots',splgroup)); obj.spl_bound.boundary(i).Dval=h5readatt(obj.fullpath,splgroup,'Dirichlet_val'); obj.spl_bound.boundary(i).coefs=reshape(h5read(obj.fullpath,sprintf('%s/pos',splgroup)),2,[])'; obj.spl_bound.boundary(i).order=h5readatt(obj.fullpath,splgroup,'order'); obj.spl_bound.boundary(i).kind=h5readatt(obj.fullpath,splgroup,'kind'); try obj.spl_bound.boundary(i).type=h5readatt(obj.fullpath,splgroup,'type'); catch end obj.spl_bound.boundary(i).fun=spmak(obj.spl_bound.boundary(i).knots,obj.spl_bound.boundary(i).coefs'); end catch obj.spl_bound.nbsplines=0; end - %% load non ideal power supply parameters + %% load non ideal power supply parameters try obj.psupply.targetbias=h5readatt(obj.fullpath, '/data/input.00/psupply','targetbias'); obj.psupply.expdens=h5readatt(obj.fullpath, '/data/input.00/psupply','expdens'); obj.psupply.PSresistor=h5readatt(obj.fullpath, '/data/input.00/psupply','PSresistor'); obj.psupply.geomcapacitor=h5readatt(obj.fullpath, '/data/input.00/psupply','geomcapacitor'); obj.psupply.nbhdt=h5readatt(obj.fullpath, '/data/input.00/psupply','nbhdt'); obj.psupply.biases=h5read(obj.fullpath, '/data/var0d/biases'); obj.psupply.current=h5read(obj.fullpath, '/data/var0d/current'); obj.psupply.tau=obj.psupply.PSresistor*obj.psupply.geomcapacitor*obj.psupply.expdens/obj.neutcol.neutdens; obj.psupply.active=true; obj.psupply.bdpos=h5read(obj.fullpath, '/data/input.00/psupply/bdpos'); catch obj.psupply.active=false; end obj.species=h5parts.empty(obj.nbspecies,0); % Read the main particles parameters if(readparts) obj.species(1)=h5parts(obj.fullpath,'/data/part',obj,obj.normalized); - + if(obj.nbspecies >1) for i=2:obj.nbspecies obj.species(i)=h5parts(obj.fullpath,sprintf('/data/part/%2d',i),obj,true); end end end try obj.tpart = h5read(obj.fullpath,'/data/part/time'); obj.nbparts = h5read(obj.fullpath,'/data/part/Nparts'); catch obj.nbparts=obj.npart; obj.tpart=obj.dt*(0:size(obj.species(1).R,2)-1)*double(obj.it2); end if(obj.nbcelldiag > 0) obj.celldiag=h5parts.empty; j=0; for i=1:obj.nbcelldiag nbparts=h5read(obj.fullpath,sprintf('%s/Nparts',sprintf('/data/celldiag/%02d',i))); if (sum(nbparts)>0) j=j+1; obj.celldiag(j)=h5parts(obj.fullpath,sprintf('/data/celldiag/%02d',i),obj,true); obj.celldiag(j).rindex=double(h5readatt(obj.fullpath, sprintf('/data/celldiag/%02d',i),'rindex'))+(1:2); obj.celldiag(j).zindex=double(h5readatt(obj.fullpath, sprintf('/data/celldiag/%02d',i),'zindex'))+(1:2); end end end end %------------------------------------------ % Functions for accesing secondary simulation quantities function Atheta=Atheta(obj,R,Z) %% returns the magnetic vector potential at position R,Z interpolated from stored Athet in h5 file % halflz=(obj.zgrid(end)+obj.zgrid(1))/2; % Atheta=0.5*obj.B0*(R-obj.width/pi*(obj.Rcurv-1)/(obj.Rcurv+1)... % .*besseli(1,2*pi*R/obj.width).*cos(2*pi*(Z-halflz)/obj.width)); - Atheta=interp2(obj.rgrid,obj.zgrid,obj.Athet,R,Z); + Atheta=interp2(obj.rgrid,obj.zgrid,obj.Athet,R,Z,'spline'); end function quantity=H(obj,indices) %% computes the total energy for the main specie % for the particle with index indices{1} at timepart step indices{2} % which is time obj.timepart(indices{2}) if strcmp(indices{1},':') p=1:obj.species(1).VR.nparts;% if nothing is defined we load all particles else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); %if nothing is defined all time steps are considered else t=indices{2}; end % if track is true we look at specific particles with their % index and follow them in time % if it is false we just care about the distribution function % and specific particles can have different positions in the % resulting array for each timestep if size(indices,1)>2 track=indices{3}; else track=false; end quantity=0.5*obj.me*(obj.species(1).VR(p,t,track).^2+obj.species(1).VTHET(p,t,track).^2+obj.species(1).VZ(p,t,track).^2)+obj.species(1).partepot(p,t,track); end function quantity=P(obj,indices) %P computes the canonical angular momentum for the main specie % for the particle with index indices{1} at timepart step indices{2} % which is time obj.timepart(indices{2}) if strcmp(indices{1},':') p=1:obj.species(1).R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end % if track is true we look at specific particles with their % index and follow them in time % if it is false we just care about the distribution function % and specific particles can have different positions in the % resulting array for each timestep if size(indices,1)>2 track=indices{3}; else track=false; end quantity=obj.species(1).R(p,t,track).*(obj.species(1).VTHET(p,t,track)*obj.me+sign(obj.qsim)*obj.qe*obj.Atheta(obj.species(1).R(p,t,track),obj.species(1).Z(p,t,track))); end function quantity=Vpar(obj,varargin) %Vpar Computes the parallel velocity for the main specie % for the particle with index indices{1} at timepart step indices{2} % which is time obj.timepart(indices{2}) if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if strcmp(indices{1},':') p=1:obj.species(1).R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end % if track is true we look at specific particles with their % index and follow them in time % if it is false we just care about the distribution function % and specific particles can have different positions in the % resulting array for each timestep if size(indices,1)>2 track=indices{3}; else track=false; end Zp=obj.species(1).Z(p,t,track);% get the particle axial positon Rp=obj.species(1).R(p,t,track);% get the particle radial position % interpolate the magnetic field at the particle position Bzp=interp2(obj.zgrid,obj.rgrid,obj.Bz',Zp,Rp,'makima'); Brp=interp2(obj.zgrid,obj.rgrid,obj.Br',Zp,Rp,'makima'); Bp=sqrt(Bzp.^2+Brp.^2); % calculate the projection angle of the radial and axial % directions on the magnetic field line Costhet=Bzp./Bp; Sinthet=Brp./Bp; % calculate the actuale parallel velocity quantity=obj.species(1).VR(p,t,track).*Sinthet+obj.species(1).VZ(p,t,track).*Costhet; end function quantity=Vperp(obj,varargin) %Vperp Computes the perpendicular velocity in the guidind center reference frame, % for the main specie particle indices{1} at time indices{2} if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if strcmp(indices{1},':') p=1:obj.species(1).R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end % if track is true we look at specific particles with their % index and follow them in time % if it is false we just care about the distribution function % and specific particles can have different positions in the % resulting array for each timestep if size(indices,2)>2 track=indices{3}; else track=false; end % if gcs is true, gives the perpendicular velocity in the % guiding center system by substracting the EXB azimuthal % velocity % else gives the total perpendicular velocity if size(indices,2)>3 gcs=indices{4}; else gcs=false; end % get the particle position Zp=obj.species(1).Z(p,t,track); Rp=obj.species(1).R(p,t,track); % interpolate the magnetic field at the particle position Bzp=interp2(obj.zgrid,obj.rgrid,obj.Bz',Zp,Rp,'makima'); Brp=interp2(obj.zgrid,obj.rgrid,obj.Br',Zp,Rp,'makima'); Bp=sqrt(Bzp.^2+Brp.^2); % calculate the projecting angles Costhet=Bzp./Bp; Sinthet=Brp./Bp; Vdrift=zeros(size(Zp)); if gcs % for each particle and each timestep % calculate the azimuthal ExB drift velocity for j=1:length(t) [~, tfield]=min(abs(obj.t2d-obj.tpart(t(j)))); timeEr=obj.Er(:,:,tfield); timeEz=obj.Ez(:,:,tfield); %posindE=sub2ind(size(timeEr),Rind(:,j),Zind(:,j)); timeErp=interp2(obj.zgrid,obj.rgrid,timeEr,Zp(:,j),Rp(:,j)); timeEzp=interp2(obj.zgrid,obj.rgrid,timeEz,Zp(:,j),Rp(:,j)); Vdrift(:,j)=(timeEzp.*Brp(:,j)-timeErp.*Bzp(:,j))./Bp(:,j).^2; end end % calculate the perpendicular velocity quantity=sqrt((obj.species(1).VTHET(p,t,track)-Vdrift).^2+(obj.species(1).VR(p,t,track).*Costhet-obj.species(1).VZ(p,t,track).*Sinthet).^2); end function quantity=cyclphase(obj,varargin) %cyclphase Computes the cyclotronic phase for the main specie % for particles with indices{1} at time indices{2} if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if strcmp(indices{1},':') p=1:obj.species(1).R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end % if track is true we look at specific particles with their % index and follow them in time % if it is false we just care about the distribution function % and specific particles can have different positions in the % resulting array for each timestep if size(indices,2)>2 track=indices{3}; else track=false; end Zp=obj.species(1).Z(p,t,track); Rp=obj.species(1).R(p,t,track); % [~, zind(1)]=min(abs(obj.zgrid-0.005262)); % [~, zind(2)]=min(abs(obj.zgrid-0.006637)); % [~, rind(1)]=min(abs(obj.rgrid-0.0784)); % [~, rind(2)]=min(abs(obj.rgrid-0.07861)); % indices=Zp=obj.zgrid(zind(1)) &... % Rp=obj.rgrid(rind(1)); %Zp=Zp(indices); %Rp=Rp(indices); %p=indices; Bzp=interp2(obj.zgrid,obj.rgrid,obj.Bz',Zp,Rp,'makima'); Brp=interp2(obj.zgrid,obj.rgrid,obj.Br',Zp,Rp,'makima'); Bp=sqrt(Bzp.^2+Brp.^2); Costhet=Bzp./Bp; Sinthet=Brp./Bp; % compute the projection of the perpendicular velocity in the % radial direction vr=(obj.species(1).VR(p,t,track).*Costhet-obj.species(1).VZ(p,t,track).*Sinthet); % Get the perpendicular velocity vperp=obj.Vperp(p,t,track,true); vr=vr(indices); vperp=vperp(indices); cospsi=vr./vperp; quantity=acos(cospsi); end function p=borderpoints(obj,subdiv) %borderpoints Return a cell array containing the curves %defining the boundary of the domain % for each boundary p(1,:) and p(2,:) give axial and radial position % for each boundary p(3,:) and p(4,:) give axial and radial normals %gw= contourc(obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0]) p=cell(0,0); if nargin<2 subdiv=1; end ndiv=sum(subdiv); %outer cylinder if any(obj.geomweight(end,:,1)>=0) idp=ceil(length(obj.zgrid)/ndiv); imin=1; for j=1:length(subdiv) imax=min(imin+subdiv(j)*idp-1,length(obj.zgrid)); p{end+1}=[obj.zgrid(imin:imax)';obj.rgrid(end)*ones(imax-imin+1,1)'; zeros(imax-imin+1,1)';ones(imax-imin+1,1)']; imin=imax; end end %inner cylinder if any(obj.geomweight(1,:,1)>=0) idp=ceil(length(obj.zgrid)/ndiv); imin=1; for j=1:length(subdiv) imax=min(imin+subdiv(j)*idp-1,length(obj.zgrid)); p{end+1}=[obj.zgrid(imin:imax)';obj.rgrid(1)*ones(imax-imin+1,1)'; zeros(imax-imin+1,1)';-ones(imax-imin+1,1)']; imin=imax; end end if obj.walltype==2 % We have an elliptic insert that we want to isolate gw=obj.ellipseborder; zpos=obj.zgrid(obj.zgrid<(min(gw(1,:))) | obj.zgrid>(max(gw(1,:)))); p{2}=[zpos,obj.rgrid(end)*ones(size(zpos))]'; p{1}=[obj.zgrid';obj.rgrid(1)*ones(size(obj.zgrid))']; gw=obj.ellipseborder; p{3}=gw; elseif obj.walltype~=0 % extract all the walls gw=contourc(obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0]); [x,y,~]=C2xyz(gw); for i=1:length(x) %subdiv=[4,1,2]; ndiv=sum(subdiv); idp=ceil(length(x{i})/ndiv); imin=1; for j=1:length(subdiv) imax=min(imin+subdiv(j)*idp-1,length(x{i})); p{end+1}=[x{i}(imin:imax);y{i}(imin:imax)]; imin=imax; end end end % figure % for i=1:length(p) % plot(p{i}(1,:),p{i}(2,:)) % hold on % end end function p=ellipseborder(obj) %ellipseborder returns the boundary points defining the %elliptic insert z=linspace(-0.998,0.998,1000)*obj.z_r; p=zeros(4,length(z)); for i=1:length(z) p(1,i)=z(i)+obj.z_0; p(2,i)=obj.r_0-obj.r_r*sqrt(1-(z(i)/obj.z_r)^2); p(3,i)=2/(obj.z_r^2)*(z(i)); p(4,i)=2/(obj.r_r^2)*(p(2,i)-obj.r_0); end norm=sqrt(p(3,:).^2+p(4,:).^2); p(3,:)=double(obj.interior)*p(3,:)./norm; p(4,:)=double(obj.interior)*p(4,:)./norm; end function charge=totcharge(obj,fieldstep) % Integrates the density profile over the full volume to obtain % the total number of electrons in the volume n=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder,ones(size(obj.CellVol)), 1, 1); charge=sum(sum(n(:,:,fieldstep))); end function Gamma=Axialflux(obj,timestep,zpos) % Computes the axial particle flux n*Uz at timestep timestep and axial position zpos Gamma=obj.fluidUZ(:,zpos,timestep).*obj.N(:,zpos,timestep); end function Gamma=Metallicflux(obj,timestep,subdiv) % Computes the particle flux at time obj.t2d(timestep) on the % metallic boundaries if nargin<3 subdiv=1; end % We find the borderpoints p=obj.borderpoints(subdiv); gamma=cell(size(p)); Nr=cell(size(p)); Nz=cell(size(p)); for i=1:length(p) bp=p{i}; if size(bp,1)==2 % We get the normals at these positions and normalise them Nr{i}=-interp2(obj.zgrid,obj.rgrid,obj.geomweight(:,:,3),bp(1,:),bp(2,:)); Nz{i}=-interp2(obj.zgrid,obj.rgrid,obj.geomweight(:,:,2),bp(1,:),bp(2,:)); norm=sqrt(Nr{i}.^2+Nz{i}.^2); Nr{i}=Nr{i}./norm; Nz{i}=Nz{i}./norm; else Nr{i}=bp(4,:); Nz{i}=bp(3,:); end gamma{i}=zeros(size(bp,2),length(timestep)); end [z,r]=ndgrid(obj.zgrid,obj.rgrid); N=obj.N(:,:,timestep(1)); n=griddedInterpolant(z,r,N'); uz=griddedInterpolant(z,r,obj.fluidUZ(:,:,timestep(1))'); ur=griddedInterpolant(z,r,obj.fluidUR(:,:,timestep(1))'); % we get the density and fluid velocities at the desired time % steps and interpolate them at the boundary position for j=1:length(timestep) n.Values=obj.N(:,:,timestep(j))'; uz.Values=obj.fluidUZ(:,:,timestep(j))'; ur.Values=obj.fluidUR(:,:,timestep(j))'; for i=1:length(p) bp=p{i}; gamma{i}(:,j)=n(bp(1:2,:)').*(ur(bp(1:2,:)').*Nr{i}'+uz(bp(1:2,:)').*Nz{i}'); end end % return the boundary position p and the corresponding flux % gamma Gamma.p=p; Gamma.gamma=gamma; end function [pot] = PotentialWell(obj,fieldstep,fieldaligned) %PotentialWell Computes the potential well at the given timestep on the FEM grid points % interpolates the model data on rgrid and zgrid if nargin<3 - fieldaligned=false; + fieldaligned=false; end model=obj.potentialwellmodel(fieldstep); z=model.z; modpot=model.pot; - + if fieldaligned r=model.rathet; lvls=linspace(min(obj.rAthet(:)),max(obj.rAthet(:)),400); [Zmesh,Rmesh]=meshgrid(obj.zgrid,lvls); else r=model.r; [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid); end pot=zeros(size(Zmesh,2),size(Zmesh,1),length(fieldstep)); for i=1:length(fieldstep) pot(:,:,i)=griddata(z,r,modpot(:,i),Zmesh,Rmesh)'; end end function Epar = Epar(obj,fieldstep) % Computes the electric field component parallel to the magnetic field line Epar=obj.Er(:,:,fieldstep).*(obj.Br./obj.B)' + (obj.Bz./obj.B)'.*obj.Ez(:,:,fieldstep); end function Eperp = Eperp(obj,fieldstep) % Computes the electric field component perpendicular to the magnetic field line Eperp=obj.Er(:,:,fieldstep).*(obj.Bz./obj.B)' - (obj.Br./obj.B)'.*obj.Ez(:,:,fieldstep); end function Ekin = Ekin(obj,varargin) %Ekin Computes the classical kinetic energy of particles indices{1} at % time obj.tpart(indices{2}) in Joules if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if strcmp(indices{1},':') p=1:obj.species(1).R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end % if track is true we look at specific particles with their % index and follow them in time % if it is false we just care about the distribution function % and specific particles can have different positions in the % resulting array for each timestep if size(indices,1)>2 track=indices{3}; else track=false; end Vr=obj.species(1).VR(p,t,track); Vthet= obj.species(1).VTHET(p,t,track); Vz=obj.species(1).VZ(p,t,track); Ekin=0.5*obj.msim/obj.weight*(Vr.^2+Vthet.^2+Vz.^2); end function sig=sigio(obj,E,init) %sigio returns the total ionisation cross-section in m^2 % at energy E[eV] % init is only used during the loading of the h5 file if nargin <3 init=false; end sig=zeros(size(E)); if(~init &&( ~obj.neutcol.present || isempty(obj.neutcol.io_cross_sec))) sig=zeros(size(E)); return end for i=1:length(E(:)) if(E(i)>obj.neutcol.Eion) sig(ind2sub(size(E),i))=obj.fit_cross_sec(E(ind2sub(size(E),i)),obj.neutcol.io_cross_sec); end end end function sig=sigmio(obj,E) %sigmio returns the total ionisation cross-section for momentum exchange for the incoming electron in m^2 % at energy E[eV] sig=zeros(size(E)); if(~obj.neutcol.present || isempty(obj.neutcol.iom_cross_sec)) return end for i=1:length(E(:)) if(E(i)>obj.neutcol.Eion) sig(ind2sub(size(E),i))=obj.fit_cross_sec(E(ind2sub(size(E),i)),obj.neutcol.iom_cross_sec); end end end function sigm=sigmela(obj,E) %sigmela returns the elastic collision cross-section for momentum exchange for the incoming electron in m^2 % at energy E[eV] sigm=zeros(size(E)); if(~obj.neutcol.present || isempty(obj.neutcol.ela_cross_sec)) return end for i=1:length(E(:)) sigm(ind2sub(size(E),i))=obj.fit_cross_sec(E(ind2sub(size(E),i)),obj.neutcol.ela_cross_sec); end end function sig=sigela(obj,E) %sigmela returns the elastic collision cross-section for the incoming electron in m^2 % at energy E[eV] % if used this will give the frequency of elastic collisions E0=obj.neutcol.E0; chi=E./(0.25*E0+E); sig=(2*chi.^2)./((1-chi).*((1+chi).*log((1+chi)./(1-chi))-2*chi)).*obj.sigmela(E); end function [Forces, Density]=Forcespline(obj,it,fdens,getmean) %Forcesplie calculates the fluid force terms in each direction %at time obj.t2d(it) % if fdens return the force density in N/m^3 othewise give % the force in N % if getmean return only the time averaged quanties over % time samples[it(1)...it(end] if strcmp(it,':') it=floor(0.95*size(obj.t2d)):size(obj.t2d)-1; end if nargin<3 fdens=true; end if nargin <4 getmean=false; end % To be able to calculate the centered finite difference in % time, we remove the first and last time indices it(it<2)=[]; it(it>length(obj.t2d)-1)=[]; m_e=obj.msim/obj.weight; q_e=obj.qsim/obj.weight; n=obj.N(:,:,it); [r,~]=meshgrid(obj.rgrid,obj.zgrid); Rinv=1./r'; Rinv(isinf(Rinv))=0; % get inverse of density to get the force in N Density.N=n; invn=1./n; invn(isnan(invn) | isinf(invn))=0; % Calculate electric forces Eforcer=q_e*obj.Er(:,:,it); Eforcez=q_e*obj.Ez(:,:,it); Dragforcer=zeros(size(n,1),size(n,2),size(n,3)); Dragforcethet=zeros(size(n,1),size(n,2),size(n,3)); Dragforcez=zeros(size(n,1),size(n,2),size(n,3)); time=obj.t2d(it); Forces.it=it; Forces.time=time; if getmean if ~fdens n=ones(size(n)); end % Electric forces Forces.Eforcer=mean(n.*q_e.*obj.Er(:,:,it),3); Forces.Eforcez=mean(n.*q_e.*obj.Ez(:,:,it),3); % Magnetic forces Forces.Bforcer=mean(q_e.*obj.fluidUTHET(:,:,it).*obj.Bz'.*n,3); Forces.Bforcethet=mean(q_e.*(obj.fluidUZ(:,:,it).*obj.Br'-obj.fluidUR(:,:,it).*obj.Bz').*n,3); Forces.Bforcez=mean(-q_e.*obj.fluidUTHET(:,:,it).*obj.Br'.*n,3); % Inertial forces Forces.inertforcer=mean(-m_e.*n.*(-obj.fluidUTHET(:,:,it).^2.*Rinv... +obj.fluidUR(:,:,it).*obj.fluidUR.der(:,:,it,[1 0])... +obj.fluidUZ(:,:,it).*obj.fluidUR.der(:,:,it,[0 1])),3); Forces.inertforcethet=mean(-m_e*n.*(obj.fluidUR(:,:,it).*obj.fluidUTHET(:,:,it).*Rinv... +obj.fluidUR(:,:,it).*obj.fluidUTHET.der(:,:,it,[1 0])... +obj.fluidUZ(:,:,it).*obj.fluidUTHET.der(:,:,it,[0 1])),3); Forces.inertforcez=mean(-m_e*n.*(obj.fluidUR(:,:,it).*obj.fluidUZ.der(:,:,it,[1 0])... +obj.fluidUZ(:,:,it).*obj.fluidUZ.der(:,:,it,[0 1])),3); % Pressure forces Forces.Pressforcer=mean(-n.*( squeeze(obj.Presstens.der(1,:,:,it,[1 0]))... + squeeze(obj.Presstens(1,:,:,it) - obj.Presstens(4,:,:,it)).*Rinv... + squeeze(obj.Presstens.der(3,:,:,it,[0 1])))... .*invn,3); Forces.Pressforcethet=mean(-n.*( squeeze(obj.Presstens.der(2,:,:,it,[1 0]))... + squeeze(obj.Presstens.der(5,:,:,it,[0 1])) ... + 2*squeeze(obj.Presstens(2,:,:,it)).*Rinv ... ).*invn,3); Forces.Pressforcez=mean(-n.*( squeeze(obj.Presstens.der(3,:,:,it,[1 0]))... + squeeze(obj.Presstens(3,:,:,it)).*Rinv... + squeeze(obj.Presstens.der(6,:,:,it,[0 1])) )... .*invn,3); % ellastic coll drag forces if( obj.neutcol.present) Ek=squeeze(obj.fluidEkin(1,:,:,it)+obj.fluidEkin(2,:,:,it)+obj.fluidEkin(3,:,:,it)); sigm=obj.sigmela(Ek/obj.qe)+obj.sigio(Ek/obj.qe)+obj.sigmio(Ek/obj.qe); dragfreq=obj.neutcol.neutdens.*sigm.*sqrt(2*obj.weight/obj.msim*Ek); Forces.Dragforcer=mean(-m_e*n.*dragfreq.*obj.fluidUR(:,:,it),3); Forces.Dragforcethet=mean(-m_e*n.*dragfreq.*obj.fluidUTHET(:,:,it),3); Forces.Dragforcez=mean(-m_e*n.*dragfreq.*obj.fluidUZ(:,:,it),3); else Forces.Dragforcer=0; Forces.Dragforcethet=0; Forces.Dragforcez=0; end % effective drag frequency due to the maxwellian source if( obj.maxwellsrce.present) dragfreqsrc=obj.maxwellsrce.frequency*obj.weight/(pi*diff(obj.maxwellsrce.zlim)*(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2)).*invn; dragfreqsrc(isinf(dragfreqsrc))=0; Forces.Dragforcer=Forces.Dragforcer+mean(-n.*m_e.*dragfreqsrc.*obj.fluidUR(:,:,it),3); Forces.Dragforcethet=Forces.Dragforcethet+mean(-n.*m_e.*dragfreqsrc.*obj.fluidUTHET(:,:,it),3); Forces.Dragforcez=Forces.Dragforcez+mean(-n.*m_e.*dragfreqsrc.*obj.fluidUZ(:,:,it),3); end % Time derivative for fluid accelleration cdt=(obj.t2d(it+1)-obj.t2d(it-1)); cdt=reshape(cdt,1,1,[]); Forces.durdt=mean(m_e*(obj.fluidUR(:,:,it+1)-obj.fluidUR(:,:,it-1))./cdt,3); Forces.duthetdt=mean(m_e*(obj.fluidUTHET(:,:,it+1)-obj.fluidUTHET(:,:,it-1))./cdt,3); Forces.duzdt=mean(m_e*(obj.fluidUZ(:,:,it+1)-obj.fluidUZ(:,:,it-1))./cdt,3); else % Allocate memory Bforcer=zeros(size(n,1),size(n,2),size(n,3)); Bforcez=zeros(size(n,1),size(n,2),size(n,3)); Bforcethet=zeros(size(n,1),size(n,2),size(n,3)); inertforcer=zeros(size(n,1),size(n,2),size(n,3)); inertforcez=zeros(size(n,1),size(n,2),size(n,3)); inertforcethet=zeros(size(n,1),size(n,2),size(n,3)); Pressforcer=zeros(size(n,1),size(n,2),size(n,3)); Pressforcethet=zeros(size(n,1),size(n,2),size(n,3)); Pressforcez=zeros(size(n,1),size(n,2),size(n,3)); durdt=zeros(size(n,1),size(n,2),size(n,3)); duthetdt=zeros(size(n,1),size(n,2),size(n,3)); duzdt=zeros(size(n,1),size(n,2),size(n,3)); fluiduThet=obj.fluidUTHET(:,:,it); Density.fluiduThet=fluiduThet; for j=1:size(n,3) % Magnetic forces Bforcer(:,:,j)=q_e.*fluiduThet(:,:,j).*obj.Bz'; Bforcethet(:,:,j)=q_e.*(obj.fluidUZ(:,:,it(j)).*obj.Br'-obj.fluidUR(:,:,it(j)).*obj.Bz'); Bforcez(:,:,j)=-q_e.*fluiduThet(:,:,j).*obj.Br'; % Inertial forces inertforcer(:,:,j)=-m_e.*(-fluiduThet(:,:,j).^2.*Rinv... +obj.fluidUR(:,:,it(j)).*obj.fluidUR.der(:,:,it(j),[1 0])... +obj.fluidUZ(:,:,it(j)).*obj.fluidUR.der(:,:,it(j),[0 1])); inert1=obj.fluidUR(:,:,it(j)).*fluiduThet(:,:,j).*Rinv; inert2=obj.fluidUR(:,:,it(j)).*obj.fluidUTHET.der(:,:,it(j),[1 0]); inert3=obj.fluidUZ(:,:,it(j)).*obj.fluidUTHET.der(:,:,it(j),[0 1]); inertforcethet(:,:,j)=-m_e.*(inert1... +inert2... +inert3); inertforcez(:,:,j)=-m_e.*(obj.fluidUR(:,:,it(j)).*obj.fluidUZ.der(:,:,it(j),[1 0])... +obj.fluidUZ(:,:,it(j)).*obj.fluidUZ.der(:,:,it(j),[0 1])); % Pressure forces Pr1=squeeze(obj.Presstens.der(1,:,:,it(j),[1 0])); Pr2=squeeze(obj.Presstens(1,:,:,it(j)) - obj.Presstens(4,:,:,it(j))).*Rinv; Pr3=squeeze(obj.Presstens.der(3,:,:,it(j),[0 1])); Pressforcer(:,:,j)=-( Pr1... + Pr2... + Pr3 )... .*invn(:,:,j); Pthet1=squeeze(obj.Presstens.der(2,:,:,it(j),[1 0])); Pthet2=squeeze(obj.Presstens.der(5,:,:,it(j),[0 1])); Pthet3=2*squeeze(obj.Presstens(2,:,:,it(j))).*Rinv; Pressforcethet(:,:,j)=-( Pthet1... + Pthet2 ... + Pthet3 ... ).*invn(:,:,j); Pz1=squeeze(obj.Presstens.der(3,:,:,it(j),[1 0])); Pz2=squeeze(obj.Presstens(3,:,:,it(j))).*Rinv; Pz3=squeeze(obj.Presstens.der(6,:,:,it(j),[0 1])); Pressforcez(:,:,j)=-( Pz1... + Pz2... + Pz3 )... .*invn(:,:,j); % ellastic coll drag forces if( obj.neutcol.present) Ek=squeeze(obj.fluidEkin(1,:,:,it(j))+obj.fluidEkin(2,:,:,it(j))+obj.fluidEkin(3,:,:,it(j))); sigm=obj.sigmela(Ek/obj.qe)+obj.sigio(Ek/obj.qe)+obj.sigmio(Ek/obj.qe); dragfreq=obj.neutcol.neutdens.*sigm.*sqrt(2*obj.weight/obj.msim*Ek); Dragforcer(:,:,j)=-m_e*dragfreq.*obj.fluidUR(:,:,it(j)); Dragforcethet(:,:,j)=-m_e*dragfreq.*obj.fluidUTHET(:,:,it(j)); Dragforcez(:,:,j)=-m_e*dragfreq.*obj.fluidUZ(:,:,it(j)); end % effective drag frequency due to the maxwellian source if( obj.maxwellsrce.present) dragfreqsrc=obj.maxwellsrce.frequency*obj.weight/(pi*diff(obj.maxwellsrce.zlim)*(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2))*invn(:,:,j); dragfreqsrc(isinf(dragfreqsrc))=0; Dragforcer(:,:,j)=Dragforcer(:,:,j)+-m_e*dragfreqsrc.*obj.fluidUR(:,:,it(j)); Dragforcethet(:,:,j)=Dragforcethet(:,:,j)+-m_e*dragfreqsrc.*obj.fluidUTHET(:,:,it(j)); Dragforcez(:,:,j)=Dragforcez(:,:,j)+-m_e*dragfreqsrc.*obj.fluidUZ(:,:,it(j)); end % Time derivative cdt=(obj.t2d(it(j)+1)-obj.t2d(it(j)-1)); durdt(:,:,j)=m_e*(obj.fluidUR(:,:,it(j)+1)-obj.fluidUR(:,:,it(j)-1))/cdt; duthetdt(:,:,j)=m_e*(obj.fluidUTHET(:,:,it(j)+1)-obj.fluidUTHET(:,:,it(j)-1))/cdt; duzdt(:,:,j)=m_e*(obj.fluidUZ(:,:,it(j)+1)-obj.fluidUZ(:,:,it(j)-1))/cdt; end if(~fdens) Forces.Eforcer=Eforcer; Forces.Eforcez=Eforcez; Forces.Bforcer=Bforcer; Forces.Bforcethet=Bforcethet; Forces.Bforcez=Bforcez; Forces.inertforcer=inertforcer; Forces.inertforcethet=inertforcethet; Forces.inertforcez=inertforcez; Forces.Pressforcer=Pressforcer; Forces.Pressforcethet=Pressforcethet; Forces.Pressforcez=Pressforcez; Forces.durdt=durdt; Forces.duthetdt=duthetdt; Forces.duzdt=duzdt; Forces.Dragforcer=Dragforcer; Forces.Dragforcethet=Dragforcethet; Forces.Dragforcez=Dragforcez; else % multiply by density to have force density Forces.Eforcer=Eforcer.*n; Forces.Eforcez=Eforcez.*n; Forces.Bforcer=Bforcer.*n; Forces.Bforcethet=Bforcethet.*n; Forces.Bforcez=Bforcez.*n; Forces.inertforcer=inertforcer.*n; Forces.inertforcethet=inertforcethet.*n; Forces.inertforcez=inertforcez.*n; Forces.Pressforcer=Pressforcer.*n; Forces.Pressforcethet=Pressforcethet.*n; Forces.Pressforcez=Pressforcez.*n; Forces.durdt=durdt.*n; Forces.duthetdt=duthetdt.*n; Forces.duzdt=duzdt.*n; Forces.Dragforcer=Dragforcer.*n; Forces.Dragforcethet=Dragforcethet.*n; Forces.Dragforcez=Dragforcez.*n; end end end function [lr,rb,lz,zb]= clouddims(obj,it,zpos,fracn) % clouddims return the cloud axial and radial limit at time it % and axial position zpos % fracn defines the fraction of the maximum density below which % we consider to have a vacuum if nargin<4 fracn=0.1; end % get the density n=obj.N(:,:,it); lr=cell(1,length(it)); lz=lr; rb=lr; zb=rb; for i=1:size(n,3) nthresh=fracn*max(max(n(:,:,i))); % find the points outside of the cloud outside=find(n(:,zpos,i)2) rmpos=outside(j); rppos=outside(j+1); lr{i}(k)=obj.rgrid(rppos-1)-obj.rgrid(rmpos+1); rb{i}(:,k)=[max(rmpos+1,1) min(rppos-1,sum(obj.nnr))]; k=k+1; end end maxgap=2; k=1; for I=rmpos+1:rppos-1 outside=find(n(I,:,i)maxgap) maxgap=zgap(j); zmpos=outside(j); zppos=outside(j+1); lz{i}(k)=obj.zgrid(zppos-1)-obj.zgrid(zmpos+1); zb{i}(:,k)=[max(zmpos+1,1) min(zppos-1,obj.nz)]; k=k+1; end end end end end %------------------------------------------ % Functions for plotting evolving quantities - function displaysplbound(obj,ax,rescale) + function line=displaysplbound(obj,ax,rescale,markers) %displaysplbound display on axis ax the boundary of the %simulation domain and the Dirichlet and Neumann walls defined %with spline curves if nargin<2 ax=gca; end if nargin<3 rescale=1; end + if nargin<4 + markers=true; + end hold on for i=1:obj.spl_bound.nbsplines knots=obj.spl_bound.boundary(i).knots(1:end); coeffs=obj.spl_bound.boundary(i).coefs'*rescale; pp=spmak(knots,coeffs); sizec=size(coeffs,2); order=length(knots)-sizec; s=linspace(knots(order),knots(sizec+1),1000); fittedpos=fnval(pp,s); - plot(fittedpos(1,:),fittedpos(2,:),'-') - plot(coeffs(1,:),coeffs(2,:),'rx','markersize',14) + line=plot(fittedpos(1,:),fittedpos(2,:),'-','linewidth',2); + if markers + plot(coeffs(1,:),coeffs(2,:),'rx','markersize',14) + end end end function displayraddim(obj,it,zpos,fracn) %displayraddim display the evolution of the radial dimension of the cloud in %time to find if the cloud size get below a critical radial %size at which the ionisation is not sufficient to compensate %the losses % also plot the well radial dimensions in time if nargin<3 zpos=floor(length(obj.zgrid)/2); end if nargin<4 fracn=0.1; end [lr,rb,lz,zb]=obj.clouddims(it,zpos,fracn); t=obj.t2d(it); Lr=zeros(size(lr)); er=obj.Er(:,:,it); r_min=Lr; r_minpred=r_min; well_r=Lr; nb=Lr; for i=1:length(lr) if ~isempty(lr{i}) && ~isempty(lz{i}) [Lr(i),id]=max(lr{i}); rm=rb{i}(1,id); rp=rb{i}(2,id); nb(i)=mean(obj.N(rm:rp,zpos,it(i))); Lp=min(lz{i}); Lm=mean(lz{i}); rpos=rm:rp; vperp=-er(rpos,zpos,i)./obj.Bz(zpos,rpos)'; Ek=0.5*obj.me*vperp.^2/obj.qe; sigio=obj.sigio(Ek); sigd=obj.sigmela(Ek)+obj.sigmio(Ek)+sigio; omegap2=obj.qe^2*obj.N(rpos,zpos,it(i))/obj.eps_0/obj.me; omegac2=(obj.qe*obj.Bz(zpos,rpos)'/obj.me).^2; ur=er(rpos,zpos,i)*obj.qe./((omegap2-omegac2)*obj.me).*sigd.*vperp*obj.neutcol.neutdens; r_minpred(i)=mean(obj.N(rp,zpos,it(i))*Lp*ur./(nb(i)*obj.neutcol.neutdens*sigio.*vperp*Lm));%mean(1./(-1/obj.rgrid(rm)+obj.neutcol.neutdens*sigio.*vperp./ur*(Lm/Lp)*nb(i)/obj.N(rp,zpos,it(i)))); rpos=rp; vperp=-er(rpos,zpos,i)./obj.Bz(zpos,rpos)'; Ek=0.5*obj.me*vperp.^2/obj.qe; sigio=obj.sigio(Ek); ur=obj.fluidUR(rpos,zpos,it(i)); r_min(i)=max(obj.N(rp,zpos,it(i))*Lp*ur/(nb(i)*obj.neutcol.neutdens*sigio.*vperp*Lm),0);%max(mean(1./(-1/obj.rgrid(rm)+obj.neutcol.neutdens*sigio.*vperp./ur*(Lm/Lp)*nb(i)/obj.N(rp,zpos,it(i)))),0); nb(i)=nb(i)*Lm*2*pi*obj.rgrid(rm)*Lr(i); else Lr(i)=NaN; r_min(i)=NaN; r_minpred(i)=NaN; end potwell=obj.PotentialWell(it(i))'; outside=find(isnan(potwell(:,zpos))); gap=diff(outside); for j=1:length(gap) if(gap(j)>2) rmpos=outside(j)+1; rppos=outside(j+1)-1; well_r(i)=obj.rgrid(rppos)-obj.rgrid(rmpos); end end end f=figure('Name', sprintf('%s rlims B=%f phi=%f',obj.name,obj.B0, (obj.potout-obj.potinn))); plot(t,Lr,'displayname','\Deltar_{cloud}','linewidth',1.3) hold on plot(t,r_min,'displayname','\Deltar_{min} (u_r simu)','linewidth',1.3) plot(t,r_minpred,'displayname','\Deltar_{min} (u_r pred)','linewidth',1.3) plot(t,well_r,'displayname','\Deltar_{well}','linewidth',1.3) ylabel('\Delta r [m]') yyaxis right plot(t,nb,'--','displayname','N') legend('location','eastoutside') xlabel('t [s]') ylabel('N') set(gca,'fontsize',12) yyaxis left ylimits=ylim; %ylim([ylimits(1) 1.1*max(Lr)]) title(sprintf('cloud radial limits at z=%1.2e[m]',obj.zgrid(zpos))) obj.savegraph(f,sprintf('%s/%s_%d_rlims',obj.folder,obj.name,zpos),[15 10]); end + function displayDensPsi(obj,tsteps) + %DISPLAYDENSPSI plot the fluid density averaged in time for times t2d(tsteps) + % Displays also the geometry of the metallic boundaries and the + % magnetic field lines + % Show on top the contour of the Psi function for the Davidson + % loading + f=figure('Name', sprintf('%sfluid_dens_psi',obj.name)); + ax1=gca; + if nargin<2 + tsteps=length(obj.t2d); + end + + geomw=obj.geomweight(:,:,1); + dens=mean(obj.N(:,:,tsteps),3); + dens(geomw<=0)=0; + maxdens=max(dens(:)); + [h,curve]=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,dens,30); + hold on + set(curve,'linecolor','none'); + %contour(ax1,obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0],'r-','linewidth',1.5,'Displayname','Boundaries'); + if(obj.conformgeom) + ylim(ax1,[obj.rgrid(1) obj.rgrid(sum(obj.nnr(1:2)))]*1000) + else + ylim(ax1,[obj.rgrid(1) obj.rgrid(end)]*1000) + end + xlim(ax1,[obj.zgrid(1) obj.zgrid(end)]*1000) + xlabel(ax1,'z [mm]') + ylabel(ax1,'r [mm]') + %title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',obj.t2d(fieldstart),obj.t2d(fieldend),double(maxdens))) + c = colorbar(ax1); + c.Label.String= 'electron density [m^{-3}]'; + view(ax1,2) + colormap(flipud(hot)); + + %% Magnetic field lines + Blines=obj.rAthet; + levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),15); + [~,h1]=contour(ax1,obj.zgrid*1000,obj.rgrid*1000,Blines,real(levels),'-.','color','k','linewidth',1.5,'Displayname','Magnetic field lines'); + + + %% Equipotential lines + Pot=mean(obj.pot(:,:,tsteps),3); + Pot(obj.geomweight(:,:,1)<0)=NaN; + %levels=8;%[-3.4 -5 -10 -15 -20 -25];%7; + potcolor='b';%[0.3660 0.6740 0.1880]; + [c1,h2]=contour(ax1,obj.zgrid*1000,obj.rgrid*1000,Pot,'--','color',potcolor,'ShowText','on','linewidth',2,'Displayname','Equipotentials [kV]'); + clabel(c1,h2,'Color',potcolor) + clabel(c1,h2, 'labelspacing', 100); + clabel(c1,h2, 'FontSize', 12); + + % Draw the metallic boundaries and the geometry itself + [c1,hContour]=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,-geomw,[0,0],'linewidth',1.5); + drawnow; + xlim(ax1,[obj.zgrid(1)*1000 obj.zgrid(end)*1000]) + % Change the color of the metallic boundaries to grey + hFills=hContour.FacePrims; + [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' + try + hFills(end).ColorData = uint8([150;150;150;255]); + for idx = 1 : numel(hFills)-1 + hFills(idx).ColorData(4) = 0; % default=255 + end + catch + end + + h0=obj.H0; + p0=obj.P0; + + if(obj.species(1).R.nt<2) + h0=obj.H0; + p0=obj.P0; + else + h0=mean(H(obj,{1:obj.species(1).VR.nparts,obj.species(1).VR.nt,false})); + p0=mean(P(obj,{1:obj.species(1).VR.nparts,obj.species(1).VR.nt,false})); + end + + + Mirrorratio=(obj.Rcurv-1)/(obj.Rcurv+1); + locpot=mean(obj.pot(:,:,end),3); + psi=1+obj.qe*locpot/h0-1/(2*obj.me*h0)*(p0./obj.rgrid+obj.qe*0.5*obj.B0.*(obj.rgrid-obj.width/pi*Mirrorratio*cos(2*pi*obj.zgrid'/obj.width).*besseli(1,2*pi*obj.rgrid/obj.width))).^2; + [~,hpsi]=contour(obj.zgrid*1000,obj.rgrid*1000,psi,[0 0],'m:','displayname','\Psi','linewidth',2); + + grid on + + legend([h1,h2,hpsi],{'Magnetic field lines','Equipotentials [V]','\Psi(r,z)=0'},'location','northeast') + + f.PaperOrientation='landscape'; + f.PaperUnits='centimeters'; + papsize=[16 14]; + f.PaperSize=papsize; + set(ax1,'fontsize',14) + %axis equal + obj.savegraph(f,sprintf('%sfluid_dens_psi',obj.name)) + end + + + + function displaypsi(obj,deltat) %% plot the initial and final radial profile at position z=0 and show the normalized enveloppe function Psi % relevant for Davidson annular distribution function f=figure('Name', sprintf('%s Psi',obj.name)); f.Name= sprintf('%s Psi',obj.name); zpos=floor(length(obj.zgrid)/2); tinit=1; tend=length(obj.t2d); if iscell(deltat) deltat=cell2mat(deltat); end if(obj.species(1).R.nt<2) h0=obj.H0; p0=obj.P0; else h0=mean(H(obj,{1:obj.species(1).VR.nparts,obj.species(1).VR.nt,false})); p0=mean(P(obj,{1:obj.species(1).VR.nparts,obj.species(1).VR.nt,false})); end lw=1.5; Mirrorratio=(obj.Rcurv-1)/(obj.Rcurv+1); locpot=mean(obj.pot(:,zpos,tend-deltat:tend),3); psi=1+obj.qe*locpot(:)/h0-1/(2*obj.me*h0)*(p0./obj.rgrid+obj.qe*0.5*obj.B0.*(obj.rgrid-obj.width/pi*Mirrorratio*cos(2*pi*obj.zgrid(zpos)/obj.width)*besseli(1,2*pi*obj.rgrid/obj.width))).^2; + + locdens=mean(obj.N(:,zpos,tend-deltat:tend),3); [maxn,In]=max(locdens);%M.N(:,zpos,tinit)); - plot(obj.rgrid,obj.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw) + + plot(obj.rgrid*1e3,obj.N(:,zpos,tinit),'b-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw) hold on - plot(obj.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(tend-deltat)*1e9,obj.t2d(tend)*1e9),'linewidth',lw) - plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end)*maxn*obj.rgrid(In),'DisplayName','N=c*1/r','linewidth',lw) - plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end).^2*maxn*obj.rgrid(In)^2,'DisplayName','N=c*1/r^2','linewidth',lw) - plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end).^4*maxn*obj.rgrid(In)^4,'DisplayName','N=c*1/r^4','linewidth',lw) - xlabel('r [m]') + plot(obj.rgrid*1e3,locdens,'r--','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(tend-deltat)*1e9,obj.t2d(tend)*1e9),'linewidth',lw) + %r0=0.005;% + r0=obj.rgrid(In); + %maxn=5e14; + plot(obj.rgrid(In-5:end)*1e3,1./obj.rgrid(In-5:end)*r0*maxn,'DisplayName','n_{e,fit}=c/r','linewidth',lw) + %plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end).^2*maxn*obj.rgrid(In)^2,'DisplayName','N=c*1/r^2','linewidth',lw) + %plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end).^4*maxn*obj.rgrid(In)^4,'DisplayName','N=c*1/r^4','linewidth',lw) + xlabel('r [mm]') ylabel('n_e [m^{-3}]') I=find(psi>0); if (length(I)>1) I=[I(1)-2; I(1)-1; I; I(end)+1; I(end)+2]; else I=obj.nnr(1):length(psi); end rq=linspace(obj.rgrid(max(I(1),1)),obj.rgrid(I(end)),500); psiinterp=interp1(obj.rgrid(I),psi(I),rq,'pchip'); zeroindices=find(diff(psiinterp>=0),2); maxpsiinterp=max(psiinterp); - plot(rq,maxn*psiinterp/abs(maxpsiinterp),'Displayname','normalized \Psi [a.u.]','linewidth',lw) - ylim([0 inf]) + plot(rq*1e3,maxn*psiinterp/abs(maxpsiinterp),'Displayname','normalized \Psi [a.u.]','linewidth',lw) + ylim([0 inf]); + drawnow + ylimits=ylim; for i=1:length(zeroindices) - border=plot([rq(zeroindices(i)) rq(zeroindices(i))],[0 obj.rgrid(In)/obj.rgrid(In-2)*maxn],'k--','linewidth',lw); + border=plot([rq(zeroindices(i)) rq(zeroindices(i))]*1e3,[0 1./obj.rgrid(In-5)*r0*maxn],'k--','linewidth',lw); set(get(get(border,'Annotation'),'LegendInformation'),'IconDisplayStyle','off'); end legend - xlim([0 0.02]) + xlim([0 0.02]*1e3) grid on - title(sprintf('Radial density profile at z=%1.2e[m]',obj.zgrid(zpos))) + %title(sprintf('Radial density profile at z=%1.2e[m]',obj.zgrid(zpos))) + set(gca,'fontsize',12) obj.savegraph(f,sprintf('%sPsi',obj.name),[15 10]); end function f=displayrprofile(obj,t,zpos,init) %% plot the initial and final radial profile at the axial center of the simulation space % also plot the azimuthal fluid rotation frequency profile % t: time index considered % zpos: axial position index % init: initial time considered for comparison f=figure('Name', sprintf('%s Prof',obj.name)); if nargin < 3 || length(zpos)<1 zpos=floor(length(obj.zgrid)/2); end if nargin<4 init=false; end if(iscell(t)) t=cell2mat(t); end lw=1.5; if init tinit=t(1); t=t(2:end); end locdens=mean(obj.N(:,zpos,t),3); %inverse of radius Rinv=1./obj.rgrid; Rinv(isnan(Rinv))=0; %azimuthal velocity and azimuthal rotation frequency in m/s and %1/s vthet=mean(obj.fluidUTHET(:,zpos,t),3); omegare=(vthet.*Rinv); % plot the initial density if(init) plot(obj.rgrid,obj.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw) end hold on %plot the time averaged current density plot(obj.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(t(1))*1e9,obj.t2d(t(end))*1e9),'linewidth',lw) xlabel('r [m]') ylabel('n_e [m^{-3}]') legend('location','Northwest') % limit the axis to the simulation domain if obj.conformgeom xlim([obj.rgrid(1) obj.rgrid(sum(obj.nnr(1:2)))]) else xlim([obj.rgrid(1) obj.rgrid(end)]) end grid on ylimits=ylim(); % plot the metallic walls for a constant radius coaxial % configuration if obj.conformgeom plot(obj.rgrid(1)*[1 1],ylimits,'k--') plot(obj.rgrid(end)*[1 1],ylimits,'k--') else plot(obj.r_a*[1 1],ylimits,'k--') if obj.walltype==0 plot(obj.r_b*[1 1],ylimits,'k--') elseif obj.walltype==1 rmax=obj.r_0-obj.r_r*sqrt(1-(obj.zgrid(zpos)-obj.z_0)^2/obj.z_r^2); plot(rmax*[1 1],ylimits,'k--') end end yyaxis right % plot the azimuthal fluid rotation frequency profile plot(obj.rgrid,omegare,'DisplayName',sprintf('<\\omega_{re}> t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(t(1))*1e9,obj.t2d(t(end))*1e9),'linewidth',lw) ylabel('\omega_{re} [1/s]') title(sprintf('Radial density profile at z=%1.2e[m]',obj.zgrid(zpos))) obj.savegraph(f,sprintf('%srProf',obj.name),[15 10]); end function displayenergy(obj) %% Plot the time evolution of the system energy and number of simulated macro particles - tmin=1; + tmin=2; tmax=length(obj.ekin); f=figure('Name', sprintf('%s Energy',obj.name)); - subplot(2,1,1) - plot(obj.t0d(tmin:tmax),obj.ekin(tmin:tmax),'o-',... - obj.t0d(tmin:tmax),obj.epot(tmin:tmax),'d-',... - obj.t0d(tmin:tmax),obj.etot(tmin:tmax),'h-',... - obj.t0d(tmin:tmax),obj.etot0(tmin:tmax),'h-',... - obj.t0d(tmin:tmax),obj.eerr(tmin:tmax),'x--') + subplot(1,2,1) + plot(obj.t0d(tmin:tmax)*1e9,obj.ekin(tmin:tmax),'-',... + obj.t0d(tmin:tmax)*1e9,obj.epot(tmin:tmax),':',... + obj.t0d(tmin:tmax)*1e9,obj.etot(tmin:tmax),'-.',... + obj.t0d(tmin:tmax)*1e9,obj.etot0(tmin:tmax),'--','linewidth',4) %obj.t0d(tmin:tmax),obj.ekin(tmin:tmax)-obj.epot(tmin:tmax),'--', - legend('E_{kin}', 'E_{pot}', 'E_{tot}','E_{ref}','E_{err}') - xlabel('Time [s]') + legend('E_{kin}', 'E_{pot}', 'E_{tot}','E_{ref}','location','east') + xlabel('Time [ns]') ylabel('Energies [J]') grid on xlimits=xlim(); - - subplot(2,1,2) + set(gca,'fontsize',14) + subplot(1,2,2) try - semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)./obj.etot0(tmin:tmax)),'h-') + semilogy(obj.t0d(tmin:tmax)*1e9,abs(obj.eerr(tmin:tmax)./obj.etot0(tmin:tmax)),'-','linewidth',1.5) catch - semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)/obj.etot(2)),'h-') + semilogy(obj.t0d(tmin:tmax)*1e9,abs(obj.eerr(tmin:tmax)/obj.etot(2)),'-','linewidth',1.5) end hold on - xlabel('t [s]') - ylabel('E_{err}/E_{tot}') + xlabel('Time [ns]') + ylabel('abs(E_{tot}-E_{ref})/E_{tot}') xlim(xlimits) grid on - try - yyaxis right - plot(obj.t0d(tmin:tmax),abs(obj.npart(tmin:tmax)./obj.npart(1)*100),'d--') - ylabel('Nparts %') - %ylim([0 110]) - catch - end + % try + % yyaxis right + % plot(obj.t0d(tmin:tmax),abs(obj.npart(tmin:tmax)./obj.npart(1)*100),'d--') + % ylabel('Nparts %') + % %ylim([0 110]) + % catch + % end ylimits=ylim; for i=1:length(obj.restarttimes) plot(obj.restarttimes(i)*[1 1],ylimits,'k--') end - + set(gca,'fontsize',14) + Position=f.Position; + Position(3)=1.5*Position(3); + f.Position=Position; obj.savegraph(f,sprintf('%s/%sEnergy',obj.folder,obj.name)); end function f=displaycharge(obj,f,linelegend) %% Plot the time evolution of the system charge of electrons % f: figure handle if you want to stack several such curves % on the same figure % linelegend: legend for this charge evolution tmin=1; tmax=length(obj.ekin); if nargin<2 f=figure('Name', sprintf('%s Charge',obj.name)); end if nargin<3 linelegend=''; end ax=f.CurrentAxes; if isempty(ax) ax=axes(f); end try plot(ax,obj.t0d(tmin:tmax),abs(obj.npart(tmin:tmax)*obj.qsim),'linewidth',1.5,'displayname',linelegend) hold on ylabel(ax,'Total charge [C]') xlabel(ax,'t [s]') grid on if(nargin>2) legend end set(ax,'fontsize',12) catch end if nargin < 2 obj.savegraph(f,sprintf('%s/%scharge',obj.folder,obj.name)); end end function displaySimParticles(obj) %% Plot the time evolution of the number of simulated markers in the main specie f=figure('Name', sprintf('%s Trapped particles',obj.name)); plot(obj.t0d,obj.npart,'linewidth',1.5) xlabel('t [s]') ylabel('N particles') set(gca,'fontsize',12) obj.savegraph(f,sprintf('%s/%sntrapped',obj.folder,obj.name),[10 12]); end function displayLarmorRad(obj,time2d) if nargin<2 time2d=length(obj.t2d); end % Plot the larmor radius for created particles with low energy % the larmor radius is calculated by considering that the % initial perpendicular velocity \approx the ExB velocity if time2d>0 Er=obj.Er(:,:,time2d); Ez=obj.Ez(:,:,time2d); else Er=obj.Erxt(:,:,1); Ez=obj.Ezxt(:,:,1); end rl=abs(obj.me/obj.qe*(-Er.*obj.Bz'+Ez.*obj.Br')./(obj.B.^3)'); figure rl(obj.geomweight(:,:,1)<0)=0; contourf(obj.zgrid,obj.rgrid,rl) hold on contour(obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0],'r-','linewidth',3) if time2d>0 n=obj.N(:,:,time2d); maxN=max(n (:)); n=n/maxN*mean(rl(:)); contour(obj.zgrid,obj.rgrid,n,linspace(0,1,6)*mean(rl(:)),'r:','linewidth',3) end c=colorbar; xlabel('z [m]') ylabel('r [m]') c.Label.String='r_L [m]'; end function displayHP(obj,tstart) % Plot the histogramm of the total energy and canonical angular momentum at time tstart and % end time of the simulation over the full simulation space for the main specie if(iscell(tstart)) tstart=cell2mat(tstart); end if(obj.species(1).R.nt>=2) tstart=obj.species(1).R.nt; f=figure('Name', sprintf('%s HP',obj.name)); legtext=sprintf("t=%2.1f - %2.1f [ns]",obj.tpart(tstart)*1e9,obj.tpart(end)*1e9); subplot(1,2,1) partsmax=min(obj.nbparts(end),obj.species(1).R.nparts); Hloc=H(obj,{1:obj.nbparts(1),1,false}); h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9)); hold on Hloc=H(obj,{1:partsmax,obj.species(1).R.nt,false}); %,'Binwidth',h1.BinWidth h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',legtext); ylabel('counts') xlabel('H [J]') legend subplot(1,2,2) Ploc=P(obj,{1:obj.nbparts(1),1,false}); h2=histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9)); hold on Ploc=P(obj,{1:partsmax,obj.species(1).R.nt,false}); histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',legtext); ylabel('counts') xlabel('P [kg\cdotm^2\cdots^{-1}]') %clear P %clear H legend %xlim([0.95*h2.BinLimits(1) 1.05*h2.BinLimits(2)]) obj.savegraph(f,sprintf('%s/%sParts_HP',obj.folder,obj.name)); end end function displayaveragetemp(obj) % Computes and show the particles average temperature as a function of time f=figure('Name',sprintf('%s potinn=%f part temperature',obj.name,obj.potinn)); vr2=obj.species(1).VR(:,:,false); vr2=mean(vr2.^2,1)-mean(vr2,1).^2; vz2=obj.species(1).VZ(:,:,false); vz2=mean(vz2.^2,1)-mean(vz2,1).^2; vthet2=obj.species(1).VTHET(:,:,false); vthet2=mean(vthet2.^2,1)-mean(vthet2,1).^2; plot(obj.tpart,0.5*obj.me*vr2/obj.qe,'displayname','T_r') hold on plot(obj.tpart,0.5*obj.me*vz2/obj.qe,'displayname','T_z') plot(obj.tpart,0.5*obj.me*vthet2/obj.qe,'displayname','T_{thet}') xlabel('time [s]') ylabel('T [eV]') title(sprintf('\\phi_a=%.1f kV \\phi_b=%.1f kV R=%.1f',obj.potinn/1e3,obj.potout/1e3,obj.Rcurv)) legend grid obj.savegraph(f,sprintf('%s/%s_partstemp',obj.folder,obj.name)); end function displayCurrentsevol(obj,timesteps) % Computes and display the time evolution of the outgoing currents on each domain boundary % at timesteps timesteps if nargin<2 timesteps=1:length(obj.t2d); end currents=obj.OutCurrents(timesteps); f=figure('Name',sprintf('%s Currents',obj.name)); if(obj.B(1,1)>obj.B(end,1)) lname='HFS'; rname='LFS'; else lname='LFS'; rname='HFS'; end plot(obj.t2d(timesteps),currents(1,:),'Displayname',lname,'linewidth',1.8); hold on plot(obj.t2d(timesteps),currents(2,:),'Displayname',rname,'linewidth',1.8); plot(obj.t2d(timesteps),currents(3,:),'Displayname','outer cylinder','linewidth',1.8); plot(obj.t2d(timesteps),currents(4,:),'Displayname','inner cylinder','linewidth',1.8); if size(currents,1)>=5 plot(obj.t2d(timesteps),currents(5,:),'Displayname','ellipse','linewidth',1.8); end legend('location','Northeast') xlabel('time [s]') ylabel('I [A]') grid on set(gca,'fontsize',12) title(sprintf('\\phi_b-\\phi_a=%.2g kV, R=%.1f',(obj.potout-obj.potinn)/1e3,obj.Rcurv)) obj.savegraph(f,sprintf('%s/%s_outCurrents',obj.folder,obj.name),[16 12]); end function displayChargeLossevol(obj,timesteps,toptitle,scalet,dens) % Computes and display the time evolution of the outgoing currents on each domain boundary % at time obj.t2d(timesteps) %scalet=true scales the time by the ellastic collision %frequency %dens = true plot the time evolution of the maximum electron %density in the simulation domain otherwise plot the total %number of electrons in the domain if nargin<2 timesteps=1:length(obj.t2d); end if nargin<4 scalet=true; end if nargin <5 dens=true; end if scalet if obj.neutcol.present vexb0=(obj.Ez(:,:,1).*obj.Br'-obj.Er(:,:,1).*obj.Bz')./(obj.B'.^2); vexb0(obj.geomweight(:,:,1)<=0)=0; E=0.5*obj.msim/obj.weight*mean(abs(vexb0(:)))^2/obj.qe; taucol=1/(obj.neutcol.neutdens*mean(abs(vexb0(:)))*(obj.sigio(E)+obj.sigmela(E)+obj.sigmio(E))); try Sio_S=1e17*(obj.neutcol.neutdens*mean(abs(vexb0(:)))*obj.sigio(E))/(obj.maxwellsrce.frequency*obj.weight/(pi*(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2)*diff(obj.maxwellsrce.zlim))) catch end tlabel='t/\tau_d [-]'; else taucol=2*pi/obj.omece; tlabel='t/\tau_ce [-]'; end else taucol=1e-9; tlabel='t [ns]'; end if dens N=obj.N(:,:,timesteps); geomw=obj.geomweight(:,:,1); geomw(geomw<0)=0; geomw(geomw>0)=1; N=N.*geomw; nmax=squeeze(max(max(N,[],1),[],2)); tn=(obj.t2d(timesteps)); nlabel='n_{e,max} [m^{-3}]'; ndlabel='n_{e,max}'; else t0dst=find(obj.t0d>=obj.t2d(timesteps(1)),1,'first'); t0dlst=find(obj.t0d<=obj.t2d(timesteps(end)),1,'last'); tn=obj.t0d(t0dst:t0dlst); nmax=obj.npart(t0dst:t0dlst)*obj.weight; nlabel='Nb e^-'; ndlabel='Nb e^-'; end [currents,pos]=obj.OutCurrents(timesteps); P=obj.neutcol.neutdens*obj.kb*300/100;% pressure at room temperature in mbar currents=currents/P; f=figure('Name',sprintf('%s Charges',obj.name)); % Plot the evolution of nb of particles yyaxis right p=plot(tn/taucol,nmax,'b-.','linewidth',1.8,'Displayname',ndlabel); ylabel(nlabel) ax=gca; ax.YAxis(2).Color=p.Color; ylim([0 inf]) if(obj.B(1,1)>obj.B(end,1)) lname='HFS'; rname='LFS'; else lname='LFS'; rname='HFS'; end yyaxis left mincurr=max(currents(:))*5e-3; if (max(currents(1,:)>mincurr)) plot(obj.t2d(timesteps)/taucol,currents(1,:),'r:','Displayname',lname,'linewidth',1.8); end hold on if (max(currents(2,:)>mincurr)) plot(obj.t2d(timesteps)/taucol,currents(2,:),'r--','Displayname',rname,'linewidth',1.8); end if (max(currents(3,:)>mincurr)) plot(obj.t2d(timesteps)/taucol,currents(3,:),'r-','Displayname','outer cylinder','linewidth',1.8); end if (max(currents(4,:)>mincurr)) plot(obj.t2d(timesteps)/taucol,currents(4,:),'Displayname','inner cylinder','linewidth',1.8); end if (size(currents,1)>=5 && max(currents(5,:)>mincurr)) plot(obj.t2d(timesteps)/taucol,currents(5,:),'r-','Displayname','ellipse','linewidth',1.8); end xlabel(tlabel) ylabel('I/p_n [A/mbar]') grid on set(gca,'fontsize',12) ax.YAxis(1).Color='red'; legend('Orientation','horizontal','location','south','numcolumns',3) if nargin <3 title(sprintf('\\phi_b-\\phi_a=%.2g kV, B=%f T',(obj.potout-obj.potinn)/1e3,max(obj.B(:)))) elseif ~isempty(toptitle) title(toptitle) end obj.savegraph(f,sprintf('%s/%s_ChargeEvol%i%i',obj.folder,obj.name,scalet,dens),[16 14]); end function display1Dpotentialwell(obj,timestep,rpos) % Display the potential well along the magentic field line % passing by rgrid(rpos) at the center of the simulation space if iscell(timestep) timestep=cell2mat(timestep); end f=figure('Name',sprintf('%s 1D Potential well',obj.name)); model=obj.potentialwellmodel(timestep); z=model.z; r=model.r; Pot=model.pot; rathet=model.rathet; if (mod(rpos, 1) ~= 0) [~,rpos]=min(abs(M.rgrid-rpos)); end crpos=obj.rgrid(rpos); id=find(timestep==0); timestep(id)=[]; n=obj.N(:,:,timestep); if(~isempty(timestep==0)) N0=zeros(obj.N.nr+1,obj.N.nz+1); n=cat(3,n(:,:,1:id-1),N0,n(:,:,id:end)); end n=mean(n,3); linepot=zeros(length(obj.zgrid),length(timestep)); rathetpos=obj.rAthet(rpos,ceil(length(obj.zgrid)/2)); F=scatteredInterpolant(z',rathet',Pot(:,1)); for i=1:length(timestep) F=scatteredInterpolant(z',rathet',Pot(:,i)); linepot(:,i)=F(obj.zgrid,rathetpos*ones(size(obj.zgrid))); %linepot(:,i)=griddata(z,rathet,pot(:,i),obj.zgrid,rathetpos); end linepot=mean(linepot,2); [Zinit,~]=meshgrid(obj.zgrid,obj.rAthet(:,1)); n=griddata(Zinit,obj.rAthet,n,obj.zgrid,rathetpos); plot(obj.zgrid,linepot) ylabel('Potentiel [eV]') xlabel('z [m]') xlim([obj.zgrid(1) obj.zgrid(end)]) hold(gca, 'on') yyaxis right plot(obj.zgrid,n) ylabel('n [m^{-3}]') if length(timestep)==1 title(sprintf('Potential well t=%1.2f [ns] r=%1.2f [mm]',obj.t2d(timestep)*1e9,1e3*crpos)) else title(sprintf('Potential well t=[%1.2f-%1.2f] [ns] r=%1.2f [mm]',obj.t2d(timestep(1))*1e9,obj.t2d(timestep(end))*1e9,1e3*crpos)) end obj.savegraph(f,sprintf('%s/%s_well1Dr_%d',obj.folder,obj.name,rpos)); end function displayVdistribRThetZ(obj,timestep, rpos, zpos) %displayVdistribRThetZ plot the velocity distribution function % in m/s %extracted from the markers at position window from rpos(1) %rpos(end) and zpos(1) to zpos(end) % and at time obj.tpart(timestep) %rpos and zpos are given as grid indices if(obj.species(1).R.nt>=2) if nargin<2 timesteppart=length(obj.tpart); + timestep=timesteppart; else timesteppart=timestep; end if nargin<3 || isempty(rpos) rpos=1:length(obj.rgrid); rspan=[obj.rgrid(1) obj.rgrid(end)]; else r=[obj.rgrid(1);(obj.rgrid(1:end-1)+obj.rgrid(2:end))*0.5;obj.rgrid(end)]; rspan=[r(rpos) r(rpos+1)]; end if nargin<4 || isempty(zpos) zpos=1:length(obj.zgrid); zspan=[obj.zgrid(1) obj.zgrid(end)]; else z=[obj.zgrid(1);(obj.zgrid(1:end-1)+obj.zgrid(2:end))*0.5;obj.zgrid(end)]; zspan=[z(zpos) z(zpos+1)]; end nbp=min(obj.nbparts(1),obj.species(1).R.nparts); R=obj.species(1).R(1:nbp,1,false); Z=obj.species(1).Z(1:nbp,1,false); Vr=obj.species(1).VR(1:nbp,1,false); Vz=obj.species(1).VZ(1:nbp,1,false); Vthet=obj.species(1).VTHET(1:nbp,1,false); ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2); Vr=Vr(ids); Vz=Vz(ids); Vthet=Vthet(ids); vTr=std(Vr,1); vTz=std(Vz,1); vTthet=std(Vthet,1); nbp=min(obj.nbparts(timesteppart),obj.species(1).R.nparts); Rend=obj.species(1).R(1:nbp,timesteppart,false); Zend=obj.species(1).Z(1:nbp,timesteppart,false); Vrend=obj.species(1).VR(1:nbp,timesteppart,false); Vzend=obj.species(1).VZ(1:nbp,timesteppart,false); Vthetend=obj.species(1).VTHET(1:nbp,timesteppart,false); ids=Rend>=rspan(1) & Rend<=rspan(2) & Zend>=zspan(1) & Zend<=zspan(2); nbtot=sum(ids) Vrend=Vrend(ids); Vzend=Vzend(ids); Vthetend=Vthetend(ids); vTrend=std(Vrend,1); vTzend=std(Vzend,1); vTthetend=std(Vthetend,1); binwidth=abs(max(Vrend)-min(Vrend))/sqrt(length(Vrend)); f=figure('Name',sprintf("%s vrz distrib",obj.file)); [~,time2did]=min(abs(obj.t2d-obj.tpart(timestep))); subplot(1,4,1); obj.dispV(Vr,Vrend,'V_r [m/s]',[1,timesteppart]) [~,time2did]=min(abs(obj.t2d-obj.tpart(timestep))); if length(rpos)==1 vexb=-obj.Er(rpos,zpos,time2did)/obj.Bz(zpos,rpos)'; vexb=mean(vexb(:)); if ~isempty(obj.neutcol.ela_cross_sec) % plot the radial drift velocity as nu_dE_r/(B\Omega_c) vdr=obj.neutcol.neutdens*obj.sigmela(vexb^2*obj.me*0.5/obj.qe)*vexb*-obj.Er(rpos,zpos,time2did)... ./(obj.B(zpos,rpos)'.*obj.B(zpos,rpos)'*obj.qe/obj.me); vdr=mean(vdr(:)); ylimits=ylim; plot(vdr*[1 1],ylimits,'k--','displayname',sprintf('V_{d,pred}=%1.2g [m/s]',vdr)) end end subplot(1,4,2); obj.dispV(Vthet,Vthetend,'V\theta [m/s]',[1,timesteppart]) hold on drawnow ylimits=ylim; if length(rpos)==1 if ~isempty(obj.Erxt) vexbext=-obj.Erxt(rpos,zpos)/obj.Bz(zpos,rpos)'; plot(vexbext*[1 1],ylimits,'k--','displayname',sprintf('V_{ExB,ext}=%1.2g [m/s]',vexbext)) end plot(vexb*[1 1],ylimits,'k-.','displayname',sprintf('V_{ExB,tot}=%1.2g [m/s]',vexb)) end subplot(1,4,3); obj.dispV(Vz,Vzend,'Vz [m/s]',[1,timesteppart]) subplot(1,4,4); obj.dispV(sqrt(Vr.^2+(Vthet).^2+Vz.^2),sqrt(Vrend.^2+(Vthetend).^2+Vzend.^2),'Vtot [m/s]',[1,timesteppart],'maxwell') sgtitle(sprintf('t=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm]',obj.tpart(timestep)*1e9, rspan*1e3, zspan*1e3)) obj.savegraph(f,sprintf('%s/%sParts_V_RZ',obj.folder,obj.name),[25 14]); end end function displayEkin(obj,timestep, rpos, zpos) %displayEkin plot the kinetic energy distribution function in %eV %extracted from the markers at position window from rpos(1) %rpos(end) and zpos(1) to zpos(end) % and at time obj.tpart(timestep) %rpos and zpos are given as grid indices if(obj.species(1).R.nt>=2) if nargin<2 timesteppart=[1 length(obj.tpart)]; else if length(timestep)<2 timesteppart=[1 timestep]; else timesteppart=[timestep(1) timestep(end)]; end end if nargin<3 || isempty(rpos) rspan=[obj.rgrid(1) obj.rgrid(end)]; else r=[obj.rgrid(1);(obj.rgrid(1:end-1)+obj.rgrid(2:end))*0.5;obj.rgrid(end)]; rspan=[r(rpos) r(rpos+1)]; end if nargin<4 || isempty(zpos) zspan=[obj.zgrid(1) obj.zgrid(end)]; else z=[obj.zgrid(1);(obj.zgrid(1:end-1)+obj.zgrid(2:end))*0.5;obj.zgrid(end)]; zspan=[z(zpos) z(zpos+1)]; end nbp=min(obj.nbparts(timesteppart(1)),obj.species(1).R.nparts); R=obj.species(1).R(1:nbp,timesteppart(1),false); Z=obj.species(1).Z(1:nbp,timesteppart(1),false); Ekin=obj.Ekin(1:nbp,timesteppart(1),false); ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2); Ekin=Ekin(ids)/obj.qe; nbp=min(obj.nbparts(timesteppart(2)),obj.species(1).R.nparts); Rend=obj.species(1).R(1:nbp,timesteppart(2),false); Zend=obj.species(1).Z(1:nbp,timesteppart(2),false); Ekinend=obj.Ekin(1:nbp,timesteppart(2),false); ids=Rend>=rspan(1) & Rend<=rspan(2) & Zend>=zspan(1) & Zend<=zspan(2); Ekinend=Ekinend(ids)/obj.qe; f=figure('Name',sprintf("%s E_k distrib",obj.file)); obj.dispV(Ekin,Ekinend,'E_k',timesteppart,'maxwell') sgtitle(sprintf('dt=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm]',obj.dt*1e9, rspan*1e3, zspan*1e3)) obj.savegraph(f,sprintf('%s/%sParts_E_kin',obj.folder,obj.name)); end end function displayVdistribParPer(obj,timestep, rpos, zpos, gcs) %displayVdistribParPer plot the velocity distribution function % in m/s for the parallel and perpendicular velocity %extracted from the markers at position window from rpos(1) %rpos(end) and zpos(1) to zpos(end) % and at time obj.tpart(timestep) %rpos and zpos are given as grid indices % gcs define if you give the perpendicular velocity in the % guiding center frame or in the laboratory frame if(obj.species(1).R.nt>=2) if nargin<2 timesteppart=length(obj.tpart); else timesteppart=timestep; end if nargin<3 || isempty(rpos) rspan=[obj.rgrid(1) obj.rgrid(end)]; else r=[obj.rgrid(1);(obj.rgrid(1:end-1)+obj.rgrid(2:end))*0.5;obj.rgrid(end)]; rspan=[r(rpos) r(rpos+1)]; end if nargin<4 || isempty(zpos) zspan=[obj.zgrid(1) obj.zgrid(end)]; else z=[obj.zgrid(1);(obj.zgrid(1:end-1)+obj.zgrid(2:end))*0.5;obj.zgrid(end)]; zspan=[z(zpos) z(zpos+1)]; end if nargin<5 gcs=false; % define if we look in the guiding center system end nbp=min(obj.nbparts(1),obj.species(1).R.nparts); R=obj.species(1).R(1:nbp,1,false); Z=obj.species(1).Z(1:nbp,1,false); ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2); Vperp=obj.Vperp(1:nbp,1,false,gcs); Vpar=obj.Vpar(1:nbp,1,false); Vperp=Vperp(ids); Vpar=Vpar(ids); nbp=min(obj.nbparts(timesteppart),obj.species(1).R.nparts); R=obj.species(1).R(1:nbp,timesteppart,false); Z=obj.species(1).Z(1:nbp,timesteppart,false); ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2); Vperpend=obj.Vperp(1:nbp,timesteppart,false,gcs); Vparend=obj.Vpar(1:nbp,timesteppart,false); Vperpend=Vperpend(ids); Vparend=Vparend(ids); %binwidth=abs(max(Vparend)-min(Vparend))/500; f=figure('Name',sprintf("%s v parper distrib",obj.file)); subplot(1,2,1) if gcs lgd='v_\perp gcs [m/s]'; else lgd='v_\perp [m/s]'; end obj.dispV(Vperp,Vperpend,lgd,[1,timesteppart], 'maxwell') subplot(1,2,2) obj.dispV(abs(Vpar),abs(Vparend),'v_{par} [m/s]',[1,timesteppart],'None') sgtitle(sprintf('t=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm]',obj.tpart(timestep)*1e9, rspan*1e3, zspan*1e3)) obj.savegraph(f,sprintf('%s/%sParts_V_parper',obj.folder,obj.name)); if gcs obj.savegraph(f,sprintf('%s/%sParts_V_parpergcs',obj.folder,obj.name)); else obj.savegraph(f,sprintf('%s/%sParts_V_parper',obj.folder,obj.name)); end end end function display2DVdistrib(obj,timestep, rpos, zpos, gcs) %display2DVdistrib plot the velocity distribution function % in m/s for the parallel and perpendicular velocity % and for the radial azimuthal velocity % as a 2D contour plot the show the velocity phase space distribution %extracted from the markers at position window from rpos(1) %rpos(end) and zpos(1) to zpos(end) % and at time obj.tpart(timestep) %rpos and zpos are given as grid indices % gcs define if you give the perpendicular velocity in the % guiding center frame or in the laboratory frame if(obj.species(1).R.nt>=2) if nargin<2 timesteppart=length(obj.tpart); else timesteppart=timestep; end if nargin<3 || isempty(rpos) rspan=[obj.rgrid(1) obj.rgrid(end)]; else r=[obj.rgrid(1);(obj.rgrid(1:end-1)+obj.rgrid(2:end))*0.5;obj.rgrid(end)]; rspan=[r(rpos(1)) r(rpos(end)+1)]; end if nargin<4 || isempty(zpos) zspan=[obj.zgrid(1) obj.zgrid(end)]; else z=[obj.zgrid(1);(obj.zgrid(1:end-1)+obj.zgrid(2:end))*0.5;obj.zgrid(end)]; zspan=[z(zpos(1)) z(zpos(end)+1)]; end if nargin<5 gcs=false; % define if we look in the guiding center system end nbp=min(obj.nbparts(timesteppart),obj.species(1).R.nparts); R=obj.species(1).R(1:nbp,timesteppart,false); Z=obj.species(1).Z(1:nbp,timesteppart,false); ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2); Vperp=obj.Vperp(1:nbp,timesteppart,false,gcs); Vpar=obj.Vpar(1:nbp,timesteppart,false); Vr=obj.species(1).VR(1:nbp,timesteppart,false); Vthet=obj.species(1).VTHET(1:nbp,timesteppart,false); Vper=Vperp(ids); Vpar=Vpar(ids); Vr=Vr(ids); Vthet=Vthet(ids); nbp=sum(ids(:)); f=figure('Name',sprintf("%s v parper distrib",obj.file)); subplot(2,1,1) [N,Xedges,Yedges] = histcounts2(Vpar,Vper,20); Xedges=(Xedges(1:end-1)+Xedges(2:end))/2; Yedges=(Yedges(1:end-1)+Yedges(2:end))/2; contourf(Xedges,Yedges,N') xlabel('v_{par} [m/s]') ylabel('v_{\perp} [m/s]') c=colorbar; c.Label.String='Counts'; subplot(2,1,2) [N,Xedges,Yedges] = histcounts2(Vthet,Vr,20); Xedges=(Xedges(1:end-1)+Xedges(2:end))/2; Yedges=(Yedges(1:end-1)+Yedges(2:end))/2; contourf(Xedges,Yedges,N') %histogram2(Vthet,Vr,'displaystyle','tile','binmethod','auto') %scatter(Vthet,Vr) xlabel('v_\theta [m/s]') ylabel('v_r [m/s]') c=colorbar; c.Label.String='Counts'; sgtitle(sprintf('t=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm] N=%3i',mean(obj.tpart(timestep))*1e9, rspan*1e3, zspan*1e3,nbp)) mkdir(sprintf('%s/vdist',obj.folder)) if gcs obj.savegraph(f,sprintf('%s/vdist/%sParts_V_2dparpergcs_r%iz%it%i',obj.folder,obj.name,floor(mean(rpos)),floor(mean(zpos)),floor(mean(timestep)))); else obj.savegraph(f,sprintf('%s/vdist/%sParts_V_2dparper_r%iz%it%i',obj.folder,obj.name,floor(mean(rpos)),floor(mean(zpos)),floor(mean(timestep)))); end end end function [p, maxnb, c]=displayPhaseSpace(obj,type,partsstep, Rindex, Zindex,legendtext, figtitle, f, maxnb, c, gcs) if nargin<8 f=figure; f=gca; end if nargin<7 figtitle=sprintf('r=%1.2f [mm] z=%1.2f [mm] \\Delta\\phi=%1.1f[kV] R=%1.1f',obj.rgrid(Rindex)*1e3,obj.zgrid(Zindex)*1e3,(obj.potout-obj.potinn)*obj.phinorm/1e3,obj.Rcurv); end if nargin <6 legendtext=sprintf('t=%1.3g [s]',obj.tpart(partsstep)); end fieldstep=find(obj.tpart(partsstep(end))==obj.t2d,1); if nargin>=10 ctemp=c; n=zeros(length(c{1}),length(c{2})); else nbins=15; n=zeros(nbins); end if nargin <11 gcs=true; end for i=1:length(partsstep) odstep=find(obj.tpart(partsstep(i))==obj.t0d); nbp=min(obj.species(1).R.nparts,obj.nbparts(partsstep(i))); Rp=obj.species(1).R(1:nbp,partsstep(i),false); Zp=obj.species(1).Z(1:nbp,partsstep(i),false); deltar=obj.dr(2)/2; deltarm=obj.rgrid(Rindex)-sqrt(obj.rgrid(Rindex)^2-deltar^2-2*obj.rgrid(Rindex)*deltar); deltaz=obj.dz/2; Indices=Rp>=obj.rgrid(Rindex)-deltarm & Rp=obj.zgrid(Zindex)-deltaz & Zp0)),max(Blines(obj.geomweight(:,:,1)>0)),20); Blines(obj.geomweight(:,:,1)<0)=NaN; [~,h1]=contour(ax1,obj.zgrid*1000,obj.rgrid*1000,Blines,real(levels),'-.','color','k','linewidth',1.5,'Displayname','Magnetic field lines'); % Draw the metallic boundaries and the geometry itself [c1,hContour]=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,-geomw,[0,0],'linewidth',1.5); drawnow; xlim(ax1,[obj.zgrid(1)*1000 obj.zgrid(end)*1000]) % Change the color of the metallic boundaries to grey hFills=hContour.FacePrims; [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' try hFills(end).ColorData = uint8([150;150;150;255]); for idx = 1 : numel(hFills)-1 hFills(idx).ColorData(4) = 0; % default=255 end catch end - + grid on; hold on; f.PaperOrientation='landscape'; f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; set(ax1,'fontsize',14) %axis equal obj.savegraph(f,sprintf('%sfluid_dens',obj.name)) end function displaymagfield(obj) %displaymagfield display the magnetic field lines and the %amplitude of the magnetic field using a contour % also show the domain boundaries B=obj.B'; f=figure('Name', sprintf('%s B field',obj.name)); B(obj.geomweight(:,:,1)<0)=NaN; ax1=gca; title(sprintf('Configuration')) Blv=linspace(min(B(:)),max(B(:)),50); if length(Blv)<2 Blv=min(B(:))*[1 1]; end h=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,B,Blv,'Displayname','B [T]', 'linestyle','none'); hold on; %% Magnetic field lines Blines=obj.rAthet; levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),50); Blines(obj.geomweight(:,:,1)<0)=NaN; [~,h1]=contour(ax1,obj.zgrid*1000,obj.rgrid*1000,Blines,real(levels),'r-.','linewidth',1.5,'Displayname','Magnetic field lines'); colormap(ax1,'parula') % Grey outline geomw=obj.geomweight(:,:,1); geomw(geomw>0)=NaN; geomw(geomw<0)=min(B(:)); [c1,hContour]=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,geomw, [0 0]); drawnow; xlim(ax1,[obj.zgrid(1)*1000 obj.zgrid(end)*1000]) if(obj.conformgeom) ylim([ax1 ],[obj.rgrid(1)*1000 obj.rgrid(rgridend)*1000]) else ylim([ax1],[obj.rgrid(1)*1000 obj.rgrid(end)*1000]) end legend([h1],{'Magnetic field lines'},'location','northwest') xlabel(ax1,'z [mm]') ylabel(ax1,'r [mm]') %title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),double(maxdens))) c = colorbar(ax1); c.Label.String= 'B [T]'; view(ax1,2) %set(h,'edgecolor','none'); grid on; hFills=hContour.FacePrims; [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' %caxis([min(B(:)) max(B(:))]) try hFills(1).ColorData = uint8([150;150;150;255]); for idx = 2 : numel(hFills) hFills(idx).ColorData(4) = 0; % default=255 end catch end [~, name, ~] = fileparts(obj.file); if( obj.walltype >=2 && obj.walltype<=4) rectangle('Position',[obj.zgrid(1) obj.r_b obj.zgrid(end)-obj.zgrid(1) 0.001]*1e3,'FaceColor',[150 150 150]/255,'Edgecolor','none') ylimits=ylim; ylim([ylimits(1),ylimits(2)+1]) end if(isempty(obj.spl_bound)) rectangle('Position',[obj.zgrid(1) obj.r_a-0.001 obj.zgrid(end)-obj.zgrid(1) 0.001]*1e3,'FaceColor',[150 150 150]/255,'Edgecolor','none') ylimits=ylim; ylim([ylimits(1)-1,ylimits(2)]) end f.PaperUnits='centimeters'; %axis equal papsize=[14 5 ]; pos=f.Position; pos(3)=floor(1.3*pos(3)); f.Position=pos; obj.savegraph(f,sprintf('%s/%s_Bfield',obj.folder,obj.name),papsize); end function displaySurfFlux(obj,timestep, subdiv) %displaySurfFlux plot the current densities %on the domain boundaries for time t2d(timestep) %directly on the boundaries themselves %make it easier to see where the currents are collected if nargin<3 subdiv=1; end mflux= obj.Metallicflux(timestep,subdiv); lflux= -squeeze(obj.Axialflux(timestep,1))'; rflux= squeeze(obj.Axialflux(timestep,length(obj.zgrid)))'; time=obj.t2d(timestep); if nargin<3 ids=1:length(mflux); end %% P=obj.neutcol.neutdens*obj.kb*300/100;% pressure at room temperature in mbar f=figure('name','fluxevol'); linew=5; %obj.displaysplbound(gca,1e3); contour(obj.zgrid*1e3,obj.rgrid*1e3,obj.geomweight(:,:,1),[0 0],'b-','linewidth',1.5); hold on for i=1:length(mflux.p) x=mflux.p{i}(1,:)*1000; y=mflux.p{i}(2,:)*1000; y(end)=NaN; c=mean(mflux.gamma{i},2)'*obj.qe/(100^2)/P; c(c<=0)=NaN; patch(x,y,c,'EdgeColor','interp','LineWidth',linew); hold on end x=obj.zgrid(1)*ones(size(obj.rgrid))*1000; y=obj.rgrid*1000; y(end)=NaN; c=mean(lflux,1)*obj.qe/(100^2)/P; c(c<=0)=NaN; patch(x,y,c,'EdgeColor','interp','LineWidth',linew); x=obj.zgrid(end)*ones(size(obj.rgrid))*1e3; y=obj.rgrid*1000; y(end)=NaN; c=mean(rflux,1)*obj.qe/(100^2)/P; c(c<=0)=NaN; patch(x,y,c,'EdgeColor','interp','LineWidth',linew); title(sprintf('t=%4.2f [ns]',mean(time)*1e9)) c=colorbar; c.Label.String= 'j\cdotn [A/(cm^2 mbar)]'; xlabel('z [mm]') ylabel('r [mm]') colormap(jet) set(gca,'colorscale','log') %% Magnetic field lines Blines=obj.rAthet; levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),20); Blines(obj.geomweight(:,:,1)<0)=NaN; [~,h1]=contour(obj.zgrid*1000,obj.rgrid*1000,Blines,real(levels),'m-.','linewidth',1.5,'Displayname','Magnetic field lines'); %axis equal obj.savegraph(f,sprintf('%s/%s_surfFlux_it2d_%i',obj.folder,obj.name,floor(mean(timestep))),[16 14]); end % interactive window to display the terms of the pressure tensor - dispespicPressure(obj,logdensity,showgrid,fixed,temperature) + guiPressure(obj,logdensity,showgrid,fixed,temperature) % interactive window to display the electron density, magnetic % field lines, electric potential and field at given time steps - dispespicFields(obj,logdensity,showgrid,fixed,parper) + guiFields(obj,logdensity,showgrid,fixed,parper) function displaycollfreq(obj) %displaycollfreq plot the collision frequencies in Hz/mbar for a range of %electron kinetic energies in eV for the different collision %processes considered: ionisation and elastic collisions E=logspace(1,4,1000); v=sqrt(2/obj.msim*obj.weight*E*obj.qe); P=obj.neutcol.neutdens*obj.kb*300/100;% pressure at room temperature in mbar tauio=P./(obj.neutcol.neutdens*obj.sigio(E).*v); tauiom=P./(obj.neutcol.neutdens*obj.sigmio(E).*v); tauelam=P./(obj.neutcol.neutdens*obj.sigmela(E).*v); f=figure('name','t scales coll'); loglog(E,1./tauio,'displayname','ionisation','linewidth',1.5) hold on loglog(E,1./tauiom,'displayname','ionisation momentum','linewidth',1.5) loglog(E,1./tauelam,'displayname','elastic','linewidth',1.5) loglog(E,1./tauio+1./tauiom+1./tauelam,'displayname','total drag','linewidth',1.5) loglog([E(1) E(end)],1./(2*pi/obj.omece)* [1 1],'--','displayname','cyclotronic','linewidth',1.5) xlabel('Electron kinetic energy [eV]') ylabel('\nu [Hz/mbar]') legend('location','southeast') grid on obj.savegraph(f,sprintf('%s/collfreqscales',obj.folder),[14 12]); end function displaycrosssec(obj) %displaycrosssec plot the collision crosssections in m^2 for a range of %electron kinetic energies in eV for the different collision %processes considered: ionisation and elastic collisions E=logspace(1,4,1000); sig_io=obj.sigio(E); sig_iom=obj.sigmio(E); sig_elam=obj.sigmela(E); f=figure('name','t scales coll'); loglog(E,sig_io,'displayname','ionisation','linewidth',1.5) hold on loglog(E,sig_iom,'displayname','ionisation momentum','linewidth',1.5) loglog(E,sig_elam,'displayname','elastic','linewidth',1.5) loglog(E,sig_io+sig_elam+sig_iom,'displayname','total drag','linewidth',1.5) xlabel('Energy [eV]') ylabel('\sigma [m^{2}]') legend('location','southeast') grid on obj.savegraph(f,sprintf('%s/coll_cross_sec_scales',obj.folder),[14 12]); end %------------------------------------------ % Helper functions needed for other functions function [zpos,rpos]=getpos(obj,tstep) % interactive window to return an specific axial and radial % position picked from the cloud density if nargin<2 tstep=length(obj.t2d); end n=obj.N(:,:,tstep); n(obj.geomweight(:,:,1)<0)=NaN; figure contourf(obj.zgrid,obj.rgrid,n); xlabel('z [m]') ylabel('r [m]') [x,y]=ginput(1); zpos=find(x>obj.zgrid,1,'last'); rpos=find(y>obj.rgrid,1,'last'); hold on plot(obj.zgrid(zpos),obj.rgrid(rpos),'rx') fprintf('zpos=%i rpos=%i z=%1.4f r=%1.4f\n',zpos,rpos,obj.zgrid(zpos),obj.rgrid(rpos)) end function changed=ischanged(obj) %ischanged Check if the file has been changed since the initial loading of the file %and if some data must be reloaded try filedata=dir(obj.fullpath); checkedtimestamp=filedata.date; if (max(checkedtimestamp > obj.timestamp) ) changed=true; return end changed=false; return catch changed=true; return end end function dispV(obj,V,Vend,label,t, dist, vd) %dispV generic functio to plot the velocity distribution and %comapare two timesteps V and Vend at time t(1) and t(2) if nargin<6 dist='gaussian'; end if nargin<7 vd=0; end vmean=mean(V(~isnan(V))); vtherm=std(V(~isnan(V)),1); vmeanend=mean(Vend(~isnan(Vend))); vthermend=std(Vend(~isnan(Vend)),1); if(length(V)>1) [Counts,edges]=histcounts(V,'binmethod','sqrt'); binwidth=mean(diff(edges)); plot([edges(1) 0.5*(edges(2:end)+edges(1:end-1)) edges(end)],[0 Counts 0],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(t(1))*1e9)); hold on end hold on [Counts,edges]=histcounts(Vend,'binmethod','sqrt'); plot([edges(1) 0.5*(edges(2:end)+edges(1:end-1)) edges(end)],[0 Counts 0],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(t(2))*1e9)); if strcmp(dist,'maxwell') vfit=linspace(0,edges(end),300); a=vmeanend/sqrt(2); dist=sqrt(2/pi)*vfit.^2.*exp(-((vfit).^2-vd^2)/2/a^2)/a^3; dist=dist/max(dist); plot(vfit,max(Counts)*dist,'displayname',sprintf('Maxw mu=%2.2g sigma=%2.2g',vmeanend,vthermend)) elseif strcmp(dist,'gaussian') vfit=linspace(edges(1),edges(end),300); dist=exp(-(vfit-vmeanend).^2/2/vthermend^2); plot(vfit,max(Counts)*dist,'displayname',sprintf('gauss mu=%2.2g sigma=%2.2g',vmeanend,vthermend)) end ylabel('counts') xlabel(label) grid on legend('location','southoutside','orientation','vertical') end function cross_sec=fit_cross_sec(obj,energy,crosssec_table) %Interpolate the cross-section at the given energy using the %crosssec_table and an exponential fitting cross_sec=0; if (energy<=0 || isnan(energy) || isinf(energy)) return end id=find(energy>crosssec_table(:,1),1,'last'); if(isempty(id)) id=1; end id=min(size(crosssec_table,1)-1,id); id=max(1,id); cross_sec=crosssec_table(id,2)*(energy/crosssec_table(id,1))^crosssec_table(id,3); end function fighandle=savegraph(obj, fighandle, name, papsize) %% Saves the given figure as a pdf a .fig and an eps using export_fig fighandle.PaperUnits='centimeters'; if (nargin < 4) papsize=[14 16]; end set(fighandle, 'Color', 'w'); fighandle.PaperSize=papsize; %export_fig(fighandle,name,'-png','-r300') exportgraphics(fighandle,sprintf('%s.png',name),'Resolution',300) print(fighandle,name,'-dpdf','-fillpage') savefig(fighandle,name) set(fighandle, 'Color', 'w'); exportgraphics(fighandle,sprintf('%s.eps',name)) %export_fig(fighandle,name,'-eps','-painters') end function sig=dsigmaio(obj,Ekin, Ebar, Ei, E0, chi, gamma) % calculates the integrand used for the ionisation collision % cross section for momentum exchange for the incoming electron % it is only used by obj.sigmiopre gamma=reshape(gamma,1,[],1); chi=reshape(chi,1,1,[]); siggamma=sin(gamma).*(E0^2+8*(1-chi)*(Ekin-Ei)*E0)./(E0+4*(1-chi)*(Ekin-Ei)-4*(1-chi)*(Ekin-Ei).*cos(gamma)).^2/2; sigchi=(Ekin-Ei)./(Ebar*atan((Ekin-Ei)/(2*Ebar)).*(1+(chi*(Ekin-Ei)/Ebar).^2)); dp=1- trapz(gamma,sqrt((1-chi).*(1-Ei/Ekin)).*cos(gamma).*siggamma,2);%- trapz(gamma,sqrt((1-chi).*(1-Ei/Ekin)).*cos(gamma).*siggamma,2); sig=sigchi.*dp; end function sigm=sigmiopre(obj,E, init) % returns the precalculated values used for the interpolation % of the ionisation collision cross-section for momentum % exchange for the incoming electron if nargin <3 init=false; end if(~init &&( ~obj.neutcol.present || isempty(obj.neutcol.io_cross_sec))) sigm=zeros(size(E)); return end Ebar=obj.neutcol.scatter_fac; Ei=obj.neutcol.Eion; E0=obj.neutcol.E0; nE=numel(E); nchi=300; ngamma=300; gamma=linspace(0,pi,ngamma); chi=linspace(0,0.5,nchi); %sigm2=zeros(nE,nchi); sigm=zeros(size(E)); for i=1:nE if(E(i)>=Ei) sigm2=zeros(nchi,1); for j=1:nchi %sigm2(j)=trapz(alpha,trapz(gamma,obj.dsigmaio(E(i),Ebar,Ei,E0,chi(j),alpha,gamma),2),1); sigm2(j)=obj.dsigmaio(E(i),Ebar,Ei,E0,chi(j),gamma); end sigm(i)=trapz(chi,sigm2)*obj.sigio(E(i),init); %sigm(i)=trapz(chi,trapz(alpha,trapz(gamma,dsigmaio(obj,E(i),Ebar,Ei,E0,chi,alpha,gamma),2),1),3)*obj.sigio(E(i),init); end end end end end diff --git a/matlab/@fennecshdf5/guiRProf.m b/matlab/@fennecshdf5/guiRProf.m index c62ecc1..3f99701 100644 --- a/matlab/@fennecshdf5/guiRProf.m +++ b/matlab/@fennecshdf5/guiRProf.m @@ -1,331 +1,351 @@ function guiRProf(M,logdensity,showgrid,fixed,parper) %dispespicFields Allows to display the time evolution of the density, electric potential and electric fields % M is of class fennecshdf5 and contains the simulation results fieldstep=1; zpos=1; if nargin <2 logdensity=false; showgrid=false; fixed=false; end if nargin <3 showgrid=false; fixed=false; end if nargin <4 fixed=false; end if nargin <5 parper=false; end fixed=fi(fixed); f=uifigure('Name',sprintf('Grid data %s',M.name)); mf=uipanel(f,'Position',[5 50 f.Position(3)-30 f.Position(4)-55]); mf.AutoResizeChildren='off'; m=uipanel(f,'Position',[5 5 f.Position(3)-10 40]); rpanel=uipanel(f,'Position',[f.Position(3)-28 50 25 f.Position(4)-55]); sgtitle(mf,sprintf('step=%d t=%0.5e s',fieldstep*M.it1,M.t2d(fieldstep))) sld = uislider(m,'Position',[10 30 0.6*m.Position(3) 3]); sld.Value=fieldstep; sld.Limits=[1 size(M.t2d,1)]; sld.Tag='timeslider'; sldr = uislider(rpanel,'Orientation','vertical','Position',[5 35 40 rpanel.Position(4)-40]); sldr.Value=zpos; sldr.Limits=[1 length(M.zgrid)]; sldr.Tag='axialslider'; edr = uieditfield(rpanel,'numeric','Limits',[1 length(M.zgrid)],'Value',1); edr.Position=[sldr.Position(1) sldr.Position(2)-30 40 20]; edr.RoundFractionalValues='on'; edr.Tag='axialfield'; edt = uieditfield(m,'numeric','Limits',[1 size(M.t2d,1)],'Value',1); edt.Position=[sld.Position(1)+sld.Position(3)+25 5 40 20]; edt.RoundFractionalValues='on'; edt.Tag='timefield'; MaxN=0; Printbt=uibutton(m,'Position',[edt.Position(1)+edt.Position(3)+10 5 40 20],'Text', 'Save'); Play=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play'); Pause=uibutton(m,'Position',[Play.Position(1)+Play.Position(3)+10 5 40 20],'Text', 'Pause'); %Playbt=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play/Pause'); stop=false; sld.ValueChangingFcn={@updatefigdata,edt,mf}; edt.ValueChangedFcn={@updatefigdata,sld,mf}; sldr.ValueChangingFcn={@updatefigdata,edr,mf}; edr.ValueChangedFcn={@updatefigdata,sldr,mf}; Printbt.ButtonPushedFcn={@plotGridButtonPushed}; Play.ButtonPushedFcn={@plotPlayButtonPushed}; Pause.ButtonPushedFcn={@PauseButtonPushed}; set(f,'KeyPressFcn',{ @onKeyDown,sld,edt,mf}) Plotfennecsgriddata(mf,M,fieldstep,zpos); function plotPlayButtonPushed(btn,ax) stop=false; i=sld.Value; while ~stop edt.Value=i; sld.Value=i; updatesubplotsdata(i,mf); pause(0.01) i=sld.Value; i=i+10; if(i>sld.Limits(2)) stop=true; end end end function PauseButtonPushed(btn,ax) stop = true; end function onKeyDown(src,event,slider,editfield, fig) direction=0; if strcmp(event.Key,'leftarrow') direction=-1; elseif strcmp(event.Key,'rightarrow') direction=+1; elseif strcmp(event.Key,'uparrow') direction=+10; elseif strcmp(event.Key,'downarrow') direction=-10; end if(direction~=0) currval=slider.Value; slider.Value=max(slider.Limits(1),min(currval+direction,slider.Limits(2))); updatefigdata(slider, event, editfield ,fig) end end function Plotfennecsgriddata(fig,M,fieldstep,zpos) %Plotfennecsgriddata Plot the 2d data of fennecs at time step fieldstep sgtitle(fig,sprintf('step=%d t=%0.5e s z=%0.3e mm',(fieldstep-1)*M.it1,M.t2d(fieldstep),M.zgrid(zpos)*1e3)) geomw=M.geomweight(:,zpos,1)>=0; ax1=subplot(2,2,1,'Parent',fig); p=plot(ax1,M.rgrid*1e3,M.N(:,zpos,fieldstep),'linewidth',1.5); xlim(ax1,[M.rgrid(1) M.rgrid(end)]*1e3) xlabel(ax1,'r [mm]') title(ax1,'Density') ylabel(ax1,'n[m^{-3}]'); %c.Limits=[0 max(M.N(:))]; hold(ax1, 'on') [~,id1]=min(abs(M.geomweight(1:10,zpos,1))); [~,id2]=min(abs(M.geomweight(11:end,zpos,1))); id2=id2+10; - ylimits=ylim; + ylimits=ylim(ax1); plot(ax1,M.rgrid(id1)*[1 1]*1e3,ylimits,'k--','linewidth',1.5,'Displayname','Boundaries'); plot(ax1,M.rgrid(id2)*[1 1]*1e3,ylimits,'k--','linewidth',1.5,'Displayname','Boundaries'); yyaxis(ax1,'right') hold(ax1, 'on') Er=M.Er(:,zpos,fieldstep).*geomw; Ez=M.Ez(:,zpos,fieldstep).*geomw; + Er0=M.Erxt(:,zpos,1).*geomw; + Ez0=M.Ezxt(:,zpos,1).*geomw; p1=plot(ax1,M.rgrid*1e3,Er,'linewidth',1.5); p2=plot(ax1,M.rgrid*1e3,Ez,'linewidth',1.5); + p3=plot(ax1,M.rgrid*1e3,Er0,'linewidth',1.5); + p4=plot(ax1,M.rgrid*1e3,Ez0,'linewidth',1.5); ylabel(ax1,'E [V/m]') if max(abs([Er(:); Ez(:)]))>0 ylim(ax1,[ -max(abs([Er(:); Ez(:)])) max(abs([Er(:); Ez(:)]))]) end - legend(ax1,[p p1 p2],{'n','Er','Ez'},'location','northwest') + legend(ax1,[p p1 p2 p3 p4],{'n','Er','Ez','Erxt','Ezxt'},'location','northwest') ax2=subplot(2,2,2,'Parent',fig); ur=M.fluidUR(:,zpos,fieldstep); plot(ax2,M.rgrid*1e3,ur,'linewidth',1.5); xlim(ax2,[M.rgrid(1) M.rgrid(end)]*1e3) xlabel(ax2,'r [mm]') title(ax2,'radial velocity') ylabel(ax2,'v_r [m/s]'); %c.Limits=[0 max(M.N(:))]; hold(ax2, 'on') if max(ur)>0 ylim(ax2,[ -max(ur) max(ur)]) end %ylim(ax2,[ -max(ur) max(ur)]) - ylimits=ylim; + ylimits=ylim(ax2); plot(ax2,M.rgrid(id1)*[1 1]*1e3,ylimits,'k--','linewidth',1.5,'Displayname','Boundaries'); plot(ax2,M.rgrid(id2)*[1 1]*1e3,ylimits,'k--','linewidth',1.5,'Displayname','Boundaries'); ax3=subplot(2,2,3,'Parent',fig); uthet=M.fluidUTHET(:,zpos,fieldstep); plot(ax3,M.rgrid*1e3,uthet,'linewidth',1.5); xlim(ax3,[M.rgrid(1) M.rgrid(end)]*1e3) xlabel(ax3,'r [mm]') title(ax3,'Azimuthal velocity') ylabel(ax3,'v_\theta [m/s]'); %c.Limits=[0 max(M.N(:))]; hold(ax3, 'on') if max(uthet)>0 ylim(ax3,[ -max(uthet) max(uthet)]) end - ylimits=ylim; + ylimits=ylim(ax3); plot(ax3,M.rgrid(id1)*[1 1]*1e3,ylimits,'k--','linewidth',1.5,'Displayname','Boundaries'); plot(ax3,M.rgrid(id2)*[1 1]*1e3,ylimits,'k--','linewidth',1.5,'Displayname','Boundaries'); uExb=-M.Er(:,zpos,fieldstep)./M.Bz(zpos,:)'.*(uthet~=0); plot(ax3,M.rgrid*1e3,uExb,'linewidth',1.5); ax4=subplot(2,2,4,'Parent',fig); uz=M.fluidUZ(:,zpos,fieldstep); plot(ax4,M.rgrid*1e3,uz,'linewidth',1.5); xlim(ax4,[M.rgrid(1) M.rgrid(end)]*1e3) xlabel(ax4,'r [mm]') title(ax4,'Axial velocity') ylabel(ax4,'v_z [m/s]'); %c.Limits=[0 max(M.N(:))]; hold(ax4, 'on') if max(uz)>0 ylim(ax4,[ -max(uz) max(uz)]) end - ylimits=ylim; + ylimits=ylim(ax4); plot(ax4,M.rgrid(id1)*[1 1]*1e3,ylimits,'k--','linewidth',1.5,'Displayname','Boundaries'); plot(ax4,M.rgrid(id2)*[1 1]*1e3,ylimits,'k--','linewidth',1.5,'Displayname','Boundaries'); linkaxes([ax1,ax2,ax3,ax4],'x'); grid([ax1 ax2 ax3 ax4],'minor') end function plotGridButtonPushed(btn,ax) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here f=figure(); Plotfennecsgriddata(f,M,sld.Value,edr.Value); f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); print(f,sprintf('%sGrid%d%d',name,sld.Value,edr.Value),'-dpdf','-fillpage') end function updatefigdata(control, event, Othercontrol, fig) if contains(event.Source.Tag,'time') if strcmp(event.EventName,'ValueChanged') fieldstep=floor(control.Value); control.Value=fieldstep; elseif strcmp(event.EventName,'KeyPress') fieldstep=floor(control.Value); control.Value=fieldstep; else fieldstep=floor(event.Value); end Othercontrol.Value=fieldstep; elseif contains(event.Source.Tag,'axial') if strcmp(event.EventName,'ValueChanged') zpos=floor(control.Value); control.Value=zpos; elseif strcmp(event.EventName,'KeyPress') zpos=floor(control.Value); control.Value=zpos; else zpos=floor(event.Value); end Othercontrol.Value=zpos; end updatesubplotsdata(fieldstep, zpos, fig); end function updatesubplotsdata(fieldstep, zpos, fig) sgtitle(fig,sprintf('step=%d t=%0.5e s z=%0.3e mm',(fieldstep-1)*M.it1,M.t2d(fieldstep),M.zgrid(zpos)*1e3)) [~,rcenterid]=max(M.geomweight(:,zpos,1)); [~,id1]=min(abs(M.geomweight(1:rcenterid,zpos,1))); [~,id2]=min(abs(M.geomweight(rcenterid:end,zpos,1))); id2=id2+rcenterid; + try rlim1=M.rgrid(id1)*[1 1]*1e3; rlim2=M.rgrid(id2)*[1 1]*1e3; + catch + end %% update density ax1=fig.Children(end); geomw=M.geomweight(:,zpos,1)>=0; dens=M.N(:,zpos,fieldstep).*geomw; Er=M.Er(:,zpos,fieldstep).*geomw; Ez=M.Ez(:,zpos,fieldstep).*geomw; + Er0=M.Erxt(:,zpos,1).*geomw; + Ez0=M.Ezxt(:,zpos,1).*geomw; yyaxis(ax1,'left') ax1.Children(end).YData=dens; ylimits=ylim(ax1); - ax1.Children(end-1).XData=rlim1; - ax1.Children(end-1).YData=ylimits; - ax1.Children(end-2).XData=rlim2; - ax1.Children(end-2).YData=ylimits; - + try + ax1.Children(end-1).XData=rlim1; + ax1.Children(end-1).YData=ylimits; + ax1.Children(end-2).XData=rlim2; + ax1.Children(end-2).YData=ylimits; + catch + end yyaxis(ax1,'right') ax1.Children(end).YData=Er; ax1.Children(end-1).YData=Ez; + ax1.Children(end-2).YData=Er0; + ax1.Children(end-3).YData=Ez0; if max(abs([Er; Ez]))>0 ylim(ax1,1.05*max(abs([Er; Ez]))*[ -1 1]) end % view(ax1,2) %% update Radial velocity ax2=fig.Children(end-2); ur=M.fluidUR(:,zpos,fieldstep).*geomw; ax2.Children(end).YData=ur; if max(abs(ur))>0 ylim(ax2,max(abs(ur))*[ -1 1]) end + try ax2.Children(end-1).XData=rlim1; ax2.Children(end-2).XData=rlim2; ylimits=ylim(ax2); ax2.Children(end-1).YData=ylimits; ax2.Children(end-2).YData=ylimits; - + catch + end %% update Azimuthal velocity ax3=fig.Children(end-3); uthet=M.fluidUTHET(:,zpos,fieldstep).*geomw; uExb=-M.Er(:,zpos,fieldstep)./M.Bz(zpos,:)'.*(uthet~=0); ax3.Children(end-3).YData=uExb'; ax3.Children(end).YData=uthet; if max(abs(uthet))>0 ylim(ax3,1.05*max(abs(uthet))*[-1 1]) end - + try ax3.Children(end-1).XData=rlim1; ax3.Children(end-2).XData=rlim2; ylimits=ylim(ax3); ax3.Children(end-1).YData=ylimits; ax3.Children(end-2).YData=ylimits; + catch + end %% update Axial velocity ax4=fig.Children(end-4); uz=M.fluidUZ(:,zpos,fieldstep).*geomw; ax4.Children(end).YData=uz; if max(abs(uz))>0 ylim(ax4,1.05*max(abs(uz))*[ -1 1]) end + try ax4.Children(end-1).XData=rlim1; ax4.Children(end-2).XData=rlim2; ylimits=ylim(ax4); ax4.Children(end-1).YData=ylimits; ax4.Children(end-2).YData=ylimits; + catch + end drawnow limitrate end end diff --git a/matlab/analyse_espicfields.m b/matlab/analyse_espicfields.m index 9831748..dc84c09 100644 --- a/matlab/analyse_espicfields.m +++ b/matlab/analyse_espicfields.m @@ -1,391 +1,391 @@ %% time extent of plotted variables % file='teststablegauss_13_traject2.h5'; % file='unigauss.h5'; % file='Dav_5e14fine_wr+.h5'; % file='stablegauss_8e19fine.h5' % % %file='teststable_Dav1.h5'; % % %close all hidden; % if (~exist('M','var')) % M.file=file; % end % M=load_espic2d(file,M,'all'); t2dlength=size(M.t2d,1) fieldstart=max(1,t2dlength-500);%floor(0.95*t2dlength); fieldend=t2dlength;%23; deltat=t2dlength-fieldstart; %[~, rgridend]=min(abs(M.rgrid-0.045)); rgridend=sum(M.nnr(1:2)); [~, name, ~] = fileparts(M.file); M.displayenergy %% Radial profile try M.displayrprofile([1,fieldstart:fieldend],[],true) catch end %% %fieldstart=2300; fieldend=2500; dens=mean(M.N(:,:,fieldstart:fieldend),3); geomw=M.geomweight(:,:,1); geomw(geomw<0)=0; geomw(geomw>0)=max(max(dens)); dispdens=dens; dispdens(geomw<=0)=NaN; pot=mean(M.pot(:,:,fieldstart:fieldend),3); brillratio=2*dens*M.me./(M.eps_0*M.B'.^2); [R,Z]=meshgrid(M.rgrid,M.zgrid); Rinv=1./R; Rinv(isnan(Rinv))=0; VTHET=mean(M.fluidUTHET(:,:,fieldstart:fieldend),3); omegare=(VTHET.*Rinv'); maxdens=max(max(dens)); % f=figure(); % ax3=gca; % surface(ax3,M.zgrid(1:end-1),M.rgrid(1:end-1),(squeeze(mean(M.Presstens(4,:,:,fieldstart:end),4)))*M.vlight) % xlabel(ax3,'z [m]') % ylabel(ax3,'r [m]') % xlim(ax3,[M.zgrid(1) M.zgrid(end)]) % ylim(ax3,[M.rgrid(1) M.rgrid(50)]) % colormap(ax3,'jet') % c = colorbar(ax3); % c.Label.String= 'thermal v_\theta [m/s]'; % %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]; % %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]) % title(ax3,'Azimuthal velocity') % %set(ax3,'colorscale','log') % grid(ax3, 'on') % view(ax3,2) % f.PaperOrientation='landscape'; % [~, name, ~] = fileparts(M.file); % f.PaperUnits='centimeters'; % papsize=[16 14]; % f.PaperSize=papsize; % print(f,sprintf('%sfluid_thermtheta',name),'-dpdf','-fillpage') %% Density M.displaydensity(fieldstart,fieldend); %% Fieldswide f=figure('Name', sprintf('%s fields',name)); ax1=gca; h=contourf(ax1,M.zgrid,M.rgrid,dispdens,'Displayname','n_e [m^{-3}]'); hold on; %[c1,h]=contour(ax1,M.zgrid,M.rgrid,pot./1e3,5,'m--','ShowText','off','linewidth',1.5,'Displayname','Equipotentials [kV]'); %clabel(c1,h,'Color','white') Blines=zeros(size(M.Athet)); for i=1:length(M.zgrid) Blines(i,:)=M.Athet(i,:).*M.rgrid'; end zindex=floor(length(M.zgrid/2)); levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),8); contour(ax1,M.zgrid,M.rgrid,Blines',real(levels),'r-.','linewidth',1.5,'Displayname','Magnetic field lines'); [~, hContour]=contourf(ax1,M.zgrid,M.rgrid,geomw,2); drawnow; xlim(ax1,[M.zgrid(1) M.zgrid(end)]) ylim(ax1,[M.rgrid(1) M.rgrid(end)]) legend({'n_e [m^{-3}]','Magnetic field lines'},'location','southwest') xlabel(ax1,'z [m]') ylabel(ax1,'r [m]') title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),double(maxdens))) c = colorbar(ax1); c.Label.String= 'n[m^{-3}]'; view(ax1,2) %set(h,'edgecolor','none'); grid on; hFills=hContour.FacePrims; [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' hFills(1).ColorData = uint8([255;255;255;255]); for idx = 2 : numel(hFills) hFills(idx).ColorData(4) = 0; % default=255 end %rectangle('position',[M.zgrid(5) M.rgrid(3) M.zgrid(end-5)-M.zgrid(5) (M.rgrid(rgridend)-M.rgrid(1))],'Edgecolor','white','linestyle','--','linewidth',2) f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[18 10]; f.PaperSize=papsize; print(f,sprintf('%sFieldswide',name),'-dpdf','-fillpage') savefig(f,sprintf('%sFieldswide',name)) %% Fields f=figure('Name', sprintf('%s fields',name)); ax1=gca; h=contourf(ax1,M.zgrid*1000,M.rgrid*1000,dispdens,'Displayname','n_e [m^{-3}]'); hold on; pot(M.geomweight(:,:,1)<0)=0; Blines=zeros(size(M.Athet)); for i=1:length(M.zgrid) Blines(i,:)=M.Athet(i,:).*M.rgrid'; end zindex=floor(length(M.zgrid/2)); %levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),10); levels=linspace(min(min(Blines(:,2:end))),max(max(Blines(:,:))),10); [~,h1]=contour(ax1,M.zgrid*1000,M.rgrid*1000,Blines',real(levels),'r-.','linewidth',1.5,'Displayname','Magnetic field lines'); max(abs(pot(:)))/1e3 [c1,h2]=contour(ax1,M.zgrid*1000,M.rgrid*1000,pot./1e3,'m--','ShowText','on','linewidth',1.2,'Displayname','Equipotentials [kV]'); clabel(c1,h2,'Color','white') %contour(ax1,M.zgrid*1000,M.rgrid*1000,M.geomweight(:,:,1),[0 0],'w-','linewidth',1.5); [c1,hContour]=contourf(ax1,M.zgrid*1000,M.rgrid*1000,geomw); drawnow; xlim(ax1,[M.zgrid(1)*1000 M.zgrid(end)*1000]) if(M.conformgeom) ylim(ax1,[M.rgrid(1)*1000 M.rgrid(rgridend)*1000]) else ylim(ax1,[M.rgrid(1)*1000 M.rgrid(end)*1000]) end legend([h1,h2],{'Magnetic field lines','Equipotentials [kV]'},'location','southwest') xlabel(ax1,'z [mm]') ylabel(ax1,'r [mm]') %title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),double(maxdens))) c = colorbar(ax1); c.Label.String= 'n[m^{-3}]'; view(ax1,2) %set(h,'edgecolor','none'); grid on; hFills=hContour.FacePrims; [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' try hFills(1).ColorData = uint8([150;150;150;255]); for idx = 2 : numel(hFills) hFills(idx).ColorData(4) = 0; % default=255 end catch end f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); if M.maxwellsrce.present rlen=diff(M.maxwellsrce.rlim); zlen=diff(M.maxwellsrce.zlim); rectangle('Position',[M.maxwellsrce.zlim(1) M.maxwellsrce.rlim(1) zlen rlen]*1000,'Edgecolor','g','Linewidth',3,'Linestyle','--') end %rectangle('Position',[-10 73 44 8],'Edgecolor','g','Linewidth',3,'Linestyle','--') f.PaperUnits='centimeters'; papsize=[12 8]; f.PaperSize=papsize; rectangle('Position',[M.zgrid(1) M.rgrid(1) M.zgrid(end)-M.zgrid(1) 0.064-M.rgrid(1)]*1000,'FaceColor',[150 150 150]/255) print(f,sprintf('%sFields',name),'-dpdf','-fillpage') savefig(f,sprintf('%sFields',name)) %% Brillouin f=figure('Name', sprintf('%s Brillouin',name)); ax1=gca; brillratio(geomw<=0)=NaN; h=surface(ax1,M.zgrid,M.rgrid,brillratio); set(h,'edgecolor','none'); xlim(ax1,[M.zgrid(1) M.zgrid(end)]) if(M.conformgeom) ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)]) else ylim(ax1,[M.rgrid(1) M.rgrid(end)]) end xlabel(ax1,'z [m]') ylabel(ax1,'r [m]') title(ax1,'Position') c = colorbar(ax1); c.Label.String= 'Brillouin Ratio'; view(ax1,2) caxis([0 max(1.2,max(brillratio(:)))]) title(ax1,sprintf('Brillouin Ratio t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),double(maxdens))) f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%sfluid_Brillouin',name),'-dpdf','-fillpage') savefig(f,sprintf('%sfluid_Brillouin',name)) %% Omegar f=figure('Name', sprintf('%s Omegar',name)); ax3=gca; omegare(geomw<=0)=NaN; surface(ax3,M.zgrid,M.rgrid,omegare,'edgecolor','None') xlabel(ax3,'z [m]') ylabel(ax3,'r [m]') xlim(ax3,[M.zgrid(1) M.zgrid(end)]) if(M.conformgeom) ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)]) else ylim(ax3,[M.rgrid(1) M.rgrid(end)]) end colormap(ax3,'jet') c = colorbar(ax3); c.Label.String= '|\omega_\theta| [1/s]'; %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]; %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]) title(ax3,sprintf('Azimuthal frequency t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),double(maxdens))) %set(ax3,'colorscale','log') grid(ax3, 'on') view(ax3,2) f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%sfluid_omegar',name),'-dpdf','-fillpage') savefig(f,sprintf('%sfluid_omegar',name)) %% VTHET f=figure('Name', sprintf('%s Vthet',name)); ax3=gca; VTHET(geomw<=0)=NaN; surface(ax3,M.zgrid,M.rgrid,VTHET,'edgecolor','None') xlabel(ax3,'z [m]') ylabel(ax3,'r [m]') xlim(ax3,[M.zgrid(1) M.zgrid(end)]) if(M.conformgeom) ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)]) else ylim(ax3,[M.rgrid(1) M.rgrid(end)]) end colormap(ax3,'jet') c = colorbar(ax3); c.Label.String= '|V_\theta| [m/s]'; %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]; %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]) title(ax3,sprintf('Azimuthal velocity t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),double(maxdens))) %set(ax3,'colorscale','log') grid(ax3, 'on') view(ax3,2) f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%sfluid_vthet',name),'-dpdf','-fillpage') savefig(f,sprintf('%sfluid_vthet',name)) %% vdrift f=figure('Name', sprintf('%s V drift',name)); ax3=gca; vdrift=(mean(M.Ez(:,:,fieldstart:end),3).*M.Br'-mean(M.Er(:,:,fieldstart:end),3).*M.Bz')./(M.B.^2)'; vdrift(geomw<=0)=NaN; v=vdrift; surface(ax3,M.zgrid,M.rgrid,v,'edgecolor','None') xlabel(ax3,'z [m]') ylabel(ax3,'r [m]') xlim(ax3,[M.zgrid(1) M.zgrid(end)]) if(M.conformgeom) ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)]) else ylim(ax3,[M.rgrid(1) M.rgrid(end)]) end colormap(ax3,'jet') c = colorbar(ax3); set(ax3,'zscale','log') set(ax3,'colorscale','log') c.Label.String= 'V_{ExB}'; %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]; %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]) title(ax3,sprintf('V_{ExB} velocity t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),double(maxdens))) %set(ax3,'colorscale','log') grid(ax3, 'on') view(ax3,2) f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%sfluid_EB',name),'-dpdf','-fillpage') savefig(f,sprintf('%sfluid_EB',name)) %% HP -if(M.R.nt>=2) +if(obj.species(1).R.nt>=2) taveragestart=max(floor(0.8*size(M.tpart,1)),2); M.displayHP(taveragestart); end %% VR, VTHET, VZ M.displayVdistribRThetZ() %% Vperp Vpar M.displayVdistribParPer() %% 0D diagnostics BrillouinRatioth=2*M.omepe^2/M.omece^2 BrillouinMax=max(max(brillratio)) NpartsTrapped=mean(M.npart(end-10:end))*M.msim/M.me dblengthiter=fieldstart:fieldend; Pr=squeeze(M.Presstens(1,:,:,dblengthiter)); Pz=squeeze(M.Presstens(6,:,:,dblengthiter)); enddens=M.N(:,:,dblengthiter); meanPr=mean(Pr,3); meanPz=mean(Pz,3); Tr=Pr./enddens; Tz=Pz./enddens; Debye_length=sqrt(abs(min(min(min(mean(Tr(Tr>0)./enddens(Tr>0),3))),min(min(mean(Tz(Tz>0)./enddens(Tz>0),3)))))*M.eps_0/M.qe^2) %% Debye length f=figure('Name', sprintf('%s Debye',name)); subplot(2,1,1) surface(M.zgrid,M.rgrid,sqrt(mean(abs(Tr./enddens*M.eps_0/M.qe^2),3)),'edgecolor','none') colorbar f.PaperOrientation='landscape'; xlabel('z [m]') ylabel('r [m]') c = colorbar; c.Label.String= '\lambda_r [m]'; subplot(2,1,2) surface(M.zgrid,M.rgrid,sqrt(mean(abs(Tz./enddens*M.eps_0/M.qe^2),3)),'edgecolor','none') colorbar f.PaperOrientation='landscape'; xlabel('z [m]') ylabel('r [m]') c = colorbar; c.Label.String= '\lambda_z [m]'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%s_dblength',name),'-dpdf','-fillpage') savefig(f,sprintf('%s_dblength',name)) %% Temperature f=figure('Name', sprintf('%s Temperature',name)); ax1=subplot(2,1,1); title('Radial temperature') surface(M.zgrid,M.rgrid,mean(Tr,3)/M.qe,'edgecolor','none') colormap('jet') c = colorbar; c.Label.String= 'T_{er} [eV]'; templimits1=caxis; ax2=subplot(2,1,2); title('Axial tempreature') surface(M.zgrid,M.rgrid,mean(Tz,3)/M.qe,'edgecolor','none') colormap('jet') c = colorbar; c.Label.String= 'T_{ez} [eV]'; maxtemp=max([templimits1,caxis]); caxis(ax1,[0 maxtemp]) caxis(ax2,[0 maxtemp]) f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%s_temp',name),'-dpdf','-fillpage') savefig(f,sprintf('%s_temp',name)) %% Wells M.display2Dpotentialwell(t2dlength-1); M.display2Dpotentialwell(t2dlength-1,false); diff --git a/matlab/displayParPer.m b/matlab/displayParPer.m index 25191e1..e993326 100644 --- a/matlab/displayParPer.m +++ b/matlab/displayParPer.m @@ -1,68 +1,70 @@ %% -timesteppart=length(M.tpart); +timesteppart=5; fig=figure; ax1=gca; -timestepN=find(M.tpart(end)==M.t2d,1); +timestepN=find(M.tpart(timesteppart)==M.t2d,1); Ndistrib=M.N(:,:,timestepN); h=contourf(ax1,M.zgrid,M.rgrid,Ndistrib); %set(h, 'EdgeColor', 'none'); hold on [r,z]=find(Ndistrib~=0); xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))]) ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))]) xlabel(ax1,'Z [m]') ylabel(ax1,'R [m]') c = colorbar(ax1); c.Label.String= 'n [m^{-3}]'; view(ax1,2) +hold on +contour(M.zgrid,M.rgrid,M.geomweight(:,:,1),[0 0],'r-') %set(ax1,'colorscale','log') [x,y]=ginput(1); Zindex=find(x>M.zgrid,1,'last'); Rindex=find(y>M.rgrid,1,'last'); % Rindex=155; % Zindex=344; % Rindex=123; % Zindex=658; gcs=false; f=figure(); ax1=subplot(1,2,1); ax2=subplot(1,2,2); %[p,maxnb,c]=M.displayPhaseSpace(1:2,Rindex,Zindex,'',sprintf('t=%1.3g [s]',M.tpart(1)),ax1); % ax1=subplot(1,2,1); % hold(ax1,'off') c=cell(2,1); c{1}=linspace(-2e7,2e7,21); c{2}=linspace(-0e6,3e7,51); [p,maxnb,c]=M.displayPhaseSpace('parper',2,Rindex,Zindex,'',sprintf('t=%1.3g [s]',M.tpart(2)),ax1,-1,c,gcs); M.displayPhaseSpace('parper',length(M.tpart)+(0),Rindex,Zindex,'',sprintf('t=%1.3g [s]',M.tpart(end)),ax2,-1,c,gcs); xlimits=xlim(ax1); xlimits=[min([xlimits,xlim(ax2)]) max([xlimits,xlim(ax2)])]; ylimits=ylim(ax1); ylimits=[min([ylimits,ylim(ax2)]) max([ylimits,ylim(ax2)])]; xlim([ax1,ax2],xlimits) ylim([ax1,ax2],ylimits) climitsmax=max([caxis(ax1),caxis(ax2)]); caxis(ax1,[0 climitsmax]); caxis(ax2,[0 climitsmax]); xtickformat(ax1,'%.3g') ytickformat(ax1,'%.3g') ax1.YAxis.Exponent = 0; ax1.XAxis.Exponent = 0; legend(ax2,'off'); ax1.Children(5).LineStyle='none'; xtickformat(ax2,'%.3g') ytickformat(ax2,'%.3g') ax2.YAxis.Exponent = 0; ax2.XAxis.Exponent = 0; sgtitle(sprintf('r=%1.2f [mm] z=%1.2f [mm] \\Delta\\phi=%1.1f[kV] R=%1.1f',M.rgrid(Rindex)*1e3,M.zgrid(Zindex)*1e3,(M.potout-M.potinn)/1e3,M.Rcurv)) M.savegraph(f,sprintf('%s/%s_phasespaceR%dZ%dbegendnogcs',M.folder,M.name,Rindex,Zindex),[16,10]); \ No newline at end of file diff --git a/matlab/distribution_function.m b/matlab/distribution_function.m index 8fd2bb8..c7ce5b6 100644 --- a/matlab/distribution_function.m +++ b/matlab/distribution_function.m @@ -1,374 +1,375 @@ %% Show the particles velocity histogram at the position set using ginput % The histogram is compared between timesteppart and end -timesteppart=length(M.tpart); +timesteppart=8; fig=figure; ax1=gca; -timestepN=find(M.tpart(end)==M.t2d,1); +timestepN=find(M.tpart(timesteppart)==M.t2d,1); Ndistrib=M.N(:,:,timestepN); h=contourf(ax1,M.zgrid,M.rgrid,Ndistrib); %set(h, 'EdgeColor', 'none'); hold on +contour(M.zgrid,M.rgrid,M.geomweight(:,:,1),[0 0],'r-') [r,z]=find(Ndistrib~=0); xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))]) ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))]) xlabel(ax1,'Z [m]') ylabel(ax1,'R [m]') c = colorbar(ax1); c.Label.String= 'n [m^{-3}]'; view(ax1,2) %set(ax1,'colorscale','log') % [x,y]=ginput(1); % Zindex=find(x>M.zgrid,1,'last'); % Rindex=find(y>M.rgrid,1,'last'); % Rindex=155; % Zindex=344; Z=(M.zgrid(Zindex)); R=(M.rgrid(Rindex)); plot(Z,R,'rx','Markersize',12); % Rindex=16; Zindex=64; % Rindex=28; Zindex=floor(size(M.zgrid,1)/2)+1; -nbp=min(obj.species(1).R.nparts,M.nbparts(1)); -Rp=obj.species(1).R(1:nbp,1,false); -Zp=obj.species(1).Z(1:nbp,1,false); +nbp=min(M.species(1).R.nparts,M.species(1).nbparts(1)); +Rp=M.species(1).R(1:nbp,1,false); +Zp=M.species(1).Z(1:nbp,1,false); Indices=Rp>=M.rgrid(Rindex) & Rp=M.zgrid(Zindex) & Zp=M.rgrid(Rindex) & Rend=M.zgrid(Zindex) & Zend1) h1=histogram(Vr,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); %h1=histfit(ax1,Vr); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); end hold on h1=histogram(Vrend,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); %h1=histfit(ax1,Vrend); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); ylabel('counts') xlabel('\beta_r') grid on ax2=subplot(1,3,2); binwidth=abs(max(Vthetend)-min(Vthetend))/sqrt(length(Vthetend)); if(length(Vthet)>1) h1=histogram(Vthet(:,1),'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); end %h1=histfit(ax2,Vthet(:,1)); hold on h1=histogram(Vthetend,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); %h1=histfit(ax2,Vthetend); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); ylabel('counts') xlabel('\beta_\theta') legend(ax2,'location','northoutside','orientation','vertical') grid on ax3=subplot(1,3,3); %h1=histogram(Vz(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); binwidth=abs(max(Vzend)-min(Vzend))/sqrt(length(Vzend)); if(length(Vz)>1) h1=histogram(ax3,Vz(:,1),'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); end hold on %h1=histogram(Vzend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); h1=histogram(ax3,Vzend,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); ylabel('counts') xlabel('\beta_z') grid on f=gcf; sgtitle(sprintf('R=%1.2e[m] Z=%1.2e[m] dt=%1.2e[ns]',M.rgrid(Rindex+1),M.zgrid(Zindex+1),M.dt*1e9)) f.PaperOrientation='landscape'; f.PaperUnits='centimeters'; papsize=[16 10]; f.PaperSize=papsize; [~, name, ~] = fileparts(M.file); print(f,sprintf('%sParts_V_RZ',name),'-dpdf','-fillpage') savefig(f,sprintf('%sParts_V_RZ',name)) %% Shows the particles phase space at position (rindex, zindex) % and times tinit and timesteppart timeinit=1; timestepNinit=find(M.tpart(timeinit)==M.t2d); timesteppart=length(M.tpart); timestepNend=find(M.tpart(timesteppart)==M.t2d); -Rp=obj.species(1).R(:,timeinit,false); -Zp=obj.species(1).Z(:,timeinit,false); +Rp=M.species(1).R(:,timeinit,false); +Zp=M.species(1).Z(:,timeinit,false); Indices=Rp>=M.rgrid(Rindex) & Rp=M.zgrid(Zindex) & Zp=M.rgrid(Rindex) & Rend=M.zgrid(Zindex) & Zend=M.rAthet(:,zleftlim),1,'last'); rleft=sqrt(levels(1)/(M.rAthet(rAthetposend,1)/M.rgrid(rAthetposend)^2)); rposleft=find(M.rgrid>=rleft,1,'first'); n2=mean(mean(N(rposleft,[1 end],:),3)); vpar2=(dN/V/n2)^2; if n2==0 vpar2=0; end b=M.rgrid(end); a=M.rgrid(1); vd1=((M.potout-M.potinn)/M.rgrid(Rindex))/M.Bz(Zindex,Rindex)/log(b/a); vd1=-M.Er(Rindex,Zindex,1)/M.Bz(Zindex,Rindex); %vd1=-M.Er(Rindex,Zindex,timestepNend)/M.Bz(Zindex,Rindex); vinit=sqrt(M.kb*10000/M.me); %vinit=sqrt(M.qe*80/M.me) %vinit=sqrt(120*M.qe/M.me) %vd2=((M.potout-M.potinn)/M.rgrid(25))/M.Bz(1,25)/log(b/a); [Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,floor(M.nz/2)+1)); Rcurv=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex); deltaphicomp=-0.5*M.me/M.qe*(vpar2-(vd1+vinit)^2*(1-M.Rcurv)) ratioparper=vpar2/((vd1+vinit)^2*(1-M.Rcurv)) Ndistrib=M.N(:,:,timestepNend); model=M.potentialwellmodel(timestepNend); z=model.z; r=model.r; pot=model.pot; rathet=model.rathet; Zeval=[M.zgrid(zleftlim) M.zgrid(Zindex)]; Psieval=M.rAthet(Rindex,Zindex)*[1 1]; phis=griddata(Zmesh,M.rAthet,M.pot(:,:,end),Zeval,Psieval,'natural'); deltaphi=-diff(phis) %deltaphi=-215; R=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex); % if vper is above line then electron is kept axially vper=sqrt(2*M.qe/M.me*deltaphi/(R-1)+vpar.^2/(R-1))/M.vlight; vpar=vpar/M.vlight; vpar=vpar(real(vper)~=0); vper=vper(real(vper)~=0); xlimits=xlim; ylimits=ylim; if(length(vper)>0) p(3)=plot(vpar,vper,'b-','displayname','Loss parabolla simul'); hold on plot(-vpar,vper,'b-') end vpar=linspace(1,M.vlight,1000); vper=sqrt(-2*M.qe/M.me*deltaphicomp/(R-1)+vpar.^2/(R-1))/M.vlight; vpar=vpar/M.vlight; vpar=vpar(real(vper)~=0); vper=vper(real(vper)~=0); if(length(vper)>0) p(4)=plot(vpar,vper,'k--','displayname','Loss parabolla prediction'); hold on plot(-vpar,vper,'k--') end xlim(xlimits) ylim(ylimits) legend(p) axis equal % subplot(2,2,2) % scatter(Zp,Vpar,'.','displayname','Init') % hold on % scatter(Zend,Vparend,'.','displayname','End') % xlabel('z [m]') % ylabel('\beta_{par}') % % subplot(2,2,3) % scatter(Vperp,Rp,'.','displayname','Init') % hold on % scatter(Vperpend,Rend,'.','displayname','End') % xlabel('\beta_\perp') % ylabel('R [m]') ax1=subplot(1,2,2); h=contourf(ax1,M.zgrid,M.rgrid,Ndistrib); hold on [r,z]=find(Ndistrib~=0); xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))]) ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))]) xlabel(ax1,'Z [m]') ylabel(ax1,'R [m]') c = colorbar(ax1); c.Label.String= 'n [m^{-3}]'; Zx=(M.zgrid(Zindex)); Rx=(M.rgrid(Rindex)); plot(ax1,Zx,Rx,'rx','Markersize',12); sgtitle(sprintf('r=%1.3g [m] z=%1.3g [m] \\Delta\\phi=%1.3g[kV] R=%1.3g',M.rgrid(Rindex),M.zgrid(Zindex),(M.potout-M.potinn),M.Rcurv)) M.savegraph(f,sprintf('%s/%s_phasespaceR%dZ%dpos',M.folder,M.name,Rindex,Zindex),[16,12]) %% f=figure; p(1)=scatter(Vpar,Vperp,'.','displayname','Initial'); hold on p(2)=scatter(Vparend,Vperpend,'.','displayname','Final'); xlabel('v_{par}/c') ylabel('v_\perp/c') zleftlim=1; vpar=linspace(1,M.vlight,1000); dN=1e13*M.weight; dN=0; rp=9.85e-3; rm=7e-3; V=2*pi*(rp^2-rm^2); N=M.N(:,:,timestepNend); levels=M.rAthet(Rindex,Zindex)*[1 1]; rAthetposend=find(levels(1)>=M.rAthet(:,zleftlim),1,'last'); rleft=sqrt(levels(1)/(M.rAthet(rAthetposend,1)/M.rgrid(rAthetposend)^2)); rposleft=find(M.rgrid>=rleft,1,'first'); n2=mean(mean(N(rposleft,[1 end],:),3)); vpar2=(dN/V/n2)^2; if n2==0 vpar2=0; end b=M.rgrid(end); a=M.rgrid(1); vd1=((M.potout-M.potinn)/M.rgrid(Rindex))/M.Bz(Zindex,Rindex)/log(b/a); vd1=-M.Er(Rindex,Zindex,1)/M.Bz(Zindex,Rindex); %vd1=-M.Er(Rindex,Zindex,timestepNend)/M.Bz(Zindex,Rindex); vinit=sqrt(M.kb*10000/M.me); %vinit=sqrt(M.qe*80/M.me) %vinit=sqrt(120*M.qe/M.me) %vd2=((M.potout-M.potinn)/M.rgrid(25))/M.Bz(1,25)/log(b/a); [Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,floor(M.nz/2)+1)); % Find the R magnetic ratio between local position and left Rcurv=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex); deltaphicomp=-0.5*M.me/M.qe*(vpar2-(vd1+vinit)^2*(1-M.Rcurv)) ratioparper=vpar2/((vd1+vinit)^2*(1-M.Rcurv)) Ndistrib=M.N(:,:,timestepNend); model=M.potentialwellmodel(timestepNend); z=model.z; r=model.r; pot=model.pot; rathet=model.rathet; Zeval=[M.zgrid(zleftlim) M.zgrid(Zindex)]; Psieval=M.rAthet(Rindex,Zindex)*[1 1]; phis=griddata(Zmesh,M.rAthet,M.pot(:,:,end),Zeval,Psieval,'natural'); % calculate the potential difference between left and local deltaphi=-diff(phis) %deltaphi=-215; R=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex); vper=sqrt(2*M.qe/M.me*deltaphi/(R-1)+vpar.^2/(R-1))/M.vlight; vpar=vpar/M.vlight; vpar=vpar(real(vper)~=0); vper=vper(real(vper)~=0); xlimits=xlim; ylimits=ylim; if(length(vper)>0) p(3)=plot(vpar,vper,'b-','displayname','Simulation'); hold on plot(-vpar,vper,'b-') end vpar=linspace(1,M.vlight,1000); vper=sqrt(-2*M.qe/M.me*deltaphicomp/(R-1)+vpar.^2/(R-1))/M.vlight; vpar=vpar/M.vlight; vpar=vpar(real(vper)~=0); vper=vper(real(vper)~=0); if(length(vper)>0) p(4)=plot(vpar,vper,'k--','displayname','Prediction'); hold on plot(-vpar,vper,'k--') end xlim(xlimits) ylim(ylimits) legend(p) axis equal legend('location','northeast') title(sprintf('r=%1.3g [m] z=%1.3g [m] \\Delta\\phi=%1.3g[kV] R=%1.3g',M.rgrid(Rindex),M.zgrid(Zindex),(M.potout-M.potinn),M.Rcurv)) M.savegraph(f,sprintf('%s/%s_phasespaceR%dZ%d',M.folder,M.name,Rindex,Zindex),[12,8]) \ No newline at end of file diff --git a/matlab/guiParts.m b/matlab/guiParts.m index b61b149..5959017 100644 --- a/matlab/guiParts.m +++ b/matlab/guiParts.m @@ -1,237 +1,237 @@ -function dispespicParts(M,fieldstep,parper) +function guiParts(M,fieldstep,parper) %dispespicParts Show on an interactive plot the phase space of the given specie in M % Detailed explanation goes here if nargin<2 fieldstep=1; end if nargin<3 parper=false; end f2=uifigure('Name','Particles data'); mf2=uipanel(f2,'Position',[5 50 f2.Position(3)-10 f2.Position(4)-55]); mf2.AutoResizeChildren='off'; m2=uipanel(f2,'Position',[5 5 f2.Position(3)-10 40]); sgtitle(mf2,sprintf('step=%d t=%0.5e s',fieldstep*M.it2,M.tpart(fieldstep))) sld2 = uislider(m2,'Position',[10 30 0.6*m2.Position(3) 3]); sld2.Value=fieldstep; sld2.Limits=[1 length(M.tpart)]; edt2 = uieditfield(m2,'numeric','Limits',[1 length(M.tpart)],'Value',1); edt2.Position=[sld2.Position(1)+sld2.Position(3)+25 5 40 20]; edt2.RoundFractionalValues='on'; Printbt2=uibutton(m2,'Position',[edt2.Position(1)+edt2.Position(3)+10 5 40 20],'Text', 'Save'); sld2.ValueChangingFcn={@updatefigpartdata,edt2,mf2}; edt2.ValueChangedFcn={@updatefigpartdata,sld2,mf2}; Printbt2.ButtonPushedFcn={@plotPartButtonPushed}; set(f2,'KeyPressFcn',{ @onKeyDown,sld2,edt2,mf2}) Plotfennecspartdata(mf2,M,fieldstep); function Plotfennecspartdata(fig,M,fieldstep) %Plotfennecsgriddata Plot the 2d data of fennecs at time step fieldstep %sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it2,M.tpart(fieldstep))) sgtitle(fig,sprintf('t=%0.5e s nbparts=%d',M.tpart(fieldstep),M.nbparts(fieldstep))) vmax=4e7; Z=M.Z(1:min(M.nbparts(fieldstep),M.Z.nparts),fieldstep); R=M.R(1:min(M.nbparts(fieldstep),M.R.nparts),fieldstep); if parper vpar=M.Vpar(1:min(M.nbparts(fieldstep),M.Z.nparts),fieldstep); vper=M.Vperp(1:min(M.nbparts(fieldstep),M.Z.nparts),fieldstep); else vpar=M.VZ(1:min(M.nbparts(fieldstep),M.Z.nparts),fieldstep); vper=M.VR(1:min(M.nbparts(fieldstep),M.Z.nparts),fieldstep); end ax1=subplot(3,2,1,'Parent',fig); if isa(M,'h5parts') contour(ax1,M.zgrid,M.rgrid,M.parent.geomweight(:,:,1),[0 0],'r--') rindex=M.rindex; zindex=M.zindex; else contour(ax1,M.zgrid,M.rgrid,M.geomweight(:,:,1),[0 0],'r--') rindex=[1 length(M.rgrid)]; zindex=[1 length(M.zgrid)]; end hold(ax1, 'on') rectangle(ax1,'Position',[M.zgrid(zindex(1)) M.rgrid(rindex(1))... M.zgrid(zindex(end))-M.zgrid(zindex(1)) M.rgrid(rindex(end))-M.rgrid(rindex(1))],'linestyle','--') plot(ax1,Z,R,'.'); if isa(M,'h5parts') xlim(ax1,[M.zgrid(zindex(1)) M.zgrid(zindex(end))]) ylim(ax1,[M.rgrid(rindex(1)) M.rgrid(rindex(end))]) end xlabel(ax1,'Z [m]') ylabel(ax1,'R [m]') title(ax1,'Position') grid(ax1, 'on') ax2=subplot(3,2,2,'Parent',fig); plot(ax2,vpar,vper,'.'); if parper xlabel(ax2,'V_{par} [m/s]') ylabel(ax2,'V_{perp} [m/s]') else xlabel(ax2,'V_Z[m/s]') ylabel(ax2,'V_R [m/s]') end ylim(ax2,vmax*[-1 1]) xlim(ax2,vmax*[-1 1]) axis(ax2, 'equal') grid(ax2, 'on') ax3=subplot(3,2,3,'Parent',fig); plot(ax3,Z,vpar,'.'); xlabel(ax3,'Z [m]') if parper ylabel(ax3,'V_{par} [m/s]') else ylabel(ax3,'V_Z[m/s]') end grid(ax3, 'on') xlim(ax3,[M.zgrid(zindex(1)) M.zgrid(zindex(end))]) ylim(ax3,vmax*[-1 1]) ax4=subplot(3,2,4,'Parent',fig); plot(ax4,R,vpar,'.') if parper ylabel(ax4,'V_{par} [m/s]') else ylabel(ax4,'V_Z[m/s]') end xlabel(ax4,'R [m]') xlim(ax4,[M.rgrid(rindex(1)) M.rgrid(rindex(end))]) ylim(ax4,vmax*[-1 1]) grid(ax4, 'on') ax5=subplot(3,2,5,'Parent',fig); plot(ax5,Z,vper,'.'); xlim(ax5,[M.zgrid(zindex(1)) M.zgrid(zindex(end))]) xlabel(ax5,'Z [m]') if parper ylabel(ax5,'V_{perp} [m/s]') else ylabel(ax5,'V_R [m/s]') end grid(ax5, 'on') ylim(ax5,vmax*[-1 1]) ax6=subplot(3,2,6,'Parent',fig); plot(ax6,R,vper,'.'); xlim(ax6,[M.rgrid(rindex(1)) M.rgrid(rindex(end))]) if parper ylabel(ax6,'V_{perp} [m/s]') else ylabel(ax6,'V_R [m/s]') end xlabel(ax6,'R [m]') ylim(ax6,vmax*[-1 1]) grid(ax6, 'on') end function onKeyDown(src,event,slider,editfield, fig) direction=0; if strcmp(event.Key,'leftarrow') direction=-1; elseif strcmp(event.Key,'rightarrow') direction=+1; elseif strcmp(event.Key,'uparrow') direction=+10; elseif strcmp(event.Key,'downarrow') direction=-10; end if(direction~=0) currval=slider.Value; slider.Value=max(slider.Limits(1),min(currval+direction,slider.Limits(2))); updatefigpartdata(slider, event, editfield ,fig) end end function updatefigpartdata(control, event, Othercontrol, fig) if strcmp(event.EventName,'ValueChanged') fieldstep=floor(control.Value); control.Value=fieldstep; elseif strcmp(event.EventName,'KeyPress') fieldstep=floor(control.Value); control.Value=fieldstep; else fieldstep=floor(event.Value); end Othercontrol.Value=fieldstep; sgtitle(fig,sprintf('t=%0.5e s nbparts=%03d',M.tpart(fieldstep),M.nbparts(fieldstep))) % if isa(M,'fennecshdf5') % [~,t0dstep]=min(abs(M.tpart(fieldstep)-M.t0d)); % else t0dstep=fieldstep; % end Z=M.Z(1:min(M.nbparts(t0dstep),M.Z.nparts),fieldstep); R=M.R(1:min(M.nbparts(t0dstep),M.R.nparts),fieldstep); if parper vpar=M.Vpar(1:min(M.nbparts(t0dstep),M.Z.nparts),fieldstep); vper=M.Vperp(1:min(M.nbparts(t0dstep),M.Z.nparts),fieldstep); else vpar=M.VZ(1:min(M.nbparts(t0dstep),M.Z.nparts),fieldstep); vper=M.VR(1:min(M.nbparts(t0dstep),M.Z.nparts),fieldstep); end %% update Position plot ax1=fig.Children(end); ax1.Children(1).XData=Z; ax1.Children(1).YData=R; view(ax1,2) %% update VPAR VPERP fig.Children(end-1).Children(1).XData=vpar; fig.Children(end-1).Children(1).YData=vper; axis(fig.Children(end-1),'equal') %% update Z VZ fig.Children(end-2).Children(1).XData=Z; fig.Children(end-2).Children(1).YData=vpar; %% update VZ VR fig.Children(end-3).Children(1).XData=R; fig.Children(end-3).Children(1).YData=vpar; %% update R VR fig.Children(end-4).Children(1).XData=Z; fig.Children(end-4).Children(1).YData=vper; %% update Z VR fig.Children(end-5).Children(1).XData=R; fig.Children(end-5).Children(1).YData=vper; %drawnow limitrate end function plotPartButtonPushed(btn,ax) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here f=figure(); Plotfennecspartdata(f,M,floor(sld2.Value)); f.PaperOrientation='portrait'; [~, name, ~] = fileparts(M.file); print(f,sprintf('%sParts%d',name,sld2.Value),'-dpdf','-fillpage') end end diff --git a/matlab/rbfindercoax.m b/matlab/rbfindercoax.m index b913f36..d9caf8b 100644 --- a/matlab/rbfindercoax.m +++ b/matlab/rbfindercoax.m @@ -1,147 +1,147 @@ %% wr+ H0=3.2e-14; P0=8.66e-25; r0=0.005; ne=5e14; B0=0.21; Rcurv=1.5; L=0.48; z=0.0; -consideredphia=-60000; +consideredphia=0; qe=1.60217662E-19; me=9.109383E-31; eps0=8.85418781762e-12; phia=-linspace(0,90000,80); phib=0; -b=0.16; +b=0.06; a=0.001; %% simulated r_sim= [5.8e-3, 5.4e-3];%, 0.0063]; rplussim=[1.11e-2, 1.20e-2];%, 0.042]; phiasim= [-45000, -30000];%, 3.49e17]; rb_phi=zeros(size(phia)); rbplusphi=zeros(size(phia)); i=1; minus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0); format long; rb_phi(i)=fzero(minus,r0); %plus=@(r) qe^2*ne(i)*rb_(i)/(eps0*H0)*(r-rb_(i)-rb_(i)*log(r/rb_(i)))+1-1/(2*me*H0)*(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2; plus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0); rbplusphi(i)=fzero(plus,5*rb_phi(i)); remainderphi=ones(size(phia)); remainderphi(i)=plus(rbplusphi(i)); for i=2:length(phia) try minus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0); rb_phi(i)=fzero(minus,rb_phi(i-1)); rb_t=rb_phi(i); plus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0); rbplusphi(i)=fzero(plus,rbplusphi(i-1)); remainderphi(i)=plus(rbplusphi(i)); catch end end f=figure; %subplot(2,1,1) plot(phia(remainderphi~=1)/1e3,rb_phi(remainderphi~=1)*1e3,'x-','displayname','r_b^-') hold on plot(phia(remainderphi~=1)/1e3,rbplusphi(remainderphi~=1)*1e3,'x-','displayname','r_b^+') plot(phiasim/1e3,rplussim*1e3,'+--','displayname','r_b^+ simulated','linewidth',1.5) plot(phiasim/1e3,r_sim*1e3,'+--','displayname','r_b^- simulated','linewidth',1.5) legend ylabel('r_b [mm]') -ylim([0 160]) +ylim([0 b*1000]) xlabel('\Phi_a [kV]') % subplot(2,1,2) % semilogx(ne(remainder~=1),remainder(remainder~=1),'x--') % ylabel('remainder [a.u.]') % xlabel('n_e(r_b^-) [m^{-3}]') f.PaperOrientation='landscape'; f.PaperUnits='centimeters'; papsize=[10 10]; f.PaperSize=papsize; print(f,'rbne','-dpdf','-fillpage') i=find(consideredphia>=phia,1,'first'); -z=linspace(0,L/2,6000); +z=linspace(0,L/2,10); rb_z=zeros(size(z)); rbplusz=rb_z; minus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z(1)/L)))^2/(2*me*H0); rb_z(1)=fzero(minus,rb_phi(i)); rb_t=rb_z(1); z_t=z(1); plus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z(1)/L)))^2/(2*me*H0); rbplusz(1)=fzero(plus,rbplusphi(i)); remainder(1)=plus(rbplusz(1)); for j=2:length(z) minus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z(j)/L)))^2/(2*me*H0); format long; try rb_z(j)=fzero(minus,rb_z(j-1)); z_t=z(j); plus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z_t/L)))^2/(2*me*H0); rbplusz(j)=fzero(plus,rbplusz(j-1)); remainder(j)=plus(rbplusz(j)); catch end end figure() z=[-flip(z(2:end)) z]; rb_z=[flip(rb_z(2:end)) rb_z]; rbplusz=[flip(rbplusz(2:end)) rbplusz]; plot(z(rb_z~=0),rb_z(rb_z~=0),'b') hold on plot(z(rb_z~=0),rbplusz(rb_z~=0),'r') ylabel('r_b [m]') -ylim([0 0.16]) +ylim([0 b]) xlim([-L/2 L/2]) xlabel('z [m]') title(sprintf('\\phi_a=%5.2f kV',phia(i)/1e3)) ind=rbplusz~=0 & rb_z~=0& ~isnan(rbplusz) & ~isnan(rb_z); ra=rb_z(ind); rb=rbplusz(ind); zvol=z(ind); ra=(ra(1:end-1)+ra(2:end))/2; rb=(rb(1:end-1)+rb(2:end))/2; VOL=2*pi*diff(zvol).*min(ra).*(rb-ra); nplasma=2116800 Nbesl=ne*VOL; Nbe=sum(Nbesl) Nmacrosl=floor(Nbesl/Nbe*nplasma); Nmacro=sum(Nmacrosl) remain=nplasma-Nmacro Nmacrosl(ceil(length(Nmacrosl)/2-remain/2):floor(length(Nmacrosl)/2+remain/2))=Nmacrosl(ceil(length(Nmacrosl)/2-remain/2):floor(length(Nmacrosl)/2+remain/2))+1; Nmacrosl(floor(length(Nmacrosl)/2))=Nmacrosl(floor(length(Nmacrosl)/2))+nplasma-sum(Nmacrosl); Nmacro=sum(Nmacrosl) remain=nplasma-Nmacro msim=me*Nbe/nplasma qsim=-qe*Nbe/nplasma datas.ra=ra; datas.rb=rb; datas.z=zvol; datas.npartsslice=Nmacrosl; datas.m=me; datas.weight=Nbe/nplasma; datas.q=-qe; datas.H0=H0; datas.P0=P0; datas.temperature=10000; datas.velocitytype=2; -datas.radialtype=1; \ No newline at end of file +datas.radialtype=2; \ No newline at end of file diff --git a/matlab/sscalling.m b/matlab/sscalling.m new file mode 100644 index 0000000..0df9edc --- /dev/null +++ b/matlab/sscalling.m @@ -0,0 +1,46 @@ +%% Strong scaling for a collisional system with 10M electrons using FENNECS +% the geometry is a coaxial configuration of constant radius + +Ncores = [ 1 2 4 8 16];% 36 72]; +time = [ (10890+10150)/2 6402 3799 2443 2128];% 1665 1155]; + +f=figure; +plot(Ncores, time(1)*Ncores(1)./time,'r-','displayname','Measured','linewidth',1.5) +hold on +plot(Ncores,Ncores,'k--','displayname','Ideal','linewidth',1.5) +ylabel('t_1/t_N') +yyaxis right +plot(Ncores,time(1).*Ncores(1)./(time.*Ncores),'bx:','linewidth',1.5,'markersize',12,'displayname','Efficiency') +xlabel('N_{cores}') +ylabel('Efficiency') +set(gca,'fontsize',12) +grid on +legend('location','southeast') +ax = gca; +ax.YAxis(1).Color = 'k'; +ax.YAxis(2).Color = 'b'; +ylim([0 1]) +M.savegraph(f,'sscalingOMP'); + +%% MPI with 4 cores per task + +Ntasks = [ 1 2 4 8]; +time = [3799 1533 976 654]; + +f=figure; +plot(Ntasks, time(1)./time,'r-','displayname','Measured','linewidth',1.5) +hold on +ylabel('t_1/t_N') +plot(Ntasks,Ntasks,'k--','displayname','Ideal','linewidth',1.5) +yyaxis right +plot(Ntasks,time(1).*Ntasks(1)./(time.*Ntasks),'bx:','linewidth',1.5,'markersize',12,'displayname','Efficiency') +xlabel('N_{tasks}') +ylabel('Efficiency') +set(gca,'fontsize',12) +grid on +legend('location','southeast') +ylim([0 2]) +ax = gca; +ax.YAxis(1).Color = 'k'; +ax.YAxis(2).Color = 'b'; +M.savegraph(f,'sscalingMPI');