# Plasmacyl

This example is compatible with 1.0.6dev (22 February 2017).

The next example simulation is a modification of the plasmacyl test
problem (*tests/plasmacyl.cpp* in the IBSimu distribution). In
this problem we are extracting a 60 mA/cm2 He+ beam from plasma. Let's
first build a source file *plasmacyl.cpp* and add the necessary
header files:

#include "epot_bicgstabsolver.hpp" #include "epot_problem.hpp" #include "particledatabase.hpp" #include "geometry.hpp" #include "func_solid.hpp" #include "epot_efield.hpp" #include "meshvectorfield.hpp" #include "ibsimu.hpp" #include "error.hpp" #include "particlediagplotter.hpp" #include "geomplotter.hpp" #include "config.hpp" #ifdef GTK3 #include "gtkplotter.hpp" #endif

The *gtkplotter.hpp* file is included only conditionally if
GTK3 support was compiled in the IBSimu library.

The main function is pretty much the same as with the previous example. The only difference is that we need to pass the argc and argv parameters to the simulation function simu() to be able to use the interactive plotter:

int main( int argc, char **argv ) { try { ibsimu.set_message_threshold( MSG_VERBOSE, 1 ); ibsimu.set_thread_count( 4 ); simu( &argc, &argv ); } catch ( Error e ) { e.print_error_message( ibsimu.message( 0 ) ); exit( 1 ); } return( 0 ); }

The geometry defining functions are needed for two electrodes in the problem:

bool solid1( double x, double y, double z ) { return( x <= 0.00187 && y >= 0.00054 && y >= 2.28*x - 0.0010 && (x >= 0.00054 || y >= 0.0015) ); } bool solid2( double x, double y, double z ) { return( x >= 0.0095 && y >= 0.0023333 && y >= 0.01283 - x ); }

And the simulation function itself is defined as:

void simu( int *argc, char ***argv )

The simulation geometry is 12x7 mm in size and we choose a mesh size of 0.05 mm, which is enough to model all the necessary details of the system. In this case the geometry has a cylindrical symmetry (round holes in electrodes). Therefore the boundary 3 isn't actually Neumann boundary, but a special non-existent boundary at radius zero. It is still marked as Neumann as Dirichlet boundary would force the axis to fixed potential (infinitely thin wire in the center). The first boundary is the plasma of ion source. That surface should be set to use the Neumann boundary condition. Even if the simulation should work in principle if you put the ion source to non-zero potential, I will still recommend always redefining the electrode potentials so that the ion source is at zero potential. In this way the calculation of the plasma is most accurate and minimal amount of roundoff errors happens with the floating point operations of the plasma calculation exponentials. The geometry definition:

Geometry geom( MODE_CYL, Int3D(241,141,1), Vec3D(0,0,0), 0.00005 ); Solid *s1 = new FuncSolid( solid1 ); geom.set_solid( 7, s1 ); Solid *s2 = new FuncSolid( solid2 ); geom.set_solid( 8, s2 ); geom.set_boundary( 1, Bound(BOUND_NEUMANN, 0.0 ) ); geom.set_boundary( 2, Bound(BOUND_DIRICHLET, -12.0e3) ); geom.set_boundary( 3, Bound(BOUND_NEUMANN, 0.0) ); geom.set_boundary( 4, Bound(BOUND_NEUMANN, 0.0) ); geom.set_boundary( 7, Bound(BOUND_DIRICHLET, 0.0) ); geom.set_boundary( 8, Bound(BOUND_DIRICHLET, -12.0e3) ); geom.build_mesh();

Next we define the Poisson problem to be solved. For the first
round in the Vlasov iteration we need to give a starting guess for the
plasma boundary. This is done with *set_initial_plasma()* with
first argument defining the plasma potential (5 volts) and the second
argument giving a pointer to a class defining the plasma volume. In
case of simple plasma volume at x<0.55 mm we can use a pre-defined
convenience class *InitialPlasma*:

EpotBiCGSTABSolver solver( geom ); InitialPlasma initp( AXIS_X, 0.55e-3 ); solver.set_initial_plasma( 5.0, &initp );

Next we define the needed objects like in the first example:

EpotField epot( geom ); MeshScalarField scharge( geom ); MeshVectorField bfield; EpotEfield efield( epot ); efield_extrpl_e efldextrpl[6] = { EFIELD_EXTRAPOLATE, EFIELD_EXTRAPOLATE, EFIELD_SYMMETRIC_POTENTIAL, EFIELD_EXTRAPOLATE, EFIELD_EXTRAPOLATE, EFIELD_EXTRAPOLATE }; efield.set_extrapolation( efldextrpl ); ParticleDataBaseCyl pdb( geom ); bool pmirror[6] = { false, false, true, false, false, false }; pdb.set_mirror( pmirror );

After this we are ready for the Vlasov iteration. The first iteration is a special case because of the initial guess for the plasma and therefore after it, the problem had to be redefined to use the non-linear plasma model.

for( size_t i = 0; i < 15; i++ ) { if( i == 1 ) { double rhoe = pdb.get_rhosum(); solver.set_pexp_plasma( -rhoe, Te, Up ); } solver.solve( epot, scharge ); efield.recalculate(); pdb.clear(); pdb.add_2d_beam_with_energy( 15000, 600.0, 1.0, 4.0, 5.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0015 ); pdb.iterate_trajectories( scharge, efield, bfield ); }

The particle database calculates the amount of space charge at the
definition point and that can be retrieved with
*get_rhosum()*. This value is then fed to the plasma model to
define the space charge of background electrons in the plasma. The
function *set_pexp_plasma()* enables the positive exponential
plasma model and defines plasma density (first argument), electron
temperature (second argument) and plasma potential (third
argument). The electric potential is solved using the solver with
*solver.solve( epot, scharge )*. After this the electric field
needs to be recalculated with *efield.recalculate()*. The rest of
the Vlasov iteration loops continue as in the first example. Also in
the plasma extraction problems it might be necessary to use
under-relaxation of space charge field in hard cases.

In this example we invoce the interactive plotter to show the results. You need to provide all the data you want to use in the plotter. Please note that you need to provide the argc and argv parameters to the plotter as the underlying GTK library needs these in initialization.

#ifdef GTK3 GTKPlotter plotter( argc, argv ); plotter.set_geometry( &geom ); plotter.set_epot( &epot ); plotter.set_scharge( &scharge ); plotter.set_particledatabase( &pdb ); plotter.new_geometry_plot_window(); plotter.run(); #endif

That is it again. The complete source is available here: plasmacyl.cpp and the Makefile: Makefile.

The default screen for the plotter looks something like this:

With the interactive tool you can make plots from different view points and you can modify the plotted features using menus. You can also make hardcopies in png, svg, eps and pdf formats. You can plot geometries, emittances, profiles (scatter or colormap style), field values, etc. As an example we have here an emittance plot from this simulation with fitted emittance ellipse: