# Vlasov2d

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

This simulation is a modification of the vlasov2d test
(*tests/vlasov2d.cpp*). It finds out how a 3 keV, 5 mA/cm2 H+
beam behaves in a system of three electrodes. For simulating this
problem a program needs to be written in C++ to control the IBSimu
library that you have installed. Therefore you need to write a C++
file, named *vlasov2d.cpp* for example.

First of all your program needs to access the the IBSimu library. For this you need to include some header files in your cpp-file:

#include "epot_bicgstabsolver.hpp" #include "geometry.hpp" #include "func_solid.hpp" #include "epot_efield.hpp" #include "meshvectorfield.hpp" #include "particledatabase.hpp" #include "geomplotter.hpp" #include "ibsimu.hpp" #include "error.hpp"

The needed include files depends on the functions and classes that are used by your program. More information can be found in the reference manual. Then you program need to have a main function as an entrance point to the program:

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

In here we are using the try-catch environment of C++ that will
catch errors within the IBSimu library and present some sensible error
message. The error message will include the source file, line number
and the name of the function where the error happened in addition to
the message itself. The *ibsimu* is a global object that can be
used to control some basic parameters of the code. The function
*set_message_threshold()* is used to set the amount of output
printed to stdout (or output file) by the library functions. Zero
means no output, one is standard output, two is extended output. The
target for library messages can be set using function
*set_message_output()* of the ibsimu object. Some parts of
IBSimu can use multiple CPUs (or cores) by splitting the calculation
into several threads. The number of threads created can be set using
*set_thread_count()*.

Next we write the *simu()* function, the actual workhorse of
the simulation. First we need to define the geometry where the
simulation is done:

Geometry geom( MODE_2D, Int3D(241,101,1), Vec3D(0,0,0), 0.0005 );

The first argument specifies that the simulation is planar two
dimensional. Other possibilities would be *MODE_1D* for one
dimensional simulations, *MODE_CYL* for cylindrically
symmetrical simulations and *MODE_3D* full three
dimensions. The second argument defines the size of the simulation
mesh. In this case it is 241 nodes along x-axis and 101 nodes along
y-axis and 1 node along z-axis. The third argument specifies the
location of the node *(0,0,0)* in meters. The fourth argument
defines the size of mesh, meaning the length from one node to the next
in meters. The size of the simulation area is therefore 240x0.5 mm by
100x0.5 mm.

Next we need to define the electrodes in the geometry. There are
several ways of defining them. In this example we will use the most
fundamental way, the *FuncSolid* objects. Other ways include
Constructive Solid Geometry definition, two dimensional DXF-file based
definition and three dimensional STL-files. With *FuncSolid*
you need to write a function, which takes in the coordinates of a
point and returns *true* if the point is inside the solid and
*false* otherwise. In this example we have three
electrodes. First we need to define the functions:

bool solid1( double x, double y, double z ) { return( x <= 0.02 && y >= 0.018 ); } bool solid2( double x, double y, double z ) { return( x >= 0.03 && x <= 0.04 && y >= 0.02 ); } bool solid3( double x, double y, double z ) { return( x >= 0.06 && y >= 0.03 && y >= 0.07 - 0.5*x ); }

The *FuncSolid* objects have to be defined and inserted in the geometry:

Solid *s1 = new FuncSolid( solid1 ); geom.set_solid( 7, s1 ); Solid *s2 = new FuncSolid( solid2 ); geom.set_solid( 8, s2 ); Solid *s3 = new FuncSolid( solid3 ); geom.set_solid( 9, s3 );

First the *FuncSolid* object is made using the *new*
operator of C++. This object is then inserted in the geometry using
*set_solid* function. The set_solid function is given the
number of the electrode, which is used to reference the
electrode. These numbers start from 7 for the user defined electrodes,
because numbers from 1 to 6 are used to reference the boundaries of
the simulation box.

Now that the geometry is defined, the potentials of the electrodes and the boundary conditions for the simulation box have to be defined. The electrode numbers will be used to reference the different boundaries. Electrode 1 is the minimum x boundary, 2 is maximum x boundary, 3 is minimum y, 4 is maximum y, 5 is minimum z and 6 is maximum z. Starting from 7 you reference to the user defined electrodes. The simulation box boundaries can have either Neumann or Dirichlet boundary condition. Neumann means constant first derivative of the potential with respect to the unit outward normal of the surface. Dirichlet means constant potential on the surface. User defined electrodes can only have Dirichlet boundary condition. In our example we define as follows:

geom.set_boundary( 1, Bound(BOUND_DIRICHLET, -3.0e3) ); geom.set_boundary( 2, Bound(BOUND_DIRICHLET, -1.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, -3.0e3) ); geom.set_boundary( 8, Bound(BOUND_DIRICHLET, -14.0e3) ); geom.set_boundary( 9, Bound(BOUND_DIRICHLET, -1.0e3) );

In the example we set a fixed -3 kV potential on the xmin boundary, -1 kV potential on the xmax boundary and standard Neumann boundary condition (zero first derivative of potential) on ymin and ymax boundaries. The user defined electrodes get -3 kV, -14 kV and -1 kV potentials. After these definitions we can build the mesh:

