You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: README.rst
+1-1
Original file line number
Diff line number
Diff line change
@@ -7,7 +7,7 @@ PUNC++
7
7
8
8
*Particles-in-UNstructured-Cells, C++ version* (PUNC++) is an open source scientific program for simulating plasmas using the *Particle-In-Cell* (PIC) method on an unstructured mesh. The field quantities are solved using the *Finite Element Method* (FEM) through the third party finite element environment FEniCS_. The focus is on plasma-object interaction, typically for space applications.
9
9
10
-
A user guide is available online on ReadTheDocs_, whereas for developers, the code is documented using Doxygen. See the installation section for building Doxygen code documentation.
10
+
A user guide is available online on ReadTheDocs_, whereas for developers, the code is documented using Doxygen. See the installation section for building developer's documentation.
While FEniCS is a very powerful problem solving environment for finite elements, it has proven notoriously troublesome to install, and we do not have the capacity to offer much guidance here. Since FEniCS and its dependencies change functionality frequently, we recommend using the version listed on this page. Several installation methods are described on the FEniCS website. We point out that depending on how FEniCS and its dependencies are installed, there may be a tenfold difference in the performance of PUNC++, since not all installations have the same level of optimization. The Anaconda_ version yields good performance, but make sure the folder containing ``libdolfin.so`` is accesible to the compiler, e.g. by exporting it to ``LD_LIBRARY_PATH``.
26
+
Since updates to FEniCS (and its dependencies) frequently breaks backwards compatibility, we recommend using the version listed on this page. Several installation methods are described on the FEniCS website. We point out that depending on how FEniCS and its dependencies are installed, there may be a tenfold difference in the performance of PUNC++, since not all installations have the same level of optimization. The Anaconda_ version yields good performance, but make sure the folder containing ``libdolfin.so`` is accesible to the compiler, e.g., by exporting it to ``LD_LIBRARY_PATH``.
Once you have installed all dependencies, download PUNC++::
30
+
For developers it may be of interest to make themselves familiar with FEniCS. The `FEniCS tutorial`_ is of relevance, especially the section called "Fundamentals: Solving the Poisson equation".
@@ -38,23 +42,23 @@ PUNC++ comes in two parts: a library (in the folder ``punc++/punc``) and an exec
38
42
cd punc
39
43
./build.sh # or ./install.sh
40
44
41
-
The library will be installed in a local folder ``punc++/punc/install``. It can also be installed in system-wide folders by running ``./install.sh`` instead of ``./build.sh``. For developers, the documentation of the code can be installed by running::
45
+
The library will be installed in a local folder ``punc++/punc/install``. It can also be installed in system-wide folders by running ``./install.sh`` instead of ``./build.sh``. Developer's documentation of the library can be generated by running::
42
46
43
47
./doc.sh
44
48
45
-
and launching ``doc.html``. Next, to install ``interaction``::
49
+
The documentation is available in the file ``doc.html`` just generated. To install ``interaction``::
46
50
47
51
cd ../interaction
48
52
./build.sh
49
53
50
-
That's it! In addition we recommend installing the following tools for pre- and post-processing:
54
+
In addition we recommend installing the following tools for pre- and post-processing:
51
55
52
-
- Gmsh_
56
+
- Gmsh_ 3
53
57
- ParaView_
54
58
- Metaplot_
55
59
56
60
.. _Gmsh: http://gmsh.info/
57
61
.. _ParaView: https://www.paraview.org/
58
62
.. _Metaplot: https://metaplot.readthedocs.io
59
63
60
-
The tutorials in this user guide will make use of them.
64
+
The tutorials in this user guide will make use of them. Note that Gmsh changed the format of its files in version 4, and FEniCS 2018.1.0 is not compatible with the new version.
In this tutorial, we will carry out a 2D simulation of the current collected by a cylindrical Langmuir probe in non-drifting, non-magnetized, Maxwellian proton-electron plasma. The physical parameters of the problem is given below:
We start by making a new directory for the project, and make a symbolic link of the ``interaction`` executable in this folder (modify paths as necessary)::
18
+
We start by making a new directory for the project, and make a symbolic link to the ``interaction`` executable in this folder (modify paths as necessary)::
19
19
20
20
cd ~
21
21
mkdir -p Projects/tutorials/basic
22
22
cd Projects/tutorials/basic
23
23
ln -s ~/punc++/interaction/build/interaction
24
24
25
+
If using a system-wide installation of PUNC++, you can omit creating a symbolic link.
26
+
25
27
Mesh generation
26
28
---------------
27
-
The first thing we need is a mesh of the simulation domain including any objects.
28
-
For a 2D simulation of a cylindrical probe this is just two concentric circles.
29
-
The inner circle, or interior boundary, is the probe itself.
30
-
The outer circle is the exterior boundary of the simulation domain.
29
+
Let's start by creating a mesh.
30
+
For a 2D simulation of a cylindrical probe the probe is represented by a circular, interior boundary.
31
+
We shall also use a circlular exterior boundary.
31
32
It is an assumption of the underlying numerical methods that the exterior boundary is in the background plasma and is not affected by any local perturbations in the electric potential.
32
33
This makes it necessary for the outer boundary to be outside the sheath of the probe.
33
34
The closer the outer boundary is to the probe, the less valid the assumption becomes, and the more it limits the accuracy of the simulation.
34
35
On the other hand, increasing the domain size increases the cost of the simulation.
35
-
The mesh gets bigger, and you need more simulation particles to maintain the same density of particles in the domain.
36
-
You may in fact also need to run the simulation for a longer physical time in order to reach a sufficiently steady state.
37
-
This is because the domain is uniformly filled with particles before the simulation starts, and you get rather violent, unphysical transient behaviour where the plasma is perturbed, e.g. where the sheath should be.
38
-
These transients may cause waves which takes a long time to leave the domain.
39
-
Moreover, in improving the accuracy of a simulation, the clue is to improve the limiting factor, which may or may not be the radius of the outer boundary.
40
-
Thus, finding the right set of simulation parameters requires some experimentation.
36
+
The mesh gets bigger, and more simulation particles are needed to maintain the same density of particles in the domain.
37
+
It may in fact also be necessary to run the simulation for a longer physical time in order to reach a sufficiently steady state.
38
+
This is because the domain is uniformly filled with particles before the simulation starts, and rather violent, unphysical transients occur where the plasma is perturbed, e.g. where the sheath should be.
39
+
These transients may cause waves which takes a long time to settle.
40
+
In improving the accuracy of a simulation, one should attempt to improve the limiting factor, which may or may not be the radius of the outer boundary.
41
+
Finding the right set of simulation parameters requires some experimentation.
41
42
We will use an outer boundary of 100mm radius.
42
43
43
44
This geometry can be created either using the Gmsh *Graphical User Interface* (GUI), or it can be specified directly in a text file for Gmsh to read.
44
-
We shall call our geometry file ``cylinder.geo`` and it looks as follows:
On the first four lines we define variables for the inner and outer radii, as well as the resolution of the mesh on the inner and outer boundaries.
49
-
The resolution must be sufficiently fine to resolve the characteristic lengths of the physical processes involved.
50
-
This means that it must be fine enough to resolve the probe on the inner boundary, and not much bigger than the electron Debye length on the outer boundary.
51
-
Insufficiently resolving the electron Debye length is known to cause numerical heating of the plasma on rectangular meshes (see Birdsall&Langdon), and it is reasonable to expect simlar effects for unstructured meshes.
52
-
Since we are most interested in what happens close to the probe, we can often get away with a resolution that is actually somewhat bigger than the Debye length at the outer boundary.
49
+
On the first four lines variables are defined for the inner and outer radii, as well as the resolution of the mesh on the inner and outer boundaries.
50
+
The resolution will vary continually between the boundaries.
51
+
The resolution must be sufficiently fine to resolve any characteristic lengths of the physics involved.
52
+
The outer boundary must therefore have a resolution not much bigger than the electron Debye length, whereas the inner boundary must in addition be sufficiently fine to resolve the circular cross-section of the probe.
53
+
Insufficiently resolving the electron Debye length is known to cause numerical heating of the plasma on rectangular meshes, and it is reasonable to expect simlar effects for unstructured meshes.
54
+
Since we are mostly interested in what happens close to the probe, we can often get away with a resolution that is actually somewhat bigger than the Debye length at the outer boundary.
53
55
54
-
Next, we add a point in the center of the domain, as well as north, east, south and west of the center on each boundary, and connect them by circle arcs. (It may be convenient to define the variables and add the points in text, and use the GUI to do the rest).
56
+
The lines starting with ``Point`` defines a point in the center of the domain, as well as to the north, east, south and west of the center, on both boundaries.
57
+
Each boundary is made up of four circle arcs (``Circle``) connected between these points.
55
58
56
-
It is important that the proper *Physical Groups* are created in Gmsh. PUNC++ needs one physical group for each boundary (*Physical Lines* in 2D or *Physical Surfaces* in 3D).
59
+
It is important that the proper *Physical Groups* are created in Gmsh. PUNC++ needs one physical group for each boundary (``Physical Lines`` in 2D or ``Physical Surfaces`` in 3D).
57
60
Note that it is *mandatory* that the exterior boundary has the lowest id number of all boundaries since this is how PUNC++ knows which boundary is the exterior one.
58
-
It is also possible to define a physical group for the domain (*Physical Surface* in 2D or *Physical Volume* in 3D), but this is optional.
61
+
It is also possible to define a physical group for the domain (``Physical Surface`` in 2D or ``Physical Volume`` in 3D), but this is optional.
62
+
63
+
Hint: For larger geometries it can be difficult to keep track of all line in a text editor. The author prefers to define the variables and the points in text, and subsequently opening it in Gmsh to connect arcs and lines between the points.
59
64
60
65
To generate the mesh from the geometry, run ``Mesh -> 2D`` from the GUI or from the terminal::
61
66
62
67
$ gmsh -2 cylinder.geo
63
68
64
-
The mesh is named ``cylinder.msh`` and should look something like this when opened in Gmsh:
69
+
For a 3D mesh it would be ``-3`` instead of ``-2``. The mesh is named ``cylinder.msh`` and should look something like this when opened in Gmsh:
65
70
66
71
.. image:: mesh.png
67
72
73
+
The finite element approximation of the fields will be better the less sliver the cells in the mesh are. A gamma factor of 1 indicates a completely equilateral/regular cell, whereas a gamma of 0 indicates a degenerate cell. The quality of the mesh can be inspected from Gmsh by clicking ``Tools -> Statistics`` and then ``Update``. The average, minimum and maximum gamma factor should now be displayed, and it is possible to plot the distribution of the gamma values. It may look something like this:
74
+
75
+
.. image:: mesh_statistics.png
76
+
77
+
Make sure most cells are above 0.3 and that the minimum value is not zero, i.e., that there are no degenerate cells.
78
+
68
79
Mesh files must be in either DOLFIN XML or HDF5 formats.
69
80
To convert the Gmsh mesh, use FEniCS's own conversion tool::
70
81
@@ -74,6 +85,76 @@ You should now have the files ``cylinder.xml``, ``cylinder_facet_region.xml`` an
74
85
75
86
Running the simulation
76
87
----------------------
88
+
The ``interaction`` executable loads simulation settings from an ``ini`` file which we shall call ``setup.ini`` and put in the same folder as the mesh:
Most settings in this file are pretty self-explanatory.
93
+
The mesh is found in ``cylinder.xml`` and the mesh is run for :math:`1\,\mathrm{\mu s}`.
94
+
Note that ``cylinder_facet_region.xml`` is also used although not explicitly mentioned.
95
+
The simulation consists of electrons and protons.
96
+
In particular, it contains :math:`5\cdot10^5` macroparticles of each.
97
+
Notice also that some settings take suffixes (units).
98
+
Mass can for instance be specified in both electron masses (``me``) or atomic mass units (``amu``).
99
+
100
+
The exterior boundary should, as mentioned, have the lowest id in the ``.geo`` file.
101
+
The second lowest boundary id is taken to be object 0, the third lowest id object 1, and so forth.
102
+
The objects are always perfect electric conductors, and by default their floating potential will be self-consistently determined from collected charges.
103
+
In our case, however, we would like to fix the potential of the cylinder to 3V.
104
+
We do this by adding a voltage source (``vsource``) between system ground (``-1``) and the object (``0``) of 3V.
105
+
Multiple voltage and current sources can be defined between two objects, or between objects and sytem ground by adding several ``vsource`` and ``isource`` entries under the ``[objects]`` section.
106
+
107
+
For this simulation we have also included an optional diagnostics, namely that the electric potential be saved every 3 nanoseconds.
108
+
109
+
Finally, to run the simulation, type::
110
+
111
+
$ ./interaction setup.ini
112
+
113
+
Omit ``./`` if using system-wide installation. After the simulation is complete, three ``dat`` files should appear. ``history.dat`` contains time-series diagnostics. ``population.dat`` contains all the particles at the end of the simulation and ``state.dat`` contains auxiliary variables necessary to continue the simulation from where it stopped. After having a look at the data, we would probably conclude that :math:`1\,\mathrm{\mu s}` was a bit short, and decide to increase it to :math:`2\,\mathrm{\mu s}`. We chould then simply change ``setup.ini`` and execute the program again to pick up the simulation where we left off. It is also possible to stop a simulation gracefully by pressing ``Ctrl+C`` once, and continuing the simulation later. Pressing ``Ctrl+C`` a second time kills the program instantly, without possibility for continuation. To restart a simulation from scratch, remove one of the ``dat`` files.
114
+
115
+
It is also possible to override a setting in an ``ini`` file by command line arguments. For instance, to expand the simulation time, we could also execute::
116
+
117
+
$ ./interaction setup.ini --time.stop "2e-6 s"
118
+
119
+
``interaction`` supports a wide range of settings, each of which can be set either in the ``ini`` file or as a command line argument. The full list is available from::
120
+
121
+
$ ./interaction --help
77
122
78
123
Inspecting the results
79
124
----------------------
125
+
126
+
Time-series are stored in a tabulated ASCII file ``history.dat``. The first few lines for our simulation is demonstrated below:
127
+
128
+
.. literalinclude:: history.dat
129
+
130
+
The first six columns shows the time-step, the physical time at the time-step, the number of negatively charged particles (electrons), the number of positively charged particles (ions), total kinetic energy of the particles, and total potential energy of the particles, respectively. Subsequently follows three columns for each objects, its voltage, collected current, and charge.
131
+
132
+
The commented lines in the header follows the syntax of Metaplot, which allows for easy visualization. To plot the current collected by the probe (which is object zero)::
133
+
134
+
$ mpl history.dat "I[0]"
135
+
136
+
.. image:: current.png
137
+
138
+
The current collected should according to OML theory be :math:`-35.3\,\mathrm{\mu A}`. Since the simulation is quite rough it is a bit hard to see how close it actually is. Let us perform an *Exponential Moving Average* (EMA) with a relaxation time of :math:`0.1\,\mathrm{\mu s}` and compare to the OML theory::
This is not too bad given that the simulation is quite rough. Let us also plot the number of electrons and ions::
145
+
146
+
$ mpl history.dat ne ni
147
+
148
+
.. image:: num_of_particles.png
149
+
150
+
The number of electrons seems to have reached steady-state. The number of ions have not had time to change very much. We not that Metaplot is a very experimental program at this stage. We do not delve further into it.
151
+
152
+
Finally, let us have a look at the electric potential we stored every :math:`3\,\mathrm{ns}`. The field quantities are stored in the ``fields`` folder in VTK format, and is easily visualized using ParaView. To open the electric potential::
153
+
154
+
$ paraview fields/phi.pvd
155
+
156
+
Inside ParaView, click ``Apply`` and then the green play symbol to start an animation. It should look something like the following:
157
+
158
+
.. image:: phi.gif
159
+
160
+
Initially the plasma is uniformly distributed and the potential is that of a cylinder in vacuum. As the sheath starts forming due to the surrounding plasma, acoustic shock waves emerges and travel outwards, eventually forming standing waves between the boundaries before they die out and the correct solution is arrived at at steady-state.
0 commit comments