geom.build_mesh();

This function maps the simulation area and finds where the
electrode boundaries are. Then we are ready to formulate the problem
of solving the Poisson equation. The problem is contained in the
*EpotSolver* object, in this case we use the EpotBiCGSTABSolver:

EpotBiCGSTABSolver solver( geom );

In this example we use the default stablized bi-conjugate gradient iterative solver. It is included with IBSimu and doesn't depend on any external libraries. For 2D problems the most efficient solver currently is the UMFPack solver but it is dependen't on the external SuiteSparse library and it's existence depends on the parameters given for the compiler during the compilation of the library.

After this we initialize some other object that will be needed in the calculation:

EpotField epot( geom ); MeshScalarField scharge( geom ); MeshVectorField bfield; EpotEfield efield( epot ); field_extrpl_e efldextrpl[6] = { FIELD_EXTRAPOLATE, FIELD_EXTRAPOLATE, FIELD_SYMMETRIC_POTENTIAL, FIELD_EXTRAPOLATE, FIELD_EXTRAPOLATE, FIELD_EXTRAPOLATE }; efield.set_extrapolation( efldextrpl );

The object names are rather self-explanatory: *MeshScalarField*
is a mesh-based scalar field in the geometry of *geom*. We use it for
storing electric potential and space charge density fields. The
particle calculation also need a magnetic field (even if it is
zero). The mesh of the magnetic field is not bound to the mesh of the
electric field and therefore it doesn't use a lot of memory if it is a
constant field. The magnetic field is of course a
*MeshVectorField*. The electric field *EpotEfield* is a
special object, which doesn't contain any data on it's own. It is used
to calculate the electric field at a location from the electric
potential data. The electric field is told what is returned for the
electric field outside the simulation geometry. This is needed for the
particle calculations close by to the boundaries. In this example the
geometry is symmetrical so we define the ymin boundary to mirror the
efield.

Before getting into the vlasov iteration loop we still need to
define the particle calculation system, the
*ParticleDataBase*:

ParticleDataBase2D pdb( geom ); bool pmirror[6] = { false, false, true, false, false, false }; pdb.set_mirror( pmirror );

The particle calculation object needs to be selected according to
the geometry of the simulation. In this case we use the
*ParticleDataBase2D* because the geometry is two
dimensional. By default, the particles, which hit the
simulation boundaries are removed from the simulation. In case of
symmetrical simulations it might be useful to have particles mirrored
at boundaries. This can be set using *set_mirror()*. The last
parameter set in the particle database is the type of interpolation
used. The iterator calculates the particle trajectory in discrete
steps.

Then we are finally ready to start the Vlasov iteration to find a
self-consistent solution for the electric potential and beam
propagation in the geometry. This is achieved by first solving the
electric potential in the geometry without the beam (zero space
charge, which is the default condition for the MeshScalarField
scharge). The electric field is then calculated from the electric
potential with *efield.recalculate()*. Then the beam particles
are iterated through the geometry and space charge from the beam is
calculated. This space charge is then used to calculate the next
iteration of electric potential field. This is repeated until
convergence. In the example we do 5 such loops:

for( size_t i = 0; i < 5; i++ ) { solver.solve( epot, scharge ); efield.recalculate(); pdb.clear(); pdb.add_2d_beam_with_energy( 1000, 50.0, 1.0, 1.0, 3.0e3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.012 ); pdb.iterate_trajectories( scharge, efield, bfield ); }

The particle beam in this example is defined using the function
*add_2d_beam_with_energy()*. The first argument for the
function specifies the number of particle trajectories
defined. Typically the number of particle trajectories has to be
enough to represent the total beam well enough within each mesh
square. In this case we have no temperature distribution (zero
emittance) in the beam so quite a small number of trajectories is
enough. In this case we have 1000 trajectories, which is about 40
trajectories per mesh unit. The rest of the function arguments specify
the current density of the beam (A/m2), particle charge (in electron
charges), particle mass (in atomic units), starting energy (eV),
parallel tempererature (eV), transverse temperature (eV), x1, y1, x2
and y2 coordinates of the starting line, where the beam starts from
(meters).

After sufficient amount of iteration loops the solution has hopefully converged. Sometimes it is necessary to use under-relaxation techniques for the space charge field to have convergence. Unstabilities are most likely encountered when beam intensity is close to being space charge limited.

After the loops we plot the results. Of course you could (or should) also make plots during the iteration to verify that the solution converges. In this case we only plot a simple geometry plot showing the geometry, the beam and the electric field using equipotential lines:

GeomPlotter geomplotter( geom ); geomplotter.set_size( 750, 750 ); geomplotter.set_epot( &epot ); geomplotter.set_particle_database( &pdb ); geomplotter.plot_png( "plot1.png" );

OK. That is it. The complete source file is available here: vlasov2d.cpp. You also a need a Makefile to build the simulation executable. The resulting plot file showing the system follows: