diff --git a/CITATION.cff b/CITATION.cff deleted file mode 100644 index 18cd6d7..0000000 --- a/CITATION.cff +++ /dev/null @@ -1,50 +0,0 @@ -# This CITATION.cff file was generated with cffinit. -# Visit https://bit.ly/cffinit to generate yours today! - -cff-version: 1.2.0 -title: Finite element Particle Kinetic Code -message: >- - If you use this software, please cite it using the - metadata from this file. -type: software -authors: - - given-names: Jorge - family-names: Gonzalez - email: jorge.gonzalez@upm.es - affiliation: Universidad Politécnica de Madrid - orcid: 'https://orcid.org/0000-0001-7905-5001' -repository-code: 'https://gitlab.com/JorgeGonz/fpakc' -abstract: >- - Welcome to fpakc (Finite element PArticle Kinetic Code), a - modern object oriented Fortran open-source code for - particle simulations of plasma and gases. This code works - by simulating charged and neutral particles, following - their trajectories, collisions and boundary conditions - imposed by the user. - - One of our aims is to make a code easy to maintain as well - as easy to use by a variety of reserchers and students. - - This code is currenlty in very early steps of development. - - The code aims to be easy to maintain and easy to use, - allowing its application from complex problems to easy - examples that can be used, for example, as teaching - exercises. - - Parallelization techniques such as OpenMP, MPI will be - used to distribute the cpu load. We aim to make fpakc GPU - compatible in the future. - - The codefpakc makes use of finite elements to generate - meshes in complex geometries. Particle properties are - deposited in the nodes and cells of the mesh. The - electromagnetic field, with the boundary conditions - imposed by the user, is solved also in this mesh. -keywords: - - particle-in-cell - - plasma - - finite elements -license: GPL-3.0 -version: beta -date-released: '2025-10-01' diff --git a/data/collisions/IO_e-Kr.dat b/data/collisions/IO_e-Kr.dat deleted file mode 100644 index 532ca4d..0000000 --- a/data/collisions/IO_e-Kr.dat +++ /dev/null @@ -1,52 +0,0 @@ -# D108525 "refs": {"B56": {"note": "CLM-R294 (1989)"}} -# Relative energy (eV) cross section (m^2) -1.40E+01 0 -1.62E+01 7.249E-21 -1.88E+01 1.199E-20 -2.18E+01 1.644E-20 -2.53E+01 2.1E-20 -2.94E+01 2.542E-20 -3.41E+01 2.937E-20 -3.95E+01 3.26E-20 -4.58E+01 3.499E-20 -5.32E+01 3.653E-20 -6.17E+01 3.726E-20 -7.15E+01 3.728E-20 -8.29E+01 3.671E-20 -9.62E+01 3.566E-20 -1.12E+02 3.426E-20 -1.29E+02 3.259E-20 -1.50E+02 3.075E-20 -1.74E+02 2.881E-20 -2.02E+02 2.682E-20 -2.34E+02 2.484E-20 -2.72E+02 2.289E-20 -3.15E+02 2.101E-20 -3.65E+02 1.922E-20 -4.24E+02 1.751E-20 -4.91E+02 1.592E-20 -5.70E+02 1.443E-20 -6.61E+02 1.305E-20 -7.67E+02 1.177E-20 -8.89E+02 1.06E-20 -1.03E+03 9.526E-21 -1.20E+03 8.547E-21 -1.39E+03 7.658E-21 -1.61E+03 6.851E-21 -1.87E+03 6.121E-21 -2.16E+03 5.462E-21 -2.51E+03 4.868E-21 -2.91E+03 4.334E-21 -3.38E+03 3.855E-21 -3.92E+03 3.426E-21 -4.54E+03 3.041E-21 -5.27E+03 2.698E-21 -6.11E+03 2.391E-21 -7.09E+03 2.118E-21 -8.22E+03 1.875E-21 -9.53E+03 1.658E-21 -1.11E+04 1.466E-21 -1.28E+04 1.295E-21 -1.49E+04 1.143E-21 -1.72E+04 1.009E-21 -2.00E+04 8.898E-22 diff --git a/data/collisions/IO_e-Xe.dat b/data/collisions/IO_e-Xe.dat deleted file mode 100644 index 3067578..0000000 --- a/data/collisions/IO_e-Xe.dat +++ /dev/null @@ -1,52 +0,0 @@ -# EL cross sections extracted from PROGRAM MAGBOLTZ, VERSION 7.1 JUNE 2004 www.lxcat.net/Biagi-v7.1 -# Relative energy (eV) cross section (m^2) -1.21E+01 0 -1.41E+01 3.923E-21 -1.64E+01 1.194E-20 -1.91E+01 2.1E-20 -2.22E+01 2.946E-20 -2.58E+01 3.65E-20 -3.00E+01 4.185E-20 -3.49E+01 4.552E-20 -4.06E+01 4.766E-20 -4.72E+01 4.85E-20 -5.49E+01 4.828E-20 -6.39E+01 5.031E-20 -7.43E+01 5.1E-20 -8.64E+01 5.1E-20 -1.01E+02 5.032E-20 -1.17E+02 4.906E-20 -1.36E+02 4.732E-20 -1.58E+02 4.521E-20 -1.84E+02 4.283E-20 -2.14E+02 4.029E-20 -2.49E+02 3.764E-20 -2.90E+02 3.497E-20 -3.37E+02 3.233E-20 -3.92E+02 2.975E-20 -4.56E+02 2.726E-20 -5.31E+02 2.489E-20 -6.17E+02 2.266E-20 -7.18E+02 2.056E-20 -8.35E+02 1.861E-20 -9.72E+02 1.68E-20 -1.13E+03 1.514E-20 -1.32E+03 1.361E-20 -1.53E+03 1.221E-20 -1.78E+03 1.094E-20 -2.07E+03 9.781E-21 -2.41E+03 8.735E-21 -2.80E+03 7.789E-21 -3.26E+03 6.938E-21 -3.79E+03 6.171E-21 -4.41E+03 5.484E-21 -5.13E+03 4.868E-21 -5.97E+03 4.316E-21 -6.94E+03 3.824E-21 -8.07E+03 3.385E-21 -9.39E+03 2.994E-21 -1.09E+04 2.646E-21 -1.27E+04 2.336E-21 -1.48E+04 2.062E-21 -1.72E+04 1.818E-21 -2.00E+04 1.602E-21 diff --git a/data/see/constant_yield.dat b/data/see/constant_yield.dat new file mode 100644 index 0000000..060c8ad --- /dev/null +++ b/data/see/constant_yield.dat @@ -0,0 +1,4 @@ +#Relative energy (eV) yield () +0.000 1.000E-01 +1.000 1.000E-01 + diff --git a/doc/logos/fpakc.pdf b/doc/logos/fpakc.pdf deleted file mode 100644 index 36bc5ab..0000000 Binary files a/doc/logos/fpakc.pdf and /dev/null differ diff --git a/doc/user-manual/.gitignore b/doc/user-manual/.gitignore index a8475d7..61ad7e0 100644 --- a/doc/user-manual/.gitignore +++ b/doc/user-manual/.gitignore @@ -6,7 +6,6 @@ *.aux *.ps bibliography.bib.bak -bibliography.bib.sav *.bbl *.blg *.out @@ -20,3 +19,4 @@ fpakc_UserManual-blx.bib *.gls *.ist fpakc_UserManual.run.xml +*.bib.sav diff --git a/doc/user-manual/figures/logos/PPartiC.eps b/doc/user-manual/figures/logos/PPartiC.eps new file mode 100644 index 0000000..5517180 --- /dev/null +++ b/doc/user-manual/figures/logos/PPartiC.eps @@ -0,0 +1,177 @@ +%!PS-Adobe-3.0 EPSF-3.0 +%%Creator: cairo 1.16.0 (https://cairographics.org) +%%CreationDate: Wed Oct 07 12:01:13 2020 +%%Pages: 1 +%%DocumentData: Clean7Bit +%%LanguageLevel: 2 +%%BoundingBox: 6 7 205 57 +%%EndComments +%%BeginProlog +50 dict begin +/q { gsave } bind def +/Q { grestore } bind def +/cm { 6 array astore concat } bind def +/w { setlinewidth } bind def +/J { setlinecap } bind def +/j { setlinejoin } bind def +/M { setmiterlimit } bind def +/d { setdash } bind def +/m { moveto } bind def +/l { lineto } bind def +/c { curveto } bind def +/h { closepath } bind def +/re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto + 0 exch rlineto 0 rlineto closepath } bind def +/S { stroke } bind def +/f { fill } bind def +/f* { eofill } bind def +/n { newpath } bind def +/W { clip } bind def +/W* { eoclip } bind def +/BT { } bind def +/ET { } bind def +/BDC { mark 3 1 roll /BDC pdfmark } bind def +/EMC { mark /EMC pdfmark } bind def +/cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def +/Tj { show currentpoint cairo_store_point } bind def +/TJ { + { + dup + type /stringtype eq + { show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse + } forall + currentpoint cairo_store_point +} bind def +/cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore + cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def +/Tf { pop /cairo_font exch def /cairo_font_matrix where + { pop cairo_selectfont } if } bind def +/Td { matrix translate cairo_font_matrix matrix concatmatrix dup + /cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point + /cairo_font where { pop cairo_selectfont } if } bind def +/Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def + cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def +/g { setgray } bind def +/rg { setrgbcolor } bind def +/d1 { setcachedevice } bind def +/cairo_data_source { + CairoDataIndex CairoData length lt + { CairoData CairoDataIndex get /CairoDataIndex CairoDataIndex 1 add def } + { () } ifelse +} def +/cairo_flush_ascii85_file { cairo_ascii85_file status { cairo_ascii85_file flushfile } if } def +/cairo_image { image cairo_flush_ascii85_file } def +/cairo_imagemask { imagemask cairo_flush_ascii85_file } def +%%EndProlog +%%BeginSetup +%%EndSetup +%%Page: 1 1 +%%BeginPageSetup +%%PageBoundingBox: 6 7 205 57 +%%EndPageSetup +q 6 7 199 50 rectclip +1 0 0 -1 0 71 cm q +0 g +49.766 27.305 m 49.766 29.898 49.078 31.813 47.703 33.039 c 46.328 34.27 + 44.297 34.883 41.609 34.883 c 36.625 34.883 l 36.625 47.945 l 32.438 47.945 + l 32.438 14.352 l 41.516 14.352 l 44.234 14.352 46.285 14.953 47.672 16.148 + c 49.066 17.348 49.766 19.242 49.766 21.836 c h +45.594 21.398 m 45.594 20.18 45.313 19.301 44.75 18.758 c 44.195 18.219 + 43.313 17.945 42.094 17.945 c 36.625 17.945 l 36.625 31.289 l 42.094 31.289 + l 43.313 31.289 44.195 31.008 44.75 30.445 c 45.313 29.883 45.594 29 45.594 + 27.789 c h +45.594 21.398 m f +75.641 27.305 m 75.641 29.898 74.953 31.813 73.578 33.039 c 72.203 34.27 + 70.172 34.883 67.484 34.883 c 62.5 34.883 l 62.5 47.945 l 58.313 47.945 + l 58.313 14.352 l 67.391 14.352 l 70.109 14.352 72.16 14.953 73.547 16.148 + c 74.941 17.348 75.641 19.242 75.641 21.836 c h +71.469 21.398 m 71.469 20.18 71.188 19.301 70.625 18.758 c 70.07 18.219 + 69.188 17.945 67.969 17.945 c 62.5 17.945 l 62.5 31.289 l 67.969 31.289 + l 69.188 31.289 70.07 31.008 70.625 30.445 c 71.188 29.883 71.469 29 71.469 + 27.789 c h +71.469 21.398 m f +99.891 47.945 m 99.473 47.945 99.016 47.82 98.516 47.57 c 98.023 47.313 + 97.602 46.957 97.25 46.508 c 96.508 47.469 95.535 47.945 94.328 47.945 +c 90.047 47.945 l 87.617 47.945 85.852 47.461 84.75 46.492 c 83.645 45.516 + 83.094 43.836 83.094 41.461 c 83.094 40.508 l 83.094 35.551 85.426 33.07 + 90.094 33.07 c 96.047 33.07 l 96.047 30.664 l 96.047 29.676 95.742 28.91 + 95.141 28.367 c 94.535 27.816 93.641 27.539 92.453 27.539 c 85.016 27.539 + l 85.016 23.945 l 91.969 23.945 l 94.75 23.945 96.801 24.555 98.125 25.773 + c 99.457 26.984 100.125 28.938 100.125 31.633 c 100.125 43.148 l 100.125 + 43.629 100.305 43.988 100.672 44.227 c 101.047 44.469 101.664 44.586 102.531 + 44.586 c 102.531 47.945 l h +93.844 44.352 m 94.801 44.352 95.406 44.129 95.656 43.68 c 95.914 43.234 + 96.047 42.656 96.047 41.945 c 96.047 36.664 l 90.047 36.664 l 89.211 36.664 + 88.523 36.941 87.984 37.492 c 87.441 38.035 87.172 38.719 87.172 39.539 + c 87.172 41.945 l 87.172 42.781 87.363 43.391 87.75 43.773 c 88.133 44.16 + 88.742 44.352 89.578 44.352 c h +93.844 44.352 m f +114.047 44.352 m 114.047 27.539 l 110.547 27.539 l 110.547 23.945 l 117.891 + 23.945 l 117.891 25.242 l 118.953 24.379 120.234 23.945 121.734 23.945 +c 126.672 23.945 l 126.672 27.789 l 121.25 27.789 l 120.227 27.789 119.461 + 28.117 118.953 28.773 c 118.441 29.43 118.172 30.094 118.141 30.758 c 118.141 + 44.352 l 123.359 44.352 l 123.359 47.945 l 110.547 47.945 l 110.547 44.352 + l h +114.047 44.352 m f +146.359 47.945 m 144.723 47.945 143.473 47.422 142.609 46.367 c 141.754 + 45.305 141.328 43.91 141.328 42.18 c 141.328 27.539 l 136.766 27.539 l +136.766 23.945 l 141.328 23.945 l 141.328 18.18 l 145.406 18.18 l 145.406 + 23.945 l 152.406 23.945 l 152.406 27.539 l 145.406 27.539 l 145.406 42.43 + l 145.406 43.129 145.578 43.625 145.922 43.914 c 146.273 44.207 146.82 +44.352 147.563 44.352 c 152.406 44.352 l 152.406 47.945 l h +146.359 47.945 m f +172.719 14.352 m 172.719 18.664 l 168.156 18.664 l 168.156 14.352 l h +161.344 27.539 m 161.344 23.945 l 172.469 23.945 l 172.469 44.352 l 179.578 + 44.352 l 179.578 47.945 l 168.391 47.945 l 168.391 27.539 l h +161.344 27.539 m f +197.25 47.945 m 194.469 47.945 192.332 47.27 190.844 45.914 c 189.352 44.551 + 188.609 42.508 188.609 39.789 c 188.609 22.508 l 188.609 19.789 189.352 + 17.754 190.844 16.398 c 192.332 15.035 194.469 14.352 197.25 14.352 c 204.063 + 14.352 l 204.063 17.992 l 196.766 17.992 l 195.578 17.992 194.613 18.348 + 193.875 19.055 c 193.145 19.754 192.781 20.598 192.781 21.586 c 192.781 + 40.695 l 192.781 41.688 193.145 42.535 193.875 43.242 c 194.613 43.953 +195.578 44.305 196.766 44.305 c 204.063 44.305 l 204.063 47.945 l h +197.25 47.945 m f +0.831373 0 0 rg +17.277 53.707 m 17.277 56.629 14.91 58.996 11.988 58.996 c 9.066 58.996 + 6.695 56.629 6.695 53.707 c 6.695 50.785 9.066 48.414 11.988 48.414 c 14.91 + 48.414 17.277 50.785 17.277 53.707 c f +0 0 0.501961 rg +23.906 61.84 m 23.906 62.969 22.992 63.883 21.863 63.883 c 20.734 63.883 + 19.82 62.969 19.82 61.84 c 19.82 60.715 20.734 59.801 21.863 59.801 c 22.992 + 59.801 23.906 60.715 23.906 61.84 c f +0.596078 g +0.751181 w +0 J +0 j +[] 0.0 d +4 M q 1 0 0 1 0 0 cm +16.941 50.758 m 29.465 47.813 l S Q +27.848 47.258 m 30.453 47.578 l 28.262 49.023 l 28.523 48.41 28.355 47.699 + 27.848 47.258 c h +27.848 47.258 m f* +0.137102 w +1 j +q -1 0.235294 -0.235294 -1 0 0 cm +-15.851 -50.987 m -18.248 -51.872 l -15.849 -52.753 l -16.234 -52.23 -16.233 + -51.519 -15.851 -50.987 c h +-15.851 -50.987 m S Q +0.6 g +0.751178 w +0 j +q 1 0 0 1 0 0 cm +23.234 59.613 m 30.398 49.301 l S Q +28.824 49.969 m 30.973 48.461 l 30.313 51 l 30.098 50.375 29.496 49.957 + 28.824 49.969 c h +28.824 49.969 m f* +0.115667 w +1 j +q -0.694805 1 -1 -0.694805 0 0 cm +20.193 -42.855 m 18.17 -43.597 l 20.191 -44.342 l 19.87 -43.904 19.87 -43.302 + 20.193 -42.855 c h +20.193 -42.855 m S Q +Q Q +showpage +%%Trailer +end +%%EOF diff --git a/doc/user-manual/figures/logos/PPartiC.svg b/doc/user-manual/figures/logos/PPartiC.svg new file mode 100644 index 0000000..9966134 --- /dev/null +++ b/doc/user-manual/figures/logos/PPartiC.svg @@ -0,0 +1,159 @@ + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + PPartiC + + + + + + diff --git a/doc/user-manual/figures/logos/fpakc.eps b/doc/user-manual/figures/logos/fpakc.eps new file mode 100644 index 0000000..667a7f0 --- /dev/null +++ b/doc/user-manual/figures/logos/fpakc.eps @@ -0,0 +1,294 @@ +%!PS-Adobe-3.0 EPSF-3.0 +%%Creator: cairo 1.16.0 (https://cairographics.org) +%%CreationDate: Sat Nov 21 21:23:19 2020 +%%Pages: 1 +%%DocumentData: Clean7Bit +%%LanguageLevel: 2 +%%BoundingBox: 8 15 191 78 +%%EndComments +%%BeginProlog +50 dict begin +/q { gsave } bind def +/Q { grestore } bind def +/cm { 6 array astore concat } bind def +/w { setlinewidth } bind def +/J { setlinecap } bind def +/j { setlinejoin } bind def +/M { setmiterlimit } bind def +/d { setdash } bind def +/m { moveto } bind def +/l { lineto } bind def +/c { curveto } bind def +/h { closepath } bind def +/re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto + 0 exch rlineto 0 rlineto closepath } bind def +/S { stroke } bind def +/f { fill } bind def +/f* { eofill } bind def +/n { newpath } bind def +/W { clip } bind def +/W* { eoclip } bind def +/BT { } bind def +/ET { } bind def +/BDC { mark 3 1 roll /BDC pdfmark } bind def +/EMC { mark /EMC pdfmark } bind def +/cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def +/Tj { show currentpoint cairo_store_point } bind def +/TJ { + { + dup + type /stringtype eq + { show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse + } forall + currentpoint cairo_store_point +} bind def +/cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore + cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def +/Tf { pop /cairo_font exch def /cairo_font_matrix where + { pop cairo_selectfont } if } bind def +/Td { matrix translate cairo_font_matrix matrix concatmatrix dup + /cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point + /cairo_font where { pop cairo_selectfont } if } bind def +/Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def + cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def +/g { setgray } bind def +/rg { setrgbcolor } bind def +/d1 { setcachedevice } bind def +/cairo_data_source { + CairoDataIndex CairoData length lt + { CairoData CairoDataIndex get /CairoDataIndex CairoDataIndex 1 add def } + { () } ifelse +} def +/cairo_flush_ascii85_file { cairo_ascii85_file status { cairo_ascii85_file flushfile } if } def +/cairo_image { image cairo_flush_ascii85_file } def +/cairo_imagemask { imagemask cairo_flush_ascii85_file } def +%%EndProlog +%%BeginSetup +%%BeginResource: font f-0-0 +%!FontType1-1.1 f-0-0 1.0 +11 dict begin +/FontName /f-0-0 def +/PaintType 0 def +/FontType 1 def +/FontMatrix [0.001 0 0 0.001 0 0] readonly def +/FontBBox {30 -240 640 735 } readonly def +/Encoding 256 array +0 1 255 {1 index exch /.notdef put} for +dup 97 /a put +dup 99 /c put +dup 102 /f put +dup 107 /k put +dup 112 /p put +readonly def +currentdict end +currentfile eexec +f983ef0097ece636fb4a96c74d26ab84185f6dfa4a16a7a1c27bbe3f1156aea698df336d20b467 +b10e7f33846656653c5ac6962759d3056cbdb3190bac614b984bf5a132dc418192443014ba63de +800d392b6fea026574bb2535fd7bb5338f35bf15a88ea328fdaa49670c7852e3d060f3c5d6b07f +2ef6d0f22646c5d18e19a2ae3ee120390f6dd96f76dcf1e127de5e9299077a00c17c0d71e36e5b +9d5ec58fceda57739a6a4214d4b79d6c48d2784b60c320323c7acddddf34db833cac0cf109f799 +69d114a330d372e5c978a66acc84e3fe5557f6240856a013ffaa0199444e5c5036f775eba4a5c5 +8cde66cf604b9aca2178431127b8a1ff7ed633a65c04600af5f573483112251caad907dcd8c61e +23500065b1568be79a17d2379b63d656e8f7715cdfdf357f0e30d9ab91d113c88231d875d60897 +1dee27eb5da34a09d2f1b1a2e7ab8be1036393cd6ae53eff0d77f452e9bf45eccc9e4836ae77c5 +b74262fa43bfccfd9a147a18dcdae6e50cc518129a64f669a5fae69c8082dec571d7d01a4e9f05 +6d3e52de0ea004acdbd1b3bf9b19aa17f839a75365a77b7275442a967093ffdf1694a72f9978a2 +304ac331d401b48e3e62a3a92cd39516dd480f088980d1ad8993a1f6fefb07e0d86b6f0293bb41 +68ac465726267cacb7516a0e910fe67c2dbfef06d8b64a9811506650d32fa182a0adcab8e2e21e +ca6d0dc81959c25ea2d3f7ccec13e0cb4a7ef88e97c36e74fa13010220d6835ebdcbabdb507d84 +239e5483e8a8b7a52d6e1ea4ea1f5e6bef4534710c4055265aaa86fb445f3b2fc62cfdd9e283d8 +8bd083d09f0971cde00f2031b58b304d5f647f02aabf7ba9062c33979cd391f692c72ee179b7a8 +16f9c9e668d20021bd2b6a0f0114898729c6228be2895a696aaef0ebbcc842e64d5e72cc1d9b75 +44314028987a238f8fc4c18a0db3546c9ea42194b6bbdc45587e36d605fe2b7608d9292ddc0c9b +be3e420b36fd52f9aef97f13533e101f34d4f882848f4845a7a824da815a710abe11a1ec8363c1 +06daa18dffb5a451af7bf3b20d79a63c94e305050ea893b7a61cc8cfdb2e8a593a073c2021b298 +40863c70742ab1734e1d6811cf1927832da10f562f6895575b50044179588ddd9ec4c413f68c3e +3063f9594dc94115af4d9d4d6259c2ebb5afa796131772de3d297a8cc04f7f10398acc9142b1aa +2da9741ad314918ff1553dafe4751b4c0efdc9d6a549acbf1b3d209f6ebe8f6561d627f37bbce1 +7213b92bf332c27718ca9f868f1724cda0774ea4c3a5a2ba99509eb9128c456e5526f234dc3adc +37ac61ead9dafba1b5d58a9443ceb92474535cd3515e9ce357420b230fe927e81f06b2363c70aa +b6e00858a44972ad3f8759069235bba0b8ae2c65a59fe3ee5642f88a8550a765907eb4f9432ac4 +9e896114d0bc969bc2c14acd9a50c31e2095133b6b4fc11a1136dbba4b515eaabf0cba23ffe795 +3532a1fca89780a841f3a5fe2514d31fd6d41fbcb5b8caadd53c1fa7b06506963f37006269d0ae +fc1d5d6bd7f6788544e01a77bdf35aaacdd41dd4fc16237759516c60ee9b57e7b56606e0fb5a25 +9e6cb3b2e22d3ac4db73c228fda1b327cace7cbe25594ea2a9445efeda7826604e3daf18ca977a +9a2788cddcea95b5460b648c12bebc3c39302a07481fe18ceed4c9e12ac7d51f8104bd589cf9f6 +68371d3aa8510b5ff08c972b84fd4e48c356f9098b9d130085f59bc7e748ac59e9715060d29b38 +d828e479d8a85e0aaaa7b4ab34b547643b9c3417305b29be240ba9514cd4656079f26f378a6af9 +f7c694b1332a14baa203faa5c6745eea7c9cc0b3fc13f722ff941857d141db5665b32ada8b7831 +06f5c01361bd4b161450546ad8812d364b394658c4d3e1e2cf1daba58315787fbe299a27a387f2 +1545b5d4bc3022e6a6803e971698504ebc44c9f7c02256a5ae01c9fb3194baf05fb7c9cd862427 +95ca9cfd57ba76d6697fc842ee19ff35ce8b9d3299a185a4d28537058b5ea471a818f521183ffa +4b171f078577707cfbb3b2c6de4659b6dbacb952a0213b0586bd5879096ec53540ad07463e5b25 +1a38d715de6a0bb6f92d3e20e1b25be654fe7aad2e057428b68e2d6830f2ee0ca384fb862813ba +d48832ec89991fa5e6bf9d3135378f0000000000000000000000000000000000000000000000000000000000000000 +0000000000000000000000000000000000000000000000000000000000000000 +0000000000000000000000000000000000000000000000000000000000000000 +0000000000000000000000000000000000000000000000000000000000000000 +0000000000000000000000000000000000000000000000000000000000000000 +0000000000000000000000000000000000000000000000000000000000000000 +0000000000000000000000000000000000000000000000000000000000000000 +0000000000000000000000000000000000000000000000000000000000000000 +cleartomark +%%EndResource +%%EndSetup +%%Page: 1 1 +%%BeginPageSetup +%%PageBoundingBox: 8 15 191 78 +%%EndPageSetup +q 8 15 183 63 rectclip +1 0 0 -1 0 86 cm q +0 g +BT +56.000126 0 0 -56.000126 52.419201 52.791732 Tm +/f-0-0 1 Tf +[(fpak)10(c)]TJ +ET +0.749999 w +0 J +0 j +[] 0.0 d +4 M q 1 0 0 1 0 0 cm +8.996 27.238 m 20.621 9.359 l 34.656 9.363 l 34.621 14.746 l 25.941 18.195 + l 17.332 33.629 l 17.41 55.852 l 9.316 55.785 l h +8.996 27.238 m S Q +q 1 0 0 1 0 0 cm +21.156 69.262 m 21.156 38.387 l 40.305 38.387 l 40.305 53.027 l 27.457 +53.027 l 27.457 69.52 l h +21.156 69.262 m S Q +q 1 0 0 1 0 0 cm +27.898 43.988 m 27.898 47.73 l 33.426 47.73 l 33.426 43.898 l h +27.898 43.988 m S Q +0.283465 w +q 1 0 0 1 0 0 cm +9.094 50.535 m 17.414 50.297 l S Q +q 1 0 0 1 0 0 cm +9.34 43.07 m 17.094 44.93 l S Q +q 1 0 0 1 0 0 cm +9.43 34.727 m 17.348 41.363 l S Q +q 1 0 0 1 0 0 cm +17.332 33.629 m 8.996 27.238 l S Q +q 1 0 0 1 0 0 cm +13.297 31.055 m 11.32 36.062 l 12.719 43.309 l 14.023 50.113 l 12.02 55.516 + l S Q +q 1 0 0 1 0 0 cm +17.332 33.629 m 12.203 21.801 l S Q +q 1 0 0 1 0 0 cm +21.246 26.223 m 12.203 21.801 l S Q +q 1 0 0 1 0 0 cm +15.379 16.828 m 21.246 26.223 l S Q +q 1 0 0 1 0 0 cm +23.562 22.152 m 15.379 16.828 l S Q +q 1 0 0 1 0 0 cm +25.941 18.195 m 15.379 16.828 l S Q +q 1 0 0 1 0 0 cm +13.297 31.055 m 14.234 26.23 l S Q +q 1 0 0 1 0 0 cm +14.234 26.23 m 21.246 26.223 l S Q +q 1 0 0 1 0 0 cm +25.941 18.195 m 17.918 13.906 l S Q +q 1 0 0 1 0 0 cm +20.621 9.359 m 25.941 18.195 l S Q +q 1 0 0 1 0 0 cm +34.656 9.363 m 25.941 18.195 l S Q +q 1 0 0 1 0 0 cm +25.941 18.195 m 26.527 9.582 l S Q +1 0 0 rg +23.523 67.027 m 23.523 67.027 23.523 67.031 23.523 67.031 c 23.52 67.031 + 23.52 67.027 23.52 67.027 c 23.52 67.027 23.52 67.027 23.523 67.027 c h +23.523 67.027 m f +0 g +0.749999 w +q 1 0 0 1 0 0 cm +23.523 67.027 m 23.523 67.027 23.523 67.031 23.523 67.031 c 23.52 67.031 + 23.52 67.027 23.52 67.027 c 23.52 67.027 23.52 67.027 23.523 67.027 c h +23.523 67.027 m S Q +1 0 0 rg +24.66 63.297 m 24.66 63.297 24.656 63.297 24.656 63.297 c 24.656 63.293 + 24.66 63.293 24.66 63.293 c 24.66 63.293 24.66 63.293 24.66 63.297 c h +24.66 63.297 m f +0 g +q 1 0 0 1 0 0 cm +24.66 63.297 m 24.66 63.297 24.656 63.297 24.656 63.297 c 24.656 63.293 + 24.66 63.293 24.66 63.293 c 24.66 63.293 24.66 63.293 24.66 63.297 c h +24.66 63.297 m S Q +0 0 1 rg +23.797 64.039 m 23.797 64.273 23.609 64.465 23.371 64.465 c 23.137 64.465 + 22.949 64.273 22.949 64.039 c 22.949 63.805 23.137 63.613 23.371 63.613 + c 23.609 63.613 23.797 63.805 23.797 64.039 c h +23.797 64.039 m f +1 0 0 rg +37.844 42.566 m 37.844 42.945 37.535 43.254 37.156 43.254 c 36.773 43.254 + 36.465 42.945 36.465 42.566 c 36.465 42.184 36.773 41.875 37.156 41.875 + c 37.535 41.875 37.844 42.184 37.844 42.566 c h +37.844 42.566 m f +0 g +0.141732 w +q 1 0 0 1 0 0 cm +35.977 42.188 m 31.02 39.777 l S Q +31.531 40.027 m 31.91 39.895 l 30.895 39.715 l 31.66 40.406 l h +31.531 40.027 m f* +0.0679747 w +q 1 0.486446 -0.486446 1 0 0 cm +41.243 19.965 m 41.497 19.708 l 40.605 19.963 l 41.496 20.22 l h +41.243 19.965 m S Q +0.141732 w +q 1 0 0 1 0 0 cm +29.828 39.73 m 23.582 45.012 l S Q +24.016 44.645 m 24.047 44.246 l 23.473 45.102 l 24.414 44.68 l h +24.016 44.645 m f* +0.0577315 w +q 1 -0.845211 0.845211 1 0 0 cm +-8.002 37.881 m -7.787 37.664 l -8.544 37.88 l -7.787 38.098 l h +-8.002 37.881 m S Q +0.141732 w +q 1 0 0 1 0 0 cm +23.82 46.355 m 25.027 51.145 l S Q +24.887 50.594 m 24.543 50.387 l 25.059 51.281 l 25.094 50.25 l h +24.887 50.594 m f* +0.0732958 w +q -0.252177 -1 1 -0.252177 0 0 cm +-53.469 11.403 m -53.193 11.129 l -54.156 11.402 l -53.195 11.679 l h +-53.469 11.403 m S Q +0.141732 w +q 1 0 0 1 0 0 cm +23.453 63.156 m 26.16 58.465 l S Q +25.875 58.953 m 25.98 59.34 l 26.23 58.34 l 25.488 59.059 l h +25.875 58.953 m f* +0.0654792 w +q -0.576787 1 -1 -0.576787 0 0 cm +33.038 -44.931 m 33.282 -45.177 l 32.424 -44.932 l 33.284 -44.686 l h +33.038 -44.931 m S Q +0.141732 w +q 1 0 0 1 0 0 cm +26.215 57.066 m 26.281 52.215 l S Q +26.273 52.785 m 26.555 53.07 l 26.285 52.074 l 25.984 53.062 l h +26.273 52.785 m f* +0.0755832 w +q -0.0137977 1 -1 -0.0137977 0 0 cm +52.413 -26.997 m 52.694 -27.282 l 51.702 -26.999 l 52.694 -26.711 l h +52.413 -26.997 m S Q +0.141732 w +q 1 0 0 1 0 0 cm +26.5 50.707 m 34.664 48.957 l S Q +34.109 49.074 m 33.895 49.41 l 34.805 48.926 l 33.773 48.859 l h +34.109 49.074 m f* +0.0739077 w +q -1 0.214603 -0.214603 -1 0 0 cm +-22.54 -53.911 m -22.266 -54.188 l -23.235 -53.912 l -22.263 -53.637 l +h +-22.54 -53.911 m S Q +0.141732 w +q 1 0 0 1 0 0 cm +36.809 49.66 m 38.371 43.828 l S Q +38.223 44.379 m 38.426 44.727 l 38.406 43.691 l 37.879 44.578 l h +38.223 44.379 m f* +0.0730145 w +q -0.26796 1 -1 -0.26796 0 0 cm +31.85 -46.757 m 32.123 -47.034 l 31.163 -46.757 l 32.122 -46.486 l h +31.85 -46.757 m S Q +Q Q +showpage +%%Trailer +end +%%EOF diff --git a/doc/user-manual/figures/scatteringQuad.eps b/doc/user-manual/figures/scatteringQuad.eps new file mode 100644 index 0000000..ca8fd14 --- /dev/null +++ b/doc/user-manual/figures/scatteringQuad.eps @@ -0,0 +1,303 @@ +%!PS-Adobe-3.0 EPSF-3.0 +%%Creator: cairo 1.16.0 (https://cairographics.org) +%%CreationDate: Thu Oct 08 15:41:58 2020 +%%Pages: 1 +%%DocumentData: Clean7Bit +%%LanguageLevel: 2 +%%BoundingBox: 13 30 125 153 +%%EndComments +%%BeginProlog +50 dict begin +/q { gsave } bind def +/Q { grestore } bind def +/cm { 6 array astore concat } bind def +/w { setlinewidth } bind def +/J { setlinecap } bind def +/j { setlinejoin } bind def +/M { setmiterlimit } bind def +/d { setdash } bind def +/m { moveto } bind def +/l { lineto } bind def +/c { curveto } bind def +/h { closepath } bind def +/re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto + 0 exch rlineto 0 rlineto closepath } bind def +/S { stroke } bind def +/f { fill } bind def +/f* { eofill } bind def +/n { newpath } bind def +/W { clip } bind def +/W* { eoclip } bind def +/BT { } bind def +/ET { } bind def +/BDC { mark 3 1 roll /BDC pdfmark } bind def +/EMC { mark /EMC pdfmark } bind def +/cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def +/Tj { show currentpoint cairo_store_point } bind def +/TJ { + { + dup + type /stringtype eq + { show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse + } forall + currentpoint cairo_store_point +} bind def +/cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore + cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def +/Tf { pop /cairo_font exch def /cairo_font_matrix where + { pop cairo_selectfont } if } bind def +/Td { matrix translate cairo_font_matrix matrix concatmatrix dup + /cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point + /cairo_font where { pop cairo_selectfont } if } bind def +/Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def + cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def +/g { setgray } bind def +/rg { setrgbcolor } bind def +/d1 { setcachedevice } bind def +/cairo_data_source { + CairoDataIndex CairoData length lt + { CairoData CairoDataIndex get /CairoDataIndex CairoDataIndex 1 add def } + { () } ifelse +} def +/cairo_flush_ascii85_file { cairo_ascii85_file status { cairo_ascii85_file flushfile } if } def +/cairo_image { image cairo_flush_ascii85_file } def +/cairo_imagemask { imagemask cairo_flush_ascii85_file } def +%%EndProlog +%%BeginSetup +%%EndSetup +%%Page: 1 1 +%%BeginPageSetup +%%PageBoundingBox: 13 30 125 153 +%%EndPageSetup +q 13 30 112 123 rectclip +1 0 0 -1 0 171 cm q +0 0.784314 0 rg +20.836 84.285 52.273 43.371 re f +0 g +0.751178 w +0 J +0 j +[] 0.0 d +4 M q 1 0 0 1 0 0 cm +20.836 84.285 52.273 43.371 re S Q +1 0.627451 0 rg +73.109 84.285 43.184 43.371 re f +0 g +q 1 0 0 1 0 0 cm +73.109 84.285 43.184 43.371 re S Q +0.0431373 0 1 rg +73.109 32.199 43.184 52.086 re f +0 g +q 1 0 0 1 0 0 cm +73.109 32.199 43.184 52.086 re S Q +1 0 0 rg +20.836 32.199 52.273 52.086 re f +0 g +q 1 0 0 1 0 0 cm +20.836 32.199 52.273 52.086 re S Q +q 1 0 0 1 0 0 cm +20.836 32.199 95.457 95.457 re S Q +0 0 1 rg +22.902 127.77 m 22.902 128.914 21.973 129.844 20.824 129.844 c 19.68 129.844 + 18.75 128.914 18.75 127.77 c 18.75 126.621 19.68 125.691 20.824 125.691 + c 21.973 125.691 22.902 126.621 22.902 127.77 c h +22.902 127.77 m f +q 1 0 0 1 0 0 cm +22.902 127.77 m 22.902 128.914 21.973 129.844 20.824 129.844 c 19.68 129.844 + 18.75 128.914 18.75 127.77 c 18.75 126.621 19.68 125.691 20.824 125.691 + c 21.973 125.691 22.902 126.621 22.902 127.77 c h +22.902 127.77 m S Q +1 0 0 rg +118.059 127.23 m 118.059 128.379 117.129 129.309 115.98 129.309 c 114.836 + 129.309 113.906 128.379 113.906 127.23 c 113.906 126.086 114.836 125.156 + 115.98 125.156 c 117.129 125.156 118.059 126.086 118.059 127.23 c f +q 1 0 0 1 0 0 cm +118.059 127.23 m 118.059 128.379 117.129 129.309 115.98 129.309 c 114.836 + 129.309 113.906 128.379 113.906 127.23 c 113.906 126.086 114.836 125.156 + 115.98 125.156 c 117.129 125.156 118.059 126.086 118.059 127.23 c S Q +0 0.784314 0 rg +117.789 32.41 m 117.789 33.559 116.859 34.488 115.715 34.488 c 114.566 +34.488 113.637 33.559 113.637 32.41 c 113.637 31.266 114.566 30.336 115.715 + 30.336 c 116.859 30.336 117.789 31.266 117.789 32.41 c f +q 1 0 0 1 0 0 cm +117.789 32.41 m 117.789 33.559 116.859 34.488 115.715 34.488 c 114.566 +34.488 113.637 33.559 113.637 32.41 c 113.637 31.266 114.566 30.336 115.715 + 30.336 c 116.859 30.336 117.789 31.266 117.789 32.41 c S Q +1 0.627451 0 rg +23.238 32.41 m 23.238 33.559 22.309 34.488 21.16 34.488 c 20.016 34.488 + 19.086 33.559 19.086 32.41 c 19.086 31.266 20.016 30.336 21.16 30.336 c + 22.309 30.336 23.238 31.266 23.238 32.41 c f +q 1 0 0 1 0 0 cm +23.238 32.41 m 23.238 33.559 22.309 34.488 21.16 34.488 c 20.016 34.488 + 19.086 33.559 19.086 32.41 c 19.086 31.266 20.016 30.336 21.16 30.336 c + 22.309 30.336 23.238 31.266 23.238 32.41 c S Q +0 g +76.465 84.191 m 76.465 86.262 74.785 87.941 72.715 87.941 c 70.645 87.941 + 68.965 86.262 68.965 84.191 c 68.965 82.121 70.645 80.441 72.715 80.441 + c 74.785 80.441 76.465 82.121 76.465 84.191 c h +76.465 84.191 m f +q 1 0 0 1 0 0 cm +76.465 84.191 m 76.465 86.262 74.785 87.941 72.715 87.941 c 70.645 87.941 + 68.965 86.262 68.965 84.191 c 68.965 82.121 70.645 80.441 72.715 80.441 + c 74.785 80.441 76.465 82.121 76.465 84.191 c h +76.465 84.191 m S Q +0.749999 w +q 1 0 0 1 0 0 cm +20.355 84.105 m 115.98 84.105 l S Q +q 1 0 0 1 0 0 cm +72.922 127.277 m 72.922 32.008 l S Q +14.742 139.539 m 16.68 139.539 l 16.68 132.867 l 14.57 133.289 l 14.57 +132.211 l 16.664 131.789 l 17.852 131.789 l 17.852 139.539 l 19.789 139.539 + l 19.789 140.539 l 14.742 140.539 l h +14.742 139.539 m f +120.676 136.695 m 124.816 136.695 l 124.816 137.695 l 119.254 137.695 l + 119.254 136.695 l 119.699 136.238 120.309 135.617 121.082 134.836 c 121.863 + 134.047 122.352 133.535 122.551 133.305 c 122.934 132.891 123.199 132.535 + 123.348 132.242 c 123.504 131.941 123.582 131.648 123.582 131.367 c 123.582 + 130.898 123.414 130.52 123.082 130.227 c 122.758 129.938 122.336 129.789 + 121.816 129.789 c 121.441 129.789 121.043 129.852 120.629 129.977 c 120.223 + 130.102 119.785 130.301 119.316 130.57 c 119.316 129.367 l 119.793 129.18 + 120.238 129.039 120.645 128.945 c 121.059 128.844 121.441 128.789 121.785 + 128.789 c 122.691 128.789 123.414 129.02 123.957 129.477 c 124.496 129.926 + 124.77 130.531 124.77 131.289 c 124.77 131.645 124.699 131.984 124.566 +132.305 c 124.43 132.629 124.184 133.008 123.832 133.445 c 123.727 133.563 + 123.414 133.891 122.895 134.43 c 122.371 134.973 121.633 135.727 120.676 + 136.695 c h +120.676 136.695 m f +122.875 24.07 m 123.438 24.195 123.875 24.453 124.188 24.836 c 124.508 +25.211 124.672 25.68 124.672 26.242 c 124.672 27.109 124.375 27.781 123.781 + 28.258 c 123.188 28.727 122.344 28.961 121.25 28.961 c 120.883 28.961 120.504 + 28.922 120.109 28.852 c 119.723 28.781 119.328 28.672 118.922 28.523 c +118.922 27.383 l 119.242 27.57 119.598 27.719 119.984 27.82 c 120.379 27.914 + 120.789 27.961 121.219 27.961 c 121.957 27.961 122.52 27.816 122.906 27.523 + c 123.301 27.234 123.5 26.805 123.5 26.242 c 123.5 25.734 123.316 25.332 + 122.953 25.039 c 122.586 24.75 122.086 24.602 121.453 24.602 c 120.422 +24.602 l 120.422 23.633 l 121.5 23.633 l 122.07 23.633 122.516 23.52 122.828 + 23.289 c 123.141 23.051 123.297 22.711 123.297 22.273 c 123.297 21.828 +123.133 21.484 122.813 21.242 c 122.5 21.004 122.047 20.883 121.453 20.883 + c 121.117 20.883 120.766 20.922 120.391 20.992 c 120.023 21.055 119.617 + 21.16 119.172 21.305 c 119.172 20.258 l 119.629 20.133 120.051 20.039 120.438 + 19.977 c 120.832 19.914 121.203 19.883 121.547 19.883 c 122.453 19.883 +123.164 20.086 123.688 20.492 c 124.207 20.898 124.469 21.453 124.469 22.148 + c 124.469 22.641 124.328 23.051 124.047 23.383 c 123.773 23.719 123.383 + 23.945 122.875 24.07 c h +122.875 24.07 m f +17.223 19.746 m 14.238 24.418 l 17.223 24.418 l h +16.91 18.715 m 18.41 18.715 l 18.41 24.418 l 19.66 24.418 l 19.66 25.402 + l 18.41 25.402 l 18.41 27.465 l 17.223 27.465 l 17.223 25.402 l 13.285 +25.402 l 13.285 24.262 l h +16.91 18.715 m f +39.906 47.727 m 38.297 52.07 l 41.516 52.07 l h +39.234 46.555 m 40.578 46.555 l 43.906 55.305 l 42.672 55.305 l 41.875 +53.055 l 37.938 53.055 l 37.141 55.305 l 35.891 55.305 l h +39.234 46.555 m f +45.5 57.066 m 48.188 57.066 l 48.188 57.707 l 44.578 57.707 l 44.578 57.066 + l 44.867 56.766 45.266 56.359 45.766 55.848 c 46.266 55.34 46.582 55.012 + 46.719 54.863 c 46.957 54.582 47.125 54.348 47.219 54.16 c 47.32 53.965 + 47.375 53.777 47.375 53.598 c 47.375 53.297 47.266 53.051 47.047 52.863 + c 46.836 52.668 46.566 52.566 46.234 52.566 c 45.992 52.566 45.738 52.609 + 45.469 52.691 c 45.195 52.777 44.91 52.902 44.609 53.066 c 44.609 52.301 + l 44.922 52.176 45.207 52.082 45.469 52.02 c 45.738 51.957 45.988 51.926 + 46.219 51.926 c 46.801 51.926 47.27 52.074 47.625 52.363 c 47.977 52.656 + 48.156 53.047 48.156 53.535 c 48.156 53.777 48.109 54 48.016 54.207 c 47.93 + 54.418 47.773 54.66 47.547 54.941 c 47.484 55.016 47.281 55.23 46.938 55.582 + c 46.594 55.938 46.113 56.434 45.5 57.066 c h +45.5 57.066 m f +93.52 47.684 m 91.91 52.027 l 95.129 52.027 l h +92.848 46.512 m 94.191 46.512 l 97.52 55.262 l 96.285 55.262 l 95.488 53.012 + l 91.551 53.012 l 90.754 55.262 l 89.504 55.262 l h +92.848 46.512 m f +98.582 57.02 m 99.832 57.02 l 99.832 52.676 l 98.473 52.957 l 98.473 52.254 + l 99.832 51.973 l 100.598 51.973 l 100.598 57.02 l 101.848 57.02 l 101.848 + 57.66 l 98.582 57.66 l h +98.582 57.02 m f +45.98 98.063 m 44.371 102.406 l 47.59 102.406 l h +45.309 96.891 m 46.652 96.891 l 49.98 105.641 l 48.746 105.641 l 47.949 + 103.391 l 44.012 103.391 l 43.215 105.641 l 41.965 105.641 l h +45.309 96.891 m f +53.23 104.98 m 53.605 105.055 53.895 105.215 54.105 105.465 c 54.313 105.715 + 54.418 106.023 54.418 106.387 c 54.418 106.949 54.219 107.387 53.824 107.699 + c 53.438 108.004 52.891 108.152 52.184 108.152 c 51.941 108.152 51.699 +108.125 51.449 108.074 c 51.199 108.035 50.938 107.965 50.668 107.871 c +50.668 107.137 l 50.875 107.254 51.105 107.348 51.355 107.418 c 51.613 107.48 + 51.887 107.512 52.168 107.512 c 52.645 107.512 53.012 107.418 53.262 107.23 + c 53.52 107.035 53.652 106.754 53.652 106.387 c 53.652 106.055 53.531 105.793 + 53.293 105.605 c 53.051 105.418 52.723 105.324 52.309 105.324 c 51.652 +105.324 l 51.652 104.684 l 52.34 104.684 l 52.723 104.684 53.016 104.613 + 53.215 104.465 c 53.41 104.309 53.512 104.09 53.512 103.809 c 53.512 103.52 + 53.406 103.293 53.199 103.137 c 52.988 102.98 52.691 102.902 52.309 102.902 + c 52.098 102.902 51.875 102.93 51.637 102.98 c 51.395 103.023 51.129 103.09 + 50.84 103.184 c 50.84 102.496 l 51.129 102.414 51.402 102.355 51.652 102.324 + c 51.91 102.285 52.156 102.262 52.387 102.262 c 52.969 102.262 53.426 102.398 + 53.762 102.668 c 54.105 102.93 54.277 103.285 54.277 103.73 c 54.277 104.043 + 54.184 104.309 53.996 104.527 c 53.816 104.746 53.563 104.898 53.23 104.98 + c h +53.23 104.98 m f +93.707 97.875 m 92.098 102.219 l 95.316 102.219 l h +93.035 96.703 m 94.379 96.703 l 97.707 105.453 l 96.473 105.453 l 95.676 + 103.203 l 91.738 103.203 l 90.941 105.453 l 89.691 105.453 l h +93.035 96.703 m f +100.754 102.836 m 98.801 105.867 l 100.754 105.867 l h +100.551 102.164 m 101.52 102.164 l 101.52 105.867 l 102.316 105.867 l 102.316 + 106.508 l 101.52 106.508 l 101.52 107.852 l 100.754 107.852 l 100.754 106.508 + l 98.176 106.508 l 98.176 105.773 l h +100.551 102.164 m f +27.949 77.395 m 27.949 80.301 l 27.043 80.301 l 27.043 72.754 l 27.949 +72.754 l 27.949 73.582 l 28.137 73.262 28.371 73.02 28.652 72.863 c 28.941 + 72.707 29.293 72.629 29.699 72.629 c 30.363 72.629 30.902 72.895 31.309 + 73.426 c 31.723 73.949 31.934 74.637 31.934 75.488 c 31.934 76.355 31.723 + 77.051 31.309 77.582 c 30.902 78.105 30.363 78.363 29.699 78.363 c 29.293 + 78.363 28.941 78.285 28.652 78.129 c 28.371 77.973 28.137 77.73 27.949 +77.395 c h +31.012 75.488 m 31.012 74.832 30.871 74.316 30.59 73.941 c 30.316 73.566 + 29.949 73.379 29.48 73.379 c 29 73.379 28.625 73.566 28.355 73.941 c 28.082 + 74.316 27.949 74.832 27.949 75.488 c 27.949 76.156 28.082 76.676 28.355 + 77.051 c 28.625 77.426 29 77.613 29.48 77.613 c 29.949 77.613 30.316 77.426 + 30.59 77.051 c 30.871 76.676 31.012 76.156 31.012 75.488 c h +31.012 75.488 m f +35.578 70.629 m 35.148 71.379 34.828 72.125 34.609 72.863 c 34.398 73.594 + 34.297 74.332 34.297 75.082 c 34.297 75.832 34.398 76.578 34.609 77.316 + c 34.828 78.059 35.148 78.793 35.578 79.535 c 34.797 79.535 l 34.316 78.773 + 33.953 78.027 33.703 77.285 c 33.461 76.547 33.344 75.813 33.344 75.082 + c 33.344 74.355 33.461 73.625 33.703 72.895 c 33.941 72.156 34.305 71.402 + 34.797 70.629 c h +35.578 70.629 m f +41.875 72.754 m 39.906 75.41 l 41.984 78.223 l 40.922 78.223 l 39.328 76.066 + l 37.734 78.223 l 36.672 78.223 l 38.797 75.363 l 36.859 72.754 l 37.922 + 72.754 l 39.375 74.707 l 40.813 72.754 l h +41.875 72.754 m f +43.121 79.676 m 44.168 79.676 l 44.168 76.066 l 43.027 76.301 l 43.027 +75.707 l 44.152 75.488 l 44.793 75.488 l 44.793 79.676 l 45.84 79.676 l +45.84 80.223 l 43.121 80.223 l h +43.121 79.676 m f +47.621 76.988 m 48.652 76.988 l 48.652 77.816 l 47.855 79.379 l 47.215 +79.379 l 47.621 77.816 l h +47.621 76.988 m f +55.117 72.754 m 53.148 75.41 l 55.227 78.223 l 54.164 78.223 l 52.57 76.066 + l 50.977 78.223 l 49.914 78.223 l 52.039 75.363 l 50.102 72.754 l 51.164 + 72.754 l 52.617 74.707 l 54.055 72.754 l h +55.117 72.754 m f +56.801 79.676 m 59.035 79.676 l 59.035 80.223 l 56.02 80.223 l 56.02 79.676 + l 56.27 79.426 56.602 79.094 57.02 78.676 c 57.434 78.25 57.699 77.973 +57.816 77.848 c 58.023 77.621 58.164 77.426 58.238 77.27 c 58.32 77.105 +58.363 76.941 58.363 76.785 c 58.363 76.535 58.273 76.332 58.098 76.176 +c 57.918 76.02 57.691 75.941 57.41 75.941 c 57.211 75.941 56.996 75.98 56.77 + 76.051 c 56.551 76.113 56.316 76.219 56.066 76.363 c 56.066 75.707 l 56.316 + 75.605 56.551 75.527 56.77 75.473 c 56.996 75.422 57.207 75.395 57.395 +75.395 c 57.883 75.395 58.273 75.52 58.566 75.77 c 58.855 76.012 59.004 +76.34 59.004 76.754 c 59.004 76.941 58.965 77.125 58.895 77.301 c 58.82 +77.48 58.691 77.688 58.504 77.926 c 58.449 77.988 58.277 78.168 57.988 78.457 + c 57.707 78.75 57.309 79.156 56.801 79.676 c h +56.801 79.676 m f +60.484 70.629 m 61.266 70.629 l 61.754 71.402 62.117 72.156 62.359 72.895 + c 62.609 73.625 62.734 74.355 62.734 75.082 c 62.734 75.813 62.609 76.547 + 62.359 77.285 c 62.117 78.027 61.754 78.773 61.266 79.535 c 60.484 79.535 + l 60.922 78.793 61.242 78.059 61.453 77.316 c 61.672 76.578 61.781 75.832 + 61.781 75.082 c 61.781 74.332 61.672 73.594 61.453 72.863 c 61.242 72.125 + 60.922 71.379 60.484 70.629 c h +60.484 70.629 m f +Q Q +showpage +%%Trailer +end +%%EOF diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index 5321ba4..c099b01 100644 Binary files a/doc/user-manual/fpakc_UserManual.pdf and b/doc/user-manual/fpakc_UserManual.pdf differ diff --git a/doc/user-manual/fpakc_UserManual.tex b/doc/user-manual/fpakc_UserManual.tex index 845cc9a..fccf4a0 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -1,5 +1,5 @@ \documentclass[10pt,a4paper,twoside]{book} -%\usepackage[latin1]{inputenc} +\usepackage[latin1]{inputenc} \usepackage{amsmath} \usepackage{amsfonts} \usepackage{amssymb} @@ -72,7 +72,7 @@ The \Gls{fpakc} is a simulation tool that models species in plasma (ions, electrons and neutrals) following the trajectories of macro-particles as they move and interact between them and the boundaries of the domain. Particles properties are scattered into a finite element mesh in 1, 2 or three dimensions, with the possibility to choose different geometries. The official repository can be found at: \url{https://gitlab.com/JorgeGonz/fpakc.git}. - The code is currently in the very early steps of development and further refinements are expected very soon. + The code is currently in very early steps of development and further improvements are expected very soon. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Main Guidelines} @@ -86,11 +86,11 @@ \item \acrshort{fpakc} is coded in a \textit{understandable} way. This means that the code is required to be written in a clear way that is easy to understand and maintain. Variables and procedure names need to be self-understanding. - This eases the process of fixing bugs and improving the codes by a large team of developers. + This ease the process of fixing bugs and improving the codes by a large team of developers. For more information, please refer to the \acrshort{fpakc} Coding Style document. - \item \acrshort{fpakc} requires being ease to use. + \item \acrshort{fpakc} requires to be ease to use. Input files are required to be in a \textit{human} format, meaning that the different options can be easily understood without constant reference to the user guide. - \acrshort{fpakc} is aimed to be used in a wide range of applications and by various scientists: from well-established ones to newcomers to the field and also students. + \acrshort{fpakc} is aimed to be used in a wide range of applications and by a variety of scientist: from very established ones to newcomers to the field and also students. \end{enumerate} These are foundation stones of \acrshort{fpakc} and its development and should always be followed, at least for the releases in the official repository. @@ -105,16 +105,16 @@ \section{The Particle Method} \Gls{fpakc} uses macro-particles to simulate the dynamics of different plasma species (mainly ions, electrons and neutrals). These macro-particles could represent a large amount of real particles. - For now own, macro-particles will be referred as just particles by abuse of language. + For now own, macro-particles will be referred as just particles by abusing of language. During the initiation phase, the input and mesh file(s) are reading. If an initial distribution for a species is specified in the input file, particles to match that distribution are loaded into the cells. The general steps performed in each iteration are: \begin{enumerate} \item Firstly, new particles are introduced into the domain as specified in the input file. - \item Particles are then pushed, accounting for possible acceleration by external forces. - During this process, if a particle changes cell, it is found using the connectivity between elements. - If a particle encounters a boundary instead a new cell, the interaction between the boundary and the wall are computed. + \item Particles are then pushed accounting for possible acceleration by external forces. + During this process, if a particle changes cell it is found using the connectivity between elements. + If a particle encounters a boundary instead a new cell, the interaction between the boundary and the wall is computed. A particle may abandon the computational domain and is no longer accounted for. \item Next, collisions for the particles inside each cell are carried out. This may include different collision processes for each particle. @@ -124,10 +124,10 @@ \item Finally, particle properties are scattered among the mesh nodes. These properties are density, momentum and the stress tensor. \item If requested, the electromagnetic field is computed. - \item If the number of iteration requires writing output files, it is done after all steps for the particles are completed. + \item If the number of iteration requires writing output files, it is done after all steps for the particles is completed. \end{enumerate} - \Gls{fpakc} has the capability to configure all the behaviour of the simulation via the input file. + \Gls{fpakc} has the capability to configure all the behavior of the simulation via the input file. Parameters as injection, the kind of pusher used for each species, boundary conditions or collisions are user-input parameters and will be described in Chap.~\ref{ch:input_file}. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -168,8 +168,8 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Find new cell} Once the position and velocity of the particle are updated, the new cell that contains the particle is searched. - This is done by a neighbour search, starting from the previous cell containing the particle. - In the process of finding the new cell, a particle might encounter a boundary. + This is done by a neighbor search, starting from the previous cell containing the particle. + In the process of finding the new cell, it is possible that a particle encounters a boundary. When the particle interacts with the boundary, the particle may continue its life in the simulation or might be eliminated from it. Once that the new cell is found or that the particle life has been terminated, the pushing is complete. If a secondary mesh is used for the Monte-Carlo Collision method, the new cell in that mesh in which the particle reside is also found by the same method, although no interaction with the boundaries is accounted for this step. @@ -178,7 +178,7 @@ \section{Variable Weighting Scheme\label{sec:weightingScheme}} One of the issues in particle simulations, specially for axial-symmetrical cases, is that due to the disparate volume of cells, specially close to the axis, the statistics in some cells is usually poor. To try to fix that, the possibility to include a Variable Weighting Scheme in the simulations is available in \Gls{fpakc}. - These schemes detect when a particle changes cells and split it if necessary to improve statistics. + These schemes detect when a particle change cells and split it if necessary to improve statistics. The use of a Variable Weighting Scheme is defined by the user in the input file. Beware that this can increase the number of particles in the simulation and increase computational time. @@ -189,16 +189,16 @@ \Gls{fpakc} distinguish between two types of interactions: \acrfull{mcc} and \acrfull{cs}. \acrshort{mcc} refers to the process in which two particles interact in short range. These processes include, but are not limited to: elastic collisions, ionization/recombination, charge-exchange, excitation/de-excitation\ldots - A secondary mesh, with cell sizes in the range of the mean-free path, can be used for this type of collision. + A secondary mesh, with cell sizes in the range of the mean-free path, can be used for this type of collisions. \acrshort{cs} refers to the large range interaction that a charged species suffer do to the charge of other particles. The interactions between the different species is defined by the user. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{\acrlong{mcc}} - For each cell, the maximum number of collisions between particle is computed. + For each cell the maximum number of collisions between particle is computed. For each collision, a random pair of particles is chosen. A loop over all possible collisions for the pair of particles chosen is performed. - If a random number is above the probability of collision for that specific type, the collision takes place. + If a random number is above the probability of collision for that specific type, the collision take place. If not, the next type for the particle pair is checked. Below are described the type of collision process implemented in \acrshort{fpakc}: @@ -219,7 +219,7 @@ \item Recombination. When an electron and an ion interact, there is a possibility for them to be recombined into a neutral particle. - The photons emitted by this process are not modelled yet. + The photon emitted by this process is not modelled yet. \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -238,7 +238,7 @@ Once that the pushing is complete, the array of particles that remain inside the domain is copied to a new array. The new array containing only the particles inside the domain will be the one used in the next steps. In this section, particles are assigned to the list of particles inside each individual cell. - Unfortunately, this is done right now without parallelization and is very CPU consuming. + Unfortunately, this is done right now without parallelisation and is very CPU consuming. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Probing}\label{sec:probing} @@ -250,21 +250,21 @@ The user can decide the grid width and the number of points in each direction. The distribution function will be calculated and wrote with a time step decided by the user. - If a particle velocity resides outside the velocity grid (in any direction), it will not be added to the tally of the distribution function. - Due to the limitation of only considering particles in the cell, and not neighbour particles, two probes for the same species at different positions but in the same cell will output the same results. - A more advance method considering distance between the particles and the probe position as well as particles in neighbour cells could be implemented to improve the statistics of the distribution function. + If a particle velocity resides outside of the velocity grid (in any direction), it wont be added to the tally of the distribution function. + Due to the limitation of only taking into account particles in the cell, and not neighbour particles, two probes for the same species at different positions but in the same cell will output the same results. + A more advance method taking into account distance between the particles and the probe position as well as particles in neighbour cells could be implemented to improve the statistics of the distribution function. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Scattering} The properties of each particle are deposited in the nodes from the containing cell. - This process depends on the cell type, but in general, each node receives a proportional part of the particle properties as a function of the particle position inside the cell. - The figure \ref{fig:scatteringQuad} shows how a particle at a generic position $p(x_1, x_2)$ inside the cell is scattered to the four nodes. + This process depend on the cell type, but in general, each node receive a proportional part of the particle properties as a function of the particle position inside the cell. + Figure \ref{fig:scatteringQuad} shows how a particle at a generic position $p(x_1, x_2)$ inside the cell is scattered to the four nodes. \begin{wrapfigure}{l}{0.4\textwidth} \centering \includegraphics{figures/scatteringQuad} \caption{\label{fig:scatteringQuad}Example of how a particle is weighted in a quadrilateral cell.} \end{wrapfigure} - Each node receives a proportional part of the area formed by dividing the cell in for rectangles, using as an additional vertex the particle position. + Each node receives a proportional part of the area formed by dividing the cell in for rectangles using as an additional vertex the particle position. These properties are dimensionless, but they are converted to the correct units once the output is printed. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -273,11 +273,11 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Average scheme} - Particle-in-cell codes have an intrinsic statistical noise associated with them. + Particle-in-cell codes has an intrinsic statistical noise associated with them. Although this can be reduced by increasing the number of particles, this also increases the CPU requirements of the case. It is quite common that most cases reach a quasi-steady state after a number of iterations and time-average results can be obtained after to improve analysis, plotting and restarting the case using these time-average results as new species backgrounds. - Although this is possible to do once the simulation is finished with post-processing tools, this is limited to the number of iterations printed. + Although this is possible to do once the simulation is finished with post-processing tools, this is limited to the amount of iterations printed. \Gls{fpakc} implements a simple average scheme that, after a start time provided by the user, scores a mean and standard deviation of all the main species properties, and the electromagnetic field. This scheme is based on the Welford's online algorithm~\cite{welford1962note}. @@ -286,7 +286,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Installation} \section{Required Packages} - To properly compile \gls{fpakc}, the following packages are required. + In order to properly compile \gls{fpakc}, the following packages are required. \subsection{Gfortran} The \Gls{opensource} free compiler \Gls{gfortran}\cite{gfortranURL} from GCC is the basic way to compile \acrshort{fpakc}. It is distributed with all GNU/Linux distributions. @@ -369,7 +369,7 @@ make %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Case file} The required format for the case file is \Gls{json}. - \Gls{json} is a case-sensitive format, so input must be written with the correct capitalization. + \Gls{json} is a case-sensitive format, so input must be written with the correct capitalisation. The basic structure and options available for the case file are explained below. The order of the objects and variables is irrelevant, but the structure needs to be maintained. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -380,9 +380,9 @@ make \item \textbf{path}: Character. Path for the output files. This path is also used to locate the mesh input file. \item \textbf{folder}: Character. - Base name of the folder in which output files are placed. + Base name of the folder in wich output files are placed. The date and time is appended to this name. - If none is provided, only the date and time is written as the folder name. + If none is provided, only the date and time is writted as the folder name. \item \textbf{triggerOutput}: Integer. Determines the number of iterations between writing output files for macroscopic quantities. \item \textbf{cpuTime}: Logical. @@ -460,10 +460,6 @@ make \begin{itemize} \item \textbf{gmsh2}: \Gls{gmsh} file format in version 2.0. This has to be in ASCII format. \item \textbf{vtu}: \Gls{vtu} file format. This has to be in ASCII format. - \item \textbf{text}: Plain text file format only intended for 1D cases. - This has to be in ASCII format and comma separated. - The first column represents the position and the second column the physical ID of the node. - Values have to be $1$ (left boundary), $2$ (right boundary), or $0$ (no boundary.) \end{itemize} \item \textbf{meshFile}: Character. Mesh filename. @@ -519,7 +515,7 @@ make \item \textbf{absorption}: Particle is eliminated from the domain. The particle is first moved into the edge and its properties are scattered among the edge nodes. \item \textbf{transparent}: Particle abandon the numerical domain. - \item \textbf{wallTemperature}: Reflective wall with constant temperature that exchange heat with particles. + \item \textbf{wallTemperature}: Reflective wall with cosntant temperature that exchange heat with particles. Required parameters are: \begin{itemize} \item \textbf{temperature}: Real. @@ -530,8 +526,8 @@ make Specific heat capacity of the material. \end{itemize} \item \textbf{ionization}: Per each particle crossing the surface with this type of boundary, a number of ionization events are calculated. - A pair of ion-electron is generated for each ionization event, taking as a reference a neutral background. - The secondary electron is taken as the same type as the incident particle. + A pair of ion-electron is generated for each ionization event taking as a reference a neutral background. + Secondary electron is taken as same type as incident particle. The available input is: \begin{itemize} \item \textbf{neutral}: Object. @@ -544,7 +540,7 @@ make \item \textbf{mass}: Real. Units in $\unit{kg}$. Mass of neutral species. - If missing, the mass of the ion is used + If missing, the mass of the ion is ussed \item \textbf{density}: Real. Units in $\unit{m^{-3}}$. Density of neutral background. @@ -562,18 +558,18 @@ make \end{itemize} \item \textbf{effectiveTime}: Real. Units in $\unit{s}$. - As the particle is no longer simulated once it crossed the boundary, this time represents the effective time in which the particle produces ionization processes in the neutral background. + As the particle is no longer simulated once it crossed the boundary, this time represent the effective time in which the particle produces ionization processes in the neutral background. Required parameter. \item \textbf{energyThreashold}: Real. Units in $\unit{eV}$. Ionization energy threshold for the simulated process. Required parameter. \item \textbf{crossSection}: Character. - Complete path to the cross-section data for the ionization process. + Complete path to the cross section data for the ionization process. \end{itemize} \item \textbf{axis}: Identifies the symmetry axis for 2D cylindrical simulations. - If , for some reason, a particle interacts with this axis, it is reflected. + If for some reason a particle interact with this axis, it is reflected. \end{itemize} \end{itemize} @@ -589,26 +585,18 @@ make Type of boundary. Accepted values are: \begin{itemize} - \item \textbf{dirichlet}: Constant value of electric potential on the surface. - \item \textbf{dirichletTime}: Constant value of the electric potential with a time variable profile. - The value of \textbf{boundaryEM.potential} will be multiplied for the corresponding value in the file \textbf{boundaryEM.temporalProfile}. + \item \textbf{dirichlet}: Elastic reflection of particles. \end{itemize} - \item \textbf{potential}: Real. - Fixed potential for Dirichlet boundary condition. + \item \textbf{potential}: Real. + Fixed potential for Dirichlet boundary condition. \item \textbf{physicalSurface}: Integer. Identification of the edge in the mesh file. - \item \textbf{temporalProfile}: Character. - Filename of the 2 column file containing the time variable profile. - File must be located in \textbf{output.path}. - The first column is the time in $\unit{s}$. - The second column is the factor that will multiply the value of the boundary. - \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{inject} The array \textbf{inject} specifies the injection of particles from different surfaces. - The injection of particles needs to be associated to a physicalSurface in the mesh file. + The injection of particles need to be associated to a physicalSurface in the mesh file. Multiple injections can be associated to the same surface. \begin{itemize} \item \textbf{name}: Character. @@ -622,9 +610,7 @@ make Available values are: \begin{itemize} \item \textbf{A}: Ampere. - \item \textbf{Am2}: Ampere per square meter. - This value will be multiplied by the area of injection. - \item \textbf{sccm}: Standard cubic centimetre. + \item \textbf{sccm}: Standard cubic centimeter. \item \textbf{part/s}: Particles (real) per second. \end{itemize} \item \textbf{v}: Real. @@ -641,7 +627,7 @@ make \begin{itemize} \item \textbf{Maxwellian}: Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction. \item \textbf{Half-Maxwellian}: Half-Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction. - Only considers the positive part of the half-Maxwellian. + Only takes into account the positive part of the half-Maxwellian. \item \textbf{Delta}: Dirac's delta distribution function. All particles are injected with velocity \textbf{v} times the value of \textbf{n} in the specified direction. \end{itemize} \item \textbf{T}: Real. @@ -650,11 +636,6 @@ make Temperature in each direction. \item \textbf{physicalSurface}: Integer. Identification of the edge in the mesh file. - \item \textbf{particlesPerEdge}: Integer. - Optional. - Number of particles to be injected by each edge in the numerical domain. - The weight of the particles for each edge will modified by the surface of the edge to ensure the right flux is injected. - If no value is provided, the number of particles to inject per edge will be calculated with the species weight and the surface of the edge respect to the total one. \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{reference} @@ -670,7 +651,7 @@ make \item \textbf{radius}: Real. Reference atomic radius in $\unit{m}$. \item \textbf{crossSection}: Real. - Reference cross-section in $\unit{m^2}$. + Reference cross section in $\unit{m^2}$. If this value is present, radius is ignored. \end{itemize} @@ -696,8 +677,8 @@ make Indicates the type of pusher used for each species: \begin{itemize} \item \textbf{Neutral}: Pushes a particle without any external force. - \item \textbf{Electrostatic}: Pushes a particle, including the effect of the electrostatic field. - \item \textbf{Electromagnetic}: Pushes a particle, accounting for the electromagnetic field. + \item \textbf{Electrostatic}: Pushes a particle including the effect of the electrostatic field. + \item \textbf{Electromagnetic}: Pushes particles accounting for the electromagnetic field. \end{itemize} \item \textbf{WeightingScheme}: Character. Indicates the variable weighting scheme to be used in the simulation. @@ -725,15 +706,11 @@ make \begin{itemize} \item \textbf{species}: Character. Name of species as defined in the object \textbf{species}. - \item \textbf{file}: Character. - Output file from previous run used as an initial state for the species. - The file format must be the same as in \textbf{geometry.meshType} - Initial particles are assumed to have a Maxwellian distribution. - File must be located in \textbf{output.path}. - \item \textbf{particlesPerCell}: Integer. - Optional. - Initial number of particles per cell. - If not, the number of particles per cell will be assigned based on the species weight and the cell volume. + \item \textbf{file}: Character. + Output file from previous run used as an initial state for the species. + The file format must be the same as in \textbf{geometry.meshType} + Initial particles are assumed to have a Maxwellian distribution. + File must be located at \textbf{output.path}. \end{itemize} \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -749,11 +726,11 @@ make \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{interactions}\label{ssec:input_interactions} - This object determines the different interactions among species. + This object determine the different interactions among species. Acceptable values are: \begin{itemize} \item \textbf{folderCollisions}: Character. - Indicates the path to in which the cross-section tables are allocated. + Indicates the path to in which the cross section tables are allocated. \item \textbf{meshCollisions}: Character. Determines a specific mesh for \acrshort{mcc} processes. The file needs to be located in the folder \textbf{output.folder}. @@ -780,18 +757,13 @@ make Accepted values are \textbf{elastic}, \textbf{chargeExchange}, \textbf{ionization} and \textbf{recombination}. Please refer to Sec.~\ref{ssec:collisions} for a description of the different collision types. \item \textbf{crossSection}: Character. - File in \textbf{interactions.folderCollisions} that contains the cross-section data as a 1D table of relative energy (in $\unit{eV}$) and cross-section (in $\unit{m^-2}$). + File in \textbf{interactions.folderCollisions} that contains the cross section data as a 1D table of relative energy (in $\unit{eV}$) and cross section (in $\unit{m^-2}$). \item \textbf{energyThreshold}: Real. Energy threshold of the collisional process in $\unit{eV}$. Only valid for \textbf{ionization} and \textbf{recombination} processes. \item \textbf{electron}: Character. Name of species designed as electrons. Only valid for \textbf{ionization} and \textbf{recombination} processes. - \item \textbf{electronSecondary}: Character. - Optional. - Name of species designed as secondary electrons. - If none provided, \textbf{electron} is used. - Only valid for \textbf{ionization}. \end{itemize} \end{itemize} \item \textbf{Coulomb}: Array of objects. @@ -801,7 +773,7 @@ make \begin{itemize} \item \textbf{species\_i}, \textbf{species\_j}: Character. Define the two species involved in the collision processes. - Order is indifferent. + Order is indiferent. \end{itemize} \end{itemize} @@ -827,9 +799,9 @@ make %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{1D Emissive Cathode (1D\_Cathode)} - Emission from a 1D cathode in both, Cartesian and radial coordinates. - Both cases insert the same number of electrons from the minimum coordinate and have the same boundary conditions for particles and the electrostatic field. - This case is useful to illustrate how \acrshort{fpakc} can deal with different geometries by just modifying some parameters in the input file. + Emission from a 1D cathode in both, cartesian and radial coordinates. + Both cases insert the same amount of electrons from the minimum coordinate and have the same boundary conditions for particles and the electrostatic field. + This case is useful to ilustrate hoy \acrshort{fpakc} can deal with different geometries by just modifying some parameters in the input file. The same mesh file (\lstinline|mesh.msh|) is used for both cases. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/src/fpakc.f90 b/src/fpakc.f90 index ae6e7bb..e85f5f6 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -11,12 +11,12 @@ PROGRAM fpakc USE OMP_LIB IMPLICIT NONE + ! t = time step + INTEGER:: t ! arg1 = Input argument 1 (input file) CHARACTER(200):: arg1 ! inputFile = path+name of input file CHARACTER(:), ALLOCATABLE:: inputFile - ! generic integer for time step - INTEGER:: t tStep = omp_get_wtime() !Gets the input file @@ -27,18 +27,11 @@ PROGRAM fpakc !Reads the json configuration file CALL readConfig(inputFile) - - !Create output folder and initial files - CALL initOutput(inputFile) - !Do '0' iteration - timeStep = tInitial + t = tInitial !$OMP PARALLEL DEFAULT(SHARED) !$OMP SINGLE - ! Initial reset of probes - CALL resetProbes() - CALL verboseError("Initial scatter of particles...") !$OMP END SINGLE CALL doScatter() @@ -52,21 +45,19 @@ PROGRAM fpakc tStep = omp_get_wtime() - tStep !Output initial state - CALL doOutput() + CALL doOutput(t) CALL verboseError('Starting main loop...') !$OMP PARALLEL DEFAULT(SHARED) DO t = tInitial + 1, tFinal + !Insert new particles and push them !$OMP SINGLE tStep = omp_get_wtime() - ! Update global time step index - timeStep = t - !Checks if a species needs to me moved in this iteration - CALL solver%updatePushSpecies() + CALL solver%updatePushSpecies(t) !Checks if probes need to be calculated this iteration - CALL resetProbes() + CALL resetProbes(t) tPush = omp_get_wtime() !$OMP END SINGLE @@ -84,7 +75,7 @@ PROGRAM fpakc !$OMP END SINGLE IF (doMCCollisions) THEN - CALL meshForMCC%doCollisions() + CALL meshForMCC%doCollisions(t) END IF @@ -129,12 +120,12 @@ PROGRAM fpakc !$OMP SINGLE tEMField = omp_get_wtime() - tEMField - CALL doAverage() + CALL doAverage(t) tStep = omp_get_wtime() - tStep !Output data - CALL doOutput() + CALL doOutput(t) !$OMP END SINGLE END DO diff --git a/src/makefile b/src/makefile index 83cf8db..247c7d2 100644 --- a/src/makefile +++ b/src/makefile @@ -9,7 +9,6 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module $(OBJDIR)/moduleMeshInputVTU.o $(OBJDIR)/moduleMeshOutputVTU.o \ $(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \ $(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \ - $(OBJDIR)/moduleMeshInputText.o $(OBJDIR)/moduleMeshOutputText.o \ $(OBJDIR)/moduleMesh3DCart.o \ $(OBJDIR)/moduleMesh2DCyl.o \ $(OBJDIR)/moduleMesh2DCart.o \ diff --git a/src/modules/common/moduleCaseParam.f90 b/src/modules/common/moduleCaseParam.f90 index 551d867..8df3210 100644 --- a/src/modules/common/moduleCaseParam.f90 +++ b/src/modules/common/moduleCaseParam.f90 @@ -2,13 +2,8 @@ MODULE moduleCaseParam !Final and initial iterations INTEGER:: tFinal, tInitial = 0 - ! Global index of current iteration - INTEGER:: timeStep - ! Time step for all species REAL(8), ALLOCATABLE:: tau(:) - ! Minimum time step REAL(8):: tauMin - ! Time step for Monte-Carlo Collisions REAL(8):: tauColl END MODULE moduleCaseParam diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index 1b8ec67..657ca28 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -40,46 +40,47 @@ MODULE moduleRandom INTEGER:: rnd REAL(8):: rnd01 - rnd = 0 + rnd = 0.D0 CALL RANDOM_NUMBER(rnd01) - rnd = a + FLOOR((b+1-a)*rnd01) + rnd = INT(REAL(b - a) * rnd01) + 1 END FUNCTION randomIntAB - !Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method - function randomMaxwellian() result(rnd) - USE moduleConstParam, only: pi - implicit none - - real(8):: rnd - real(8):: v1, v2, Rsquare - - v1 = 0.d0 - do while (v1 <= 0.d0) - v1 = random() - - end do - v2 = random() - - rnd = sqrt(-2.d0*log(v1))*cos(2*pi*v2) - - end function randomMaxwellian - !Returns a random number in a Maxwellian distribution of mean 0 and width 1 - FUNCTION randomHalfMaxwellian() RESULT(rnd) + FUNCTION randomMaxwellian() RESULT(rnd) + USE moduleConstParam, ONLY: PI IMPLICIT NONE REAL(8):: rnd - REAL(8):: x + REAL(8):: x, y rnd = 0.D0 x = 0.D0 DO WHILE (x == 0.D0) CALL RANDOM_NUMBER(x) END DO + CALL RANDOM_NUMBER(y) - rnd = DSQRT(-DLOG(x)) + rnd = DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + + END FUNCTION randomMaxwellian + + FUNCTION randomHalfMaxwellian() RESULT(rnd) + USE moduleConstParam, ONLY: PI + IMPLICIT NONE + + REAL(8):: rnd + REAL(8):: x, y + + rnd = 0.D0 + x = 0.D0 + DO WHILE (x == 0.D0) + CALL RANDOM_NUMBER(x) + END DO + y = random(-PI, PI) + + rnd = DSQRT(-2.D0*DLOG(x))*DCOS(y) END FUNCTION randomHalfMaxwellian @@ -90,21 +91,10 @@ MODULE moduleRandom REAL(8), INTENT(in):: cumWeight(1:) REAL(8), INTENT(in):: sumWeight REAL(8):: rnd0b - INTEGER:: rnd, i + INTEGER:: rnd - rnd0b = random() - i = 1 - DO - IF (rnd0b <= cumWeight(i)/sumWeight) THEN - rnd = i - EXIT - - ELSE - i = i +1 - - END IF - END DO - ! rnd = MINLOC(DABS(rnd0b - cumWeight), 1) + rnd0b = random(0.D0, sumWeight) + rnd = MINLOC(DABS(rnd0b - cumWeight), 1) END FUNCTION randomWeighted diff --git a/src/modules/common/moduleTable.f90 b/src/modules/common/moduleTable.f90 index 2dc1c75..ac5cbc7 100644 --- a/src/modules/common/moduleTable.f90 +++ b/src/modules/common/moduleTable.f90 @@ -30,6 +30,8 @@ MODULE moduleTable READ(id, '(A)', iostat = stat) dummy !If EOF or error, exit file IF (stat /= 0) EXIT + !If empty line, exit file + IF (dummy == "") EXIT !Skip comment IF (INDEX(dummy,'#') /= 0) CYCLE !Add row @@ -55,6 +57,7 @@ MODULE moduleTable !TODO: Make this a function IF (INDEX(dummy,'#') /= 0) CYCLE IF (stat /= 0) EXIT + IF (dummy == "") EXIT !Add data !TODO: substitute with extracting information from dummy BACKSPACE(id) @@ -93,7 +96,7 @@ MODULE moduleTable f = self%fMax ELSE - i = MINLOC(ABS(x - self%x), 1) + i = MINLOC(x - self%x, 1) deltaX = x - self%x(i) IF (deltaX < 0 ) THEN i = i - 1 diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 9f50947..9f16195 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -84,6 +84,20 @@ MODULE moduleInput CALL readParallel(config) CALL checkStatus(config, "readParallel") + !If everything is correct, creates the output folder + CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder ) + !Copies input file to output folder + CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder) + !Copies particle mesh + IF (mesh%dimen > 0) THEN + CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder) + IF (doubleMesh) THEN + CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder) + + END IF + + END IF + END SUBROUTINE readConfig !Checks the status of the JSON case file and, if failed, exits the execution. @@ -259,17 +273,13 @@ MODULE moduleInput !Read BC CALL readEMBoundary(config) - CASE("ElectrostaticBoltzmann") - !Read BC - CALL readEMBoundary(config) - CASE("ConstantB") !Read BC CALL readEMBoundary(config) !Read constant magnetic field DO i = 1, 3 - WRITE(iString, '(i2)') i - CALL config%get(object // '.B(' // iString // ')', B(i), found) + WRITE(istring, '(i2)') i + CALL config%get(object // '.B(' // istring // ')', B(i), found) IF (.NOT. found) THEN CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver') @@ -312,7 +322,7 @@ MODULE moduleInput LOGICAL:: found CHARACTER(:), ALLOCATABLE:: object INTEGER:: nInitial - INTEGER:: i, p, e + INTEGER:: i, j, p, e CHARACTER(LEN=2):: iString CHARACTER(:), ALLOCATABLE:: spName INTEGER:: sp @@ -328,12 +338,10 @@ MODULE moduleInput REAL(8):: densityCen !Mean velocity and temperature at particle position REAL(8):: velocityXi(1:3), temperatureXi - INTEGER:: nNewPart = 0 - REAL(8):: weight = 0.D0 - CLASS(meshCell), POINTER:: cell + INTEGER:: nNewPart = 0.D0 TYPE(particle), POINTER:: partNew + CLASS(meshCell), POINTER:: cell REAL(8):: vTh - TYPE(lNode), POINTER:: partCurr, partNext CALL config%info('solver.initial', found, n_children = nInitial) @@ -348,9 +356,6 @@ MODULE moduleInput !Reads node values at the nodes filename = path // spFile CALL mesh%readInitial(filename, density, velocity, temperature) - !Check if initial number of particles is given - CALL config%get(object // '.particlesPerCell', nNewPart, found) - !For each volume in the node, create corresponding particles DO e = 1, mesh%numCells !Scale variables @@ -363,15 +368,14 @@ MODULE moduleInput densityCen = mesh%cells(e)%obj%gatherF((/ 0.D0, 0.D0, 0.D0 /), nNodes, sourceScalar) !Calculate number of particles - IF (.NOT. found) THEN - nNewPart = FLOOR(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) + nNewPart = INT(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) - END IF - weight = densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / REAL(nNewPart) + !Allocate array of particles for this species + ALLOCATE(partOld(sp)%p(1:nNewPart)) - !Allocate new particles + !Create new particles DO p = 1, nNewPart - ALLOCATE(partNew) + partNew = partOld(sp)%p(p) partNew%species => species(sp)%obj partNew%r = mesh%cells(e)%obj%randPos() partNew%Xi = mesh%cells(e)%obj%phy2log(partNew%r) @@ -404,10 +408,7 @@ MODULE moduleInput partNew%n_in = .TRUE. - partNew%weight = weight - - !Assign particle to temporal list of particles - CALL partInitial%add(partNew) + partNew%weight = species(sp)%obj%weight !Assign particle to list in volume IF (listInCells) THEN @@ -428,30 +429,10 @@ MODULE moduleInput END DO + nPartOld(sp) = SIZE(partOld(sp)%p) + END DO - !Convert temporal list of particles into initial partOld array - !Deallocate the list of initial particles - nNewPart = partInitial%amount - IF (nNewPart > 0) THEN - ALLOCATE(partOld(1:nNewPart)) - partCurr => partInitial%head - DO p = 1, nNewPart - partNext => partCurr%next - partOld(p) = partCurr%part - DEALLOCATE(partCurr) - partCurr => partNext - - END DO - - IF (ASSOCIATED(partInitial%head)) NULLIFY(partInitial%head) - IF (ASSOCIATED(partInitial%tail)) NULLIFY(partInitial%tail) - partInitial%amount = 0 - - END IF - - nPartOld = SIZE(partOld) - END IF END SUBROUTINE readInitial @@ -598,7 +579,10 @@ MODULE moduleInput END DO + !Allocate the wrapper array that contains particles + ALLOCATE(partOld(1:nSpecies)) !Set number of particles to 0 for init state + ALLOCATE(nPartOld(1:nSpecies)) nPartOld = 0 !Initialize the lock for the non-analogue (NA) list of particles @@ -632,7 +616,7 @@ MODULE moduleInput INTEGER:: i, k, ij INTEGER:: pt_i, pt_j REAL(8):: energyThreshold, energyBinding - CHARACTER(:), ALLOCATABLE:: electron, electronSecondary + CHARACTER(:), ALLOCATABLE:: electron INTEGER:: e CLASS(meshCell), POINTER:: cell @@ -709,16 +693,8 @@ MODULE moduleInput IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions') CALL config%get(object // '.electron', electron, found) IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions') - CALL config%get(object // '.electronSecondary', electronSecondary, found) - IF (found) THEN - CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, & - crossSecFilePath, energyThreshold, electron, electronSecondary) - - ELSE - CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, & - crossSecFilePath, energyThreshold, electron) - - END IF + CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, & + crossSecFilePath, energyThreshold, electron) CASE ('recombination') !Electorn impact ionization @@ -803,7 +779,7 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config INTEGER:: i, s - CHARACTER(2):: iString, sString + CHARACTER(2):: istring, sString CHARACTER(:), ALLOCATABLE:: object, bType REAL(8):: Tw, cw !Wall temperature and specific heat !Neutral Properties @@ -811,16 +787,16 @@ MODULE moduleInput REAL(8), DIMENSION(:), ALLOCATABLE:: v0 REAL(8):: effTime REAL(8):: eThreshold !Energy threshold - INTEGER:: speciesID, electronSecondaryID - CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, electronSecondary + INTEGER:: speciesID + CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, yield LOGICAL:: found INTEGER:: nTypes CALL config%info('boundary', found, n_children = nBoundary) ALLOCATE(boundary(1:nBoundary)) DO i = 1, nBoundary - WRITE(iString, '(i2)') i - object = 'boundary(' // TRIM(iString) // ')' + WRITE(istring, '(i2)') i + object = 'boundary(' // TRIM(istring) // ')' boundary(i)%n = i CALL config%get(object // '.name', boundary(i)%name, found) @@ -829,77 +805,71 @@ MODULE moduleInput IF (nTypes /= nSpecies) CALL criticalError('Not enough boundary types defined in ' // object, 'readBoundary') ALLOCATE(boundary(i)%bTypes(1:nSpecies)) DO s = 1, nSpecies - associate(bound => boundary(i)%bTypes(s)%obj) - WRITE(sString,'(i2)') s - object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')' - CALL config%get(object // '.type', bType, found) - SELECT CASE(bType) - CASE('reflection') - ALLOCATE(boundaryReflection:: bound) + WRITE(sString,'(i2)') s + object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')' + CALL config%get(object // '.type', bType, found) + SELECT CASE(bType) + CASE('reflection') + ALLOCATE(boundaryReflection:: boundary(i)%bTypes(s)%obj) - CASE('absorption') - ALLOCATE(boundaryAbsorption:: bound) + CASE('absorption') + ALLOCATE(boundaryAbsorption:: boundary(i)%bTypes(s)%obj) - CASE('transparent') - ALLOCATE(boundaryTransparent:: bound) + CASE('transparent') + ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj) - CASE('axis') - ALLOCATE(boundaryAxis:: bound) + CASE('ionization') + !Neutral parameters + CALL config%get(object // '.neutral.ion', speciesName, found) + IF (.NOT. found) CALL criticalError("missing parameter 'ion' for neutrals in ionization", 'readBoundary') + speciesID = speciesName2Index(speciesName) + CALL config%get(object // '.neutral.mass', m0, found) + IF (.NOT. found) THEN + m0 = species(s)%obj%m*m_ref + END IF + CALL config%get(object // '.neutral.density', n0, found) + IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary') + CALL config%get(object // '.neutral.velocity', v0, found) + IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary') + CALL config%get(object // '.neutral.temperature', T0, found) + IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary') - CASE('wallTemperature') - CALL config%get(object // '.temperature', Tw, found) - IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary') - CALL config%get(object // '.specificHeat', cw, found) - IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary') + CALL config%get(object // '.effectiveTime', effTime, found) + IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary') - CALL initWallTemperature(bound, Tw, cw) + CALL config%get(object // '.energyThreshold', eThreshold, found) + IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary') - CASE('ionization') - !Neutral parameters - CALL config%get(object // '.neutral.ion', speciesName, found) - IF (.NOT. found) CALL criticalError("missing parameter 'ion' for neutrals in ionization", 'readBoundary') - speciesID = speciesName2Index(speciesName) - CALL config%get(object // '.neutral.mass', m0, found) - IF (.NOT. found) THEN - m0 = species(s)%obj%m*m_ref - END IF - CALL config%get(object // '.neutral.density', n0, found) - IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary') - CALL config%get(object // '.neutral.velocity', v0, found) - IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary') - CALL config%get(object // '.neutral.temperature', T0, found) - IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary') + CALL config%get(object // '.crossSection', crossSection, found) + IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary') - CALL config%get(object // '.effectiveTime', effTime, found) - IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary') + CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, & + speciesID, effTime, crossSection, eThreshold) - CALL config%get(object // '.energyThreshold', eThreshold, found) - IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary') + CASE('wallTemperature') + CALL config%get(object // '.temperature', Tw, found) + IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary') + CALL config%get(object // '.specificHeat', cw, found) + IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary') - CALL config%get(object // '.crossSection', crossSection, found) - IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary') + CALL initWallTemperature(boundary(i)%bTypes(s)%obj, Tw, cw) - CALL config%get(object // '.electronSecondary', electronSecondary, found) - electronSecondaryID = speciesName2Index(electronSecondary) - IF (found) THEN - CALL initIonization(bound, species(s)%obj%m, m0, n0, v0, T0, & - speciesID, effTime, crossSection, eThreshold,electronSecondaryID) + CASE('secondaryEmission') + CALL config%get(object // '.yield', yield, found) + IF (.NOT. found) CALL criticalError("missing parameter 'yield' for secondary emission", 'readBoundary') + CALL config%get(object // '.electron', speciesName, found) + IF (.NOT. found) CALL criticalError("missing parameter 'electron' for secondary emission", 'readBoundary') + speciesID = speciesName2Index(speciesName) - ELSE - CALL initIonization(bound, species(s)%obj%m, m0, n0, v0, T0, & - speciesID, effTime, crossSection, eThreshold) + CALL initSEE(boundary(i)%bTypes(s)%obj, yield, speciesID) - END IF + CASE('axis') + ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj) - case('outflowAdaptive') - allocate(boundaryOutflowAdaptive:: bound) + CASE DEFAULT + CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary') - CASE DEFAULT - CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary') - - END SELECT - - end associate + END SELECT END DO @@ -916,7 +886,6 @@ MODULE moduleInput USE moduleMeshInputGmsh2, ONLY: initGmsh2 USE moduleMeshInputVTU, ONLY: initVTU USE moduleMeshInput0D, ONLY: init0D - USE moduleMeshInputText, ONLY: initText USE moduleMesh3DCart USE moduleMesh2DCyl USE moduleMesh2DCart @@ -934,6 +903,7 @@ MODULE moduleInput LOGICAL:: found CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile REAL(8):: volume + CHARACTER(:), ALLOCATABLE:: meshFileVTU !Temporary to test VTU OUTPUT object = 'geometry' @@ -975,9 +945,9 @@ MODULE moduleInput !Read the 0D mesh CALL mesh%readMesh(pathMeshParticle) - !Get the volume + !Get the volumne CALL config%get(object // '.volume', volume, found) - !Rescale the volume + !Rescale the volumne IF (found) THEN mesh%cells(1)%obj%volume = mesh%cells(1)%obj%volume*volume / Vol_ref mesh%nodes(1)%obj%v = mesh%cells(1)%obj%volume @@ -1065,20 +1035,6 @@ MODULE moduleInput END IF - case ("text") - !Check if the geometry is right. - if (mesh%dimen /= 1) then - call criticalError("Text mesh is only allowed for 1D geometries", 'readGeometry') - - end if - - !Read the mesh - call initText(mesh) - if (doubleMesh) then - call initText(meshColl) - - end if - CASE DEFAULT CALL criticalError('Mesh format ' // meshFormat // ' not defined.', 'readGeometry') @@ -1125,13 +1081,13 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config CHARACTER(:), ALLOCATABLE:: object LOGICAL:: found - CHARACTER(2):: iString + CHARACTER(2):: istring INTEGER:: i CHARACTER(:), ALLOCATABLE:: speciesName REAL(8), ALLOCATABLE, DIMENSION(:):: r REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3 INTEGER, ALLOCATABLE, DIMENSION(:):: points - REAL(8):: everyTimeStep + REAL(8):: timeStep CALL config%info('output.probes', found, n_children = nProbes) @@ -1139,7 +1095,7 @@ MODULE moduleInput DO i = 1, nProbes WRITE(iString, '(I2)') i - object = 'output.probes(' // trim(iString) // ')' + object = 'output.probes(' // trim(istring) // ')' CALL config%get(object // '.species', speciesName, found) CALL config%get(object // '.position', r, found) @@ -1147,14 +1103,16 @@ MODULE moduleInput CALL config%get(object // '.velocity_2', v2, found) CALL config%get(object // '.velocity_3', v3, found) CALL config%get(object // '.points', points, found) - CALL config%get(object // '.timeStep', everyTimeStep, found) + CALL config%get(object // '.timeStep', timeStep, found) IF (ANY(points < 2)) CALL criticalError("Number of points in probe " // iString // " incorrect", 'readProbes') - CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, everyTimeStep) + CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, timeStep) END DO + CALL resetProbes(tInitial) + END SUBROUTINE readProbes SUBROUTINE readEMBoundary(config) @@ -1162,6 +1120,7 @@ MODULE moduleInput USE moduleOutput USE moduleErrors USE moduleEM + USE moduleRefParam USE moduleSpecies USE json_module IMPLICIT NONE @@ -1169,72 +1128,34 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config CHARACTER(:), ALLOCATABLE:: object LOGICAL:: found - CHARACTER(:), ALLOCATABLE:: typeEM - REAL(8):: potential - INTEGER:: physicalSurface - CHARACTER(:), ALLOCATABLE:: temporalProfile, temporalProfilePath - INTEGER:: b, s, n, ni - CHARACTER(2):: bString + CHARACTER(2):: istring + INTEGER:: i, e, s INTEGER:: info EXTERNAL:: dgetrf CALL config%info('boundaryEM', found, n_children = nBoundaryEM) - IF (found) THEN - ALLOCATE(boundaryEM(1:nBoundaryEM)) - - END IF + IF (found) ALLOCATE(boundEM(1:nBoundaryEM)) - DO b = 1, nBoundaryEM - WRITE(bString, '(I2)') b - object = 'boundaryEM(' // TRIM(bString) // ')' + DO i = 1, nBoundaryEM + WRITE(istring, '(I2)') i + object = 'boundaryEM(' // trim(istring) // ')' - CALL config%get(object // '.type', typeEM, found) + CALL config%get(object // '.type', boundEM(i)%typeEM, found) - SELECT CASE(typeEM) + SELECT CASE(boundEM(i)%typeEM) CASE ("dirichlet") - CALL config%get(object // '.potential', potential, found) - IF (.NOT. found) THEN + CALL config%get(object // '.potential', boundEM(i)%potential, found) + IF (.NOT. found) & CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary') + boundEM(i)%potential = boundEM(i)%potential/Volt_ref - END IF - - CALL config%get(object // '.physicalSurface', physicalSurface, found) - IF (.NOT. found) THEN - CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', & - 'readEMBoundary') - - END IF - - CALL initDirichlet(boundaryEM(b)%obj, physicalSurface, potential) - - CASE ("dirichletTime") - CALL config%get(object // '.potential', potential, found) - IF (.NOT. found) THEN - CALL criticalError('Required parameter "potential" for Dirichlet Time boundary condition not found', & - 'readEMBoundary') - - END IF - - CALL config%get(object // '.temporalProfile', temporalProfile, found) - IF (.NOT. found) THEN - CALL criticalError('Required parameter "temporalProfile" for Dirichlet Time boundary condition not found', & - 'readEMBoundary') - - END IF - temporalProfilePath = path // temporalProfile - - CALL config%get(object // '.physicalSurface', physicalSurface, found) - IF (.NOT. found) THEN - CALL criticalError('Required parameter "physicalSurface" for Dirichlet Time boundary condition not found', & - 'readEMBoundary') - - END IF - - CALL initDirichletTime(boundaryEM(b)%obj, physicalSurface, potential, temporalProfilePath) + CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found) + IF (.NOT. found) & + CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary') CASE DEFAULT - CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary') + CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary') END SELECT @@ -1253,28 +1174,18 @@ MODULE moduleInput END DO - ! Modify K matrix due to boundary conditions - DO b = 1, nBoundaryEM - SELECT TYPE(boundary => boundaryEM(b)%obj) - TYPE IS(boundaryEMDirichlet) - DO n = 1, boundary%nNodes - ni = boundary%nodes(n)%obj%n - mesh%K(ni, :) = 0.D0 - mesh%K(ni, ni) = 1.D0 + IF (ALLOCATED(boundEM)) THEN + DO e = 1, mesh%numEdges + IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN + DO i = 1, nBoundaryEM + IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN + CALL boundEM(i)%apply(mesh%edges(e)%obj) - END DO - - TYPE IS(boundaryEMDirichletTime) - DO n = 1, boundary%nNodes - ni = boundary%nodes(n)%obj%n - mesh%K(ni, :) = 0.D0 - mesh%K(ni, ni) = 1.D0 - - END DO - - END SELECT - - END DO + END IF + END DO + END IF + END DO + END IF !Compute the PLU factorization of K once boundary conditions have been read CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info) @@ -1295,25 +1206,24 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config INTEGER:: i - CHARACTER(2):: iString + CHARACTER(2):: istring CHARACTER(:), ALLOCATABLE:: object LOGICAL:: found CHARACTER(:), ALLOCATABLE:: speciesName CHARACTER(:), ALLOCATABLE:: name REAL(8):: v - REAL(8), ALLOCATABLE:: temperature(:), normal(:) + REAL(8), ALLOCATABLE:: T(:), normal(:) REAL(8):: flow CHARACTER(:), ALLOCATABLE:: units INTEGER:: physicalSurface - INTEGER:: particlesPerEdge INTEGER:: sp CALL config%info('inject', found, n_children = nInject) ALLOCATE(inject(1:nInject)) nPartInj = 0 DO i = 1, nInject - WRITE(iString, '(i2)') i - object = 'inject(' // trim(iString) // ')' + WRITE(istring, '(i2)') i + object = 'inject(' // trim(istring) // ')' !Find species CALL config%get(object // '.species', speciesName, found) @@ -1321,7 +1231,7 @@ MODULE moduleInput CALL config%get(object // '.name', name, found) CALL config%get(object // '.v', v, found) - CALL config%get(object // '.T', temperature, found) + CALL config%get(object // '.T', T, found) CALL config%get(object // '.n', normal, found) IF (.NOT. found) THEN ALLOCATE(normal(1:3)) @@ -1330,10 +1240,8 @@ MODULE moduleInput CALL config%get(object // '.flow', flow, found) CALL config%get(object // '.units', units, found) CALL config%get(object // '.physicalSurface', physicalSurface, found) - particlesPerEdge = 0 - CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found) - CALL inject(i)%init(i, v, normal, temperature, flow, units, sp, physicalSurface, particlesPerEdge) + CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface) CALL readVelDistr(config, inject(i), object) @@ -1352,7 +1260,6 @@ MODULE moduleInput USE moduleCaseParam, ONLY: tauMin USE moduleMesh, ONLY: mesh USE moduleSpecies, ONLY: nSpecies - USE moduleRefParam, ONLY: ti_ref IMPLICIT NONE TYPE(json_file), INTENT(inout):: config @@ -1366,11 +1273,8 @@ MODULE moduleInput CALL config%get('average.startTime', tStart, found) IF (found) THEN - tAverageStart = INT(tStart / ti_ref / tauMin) + tAverageStart = INT(tStart / tauMin) - ELSE - tAverageStart = 0 - END IF ALLOCATE(averageScheme(1:mesh%numNodes)) @@ -1397,28 +1301,28 @@ MODULE moduleInput TYPE(injectGeneric), INTENT(inout):: inj CHARACTER(:), ALLOCATABLE, INTENT(in):: object INTEGER:: i - CHARACTER(2):: iString + CHARACTER(2):: istring CHARACTER(:), ALLOCATABLE:: fvType LOGICAL:: found - REAL(8):: v, temperature, m + REAL(8):: v, T, m !Reads species mass m = inj%species%m !Reads distribution functions for velocity DO i = 1, 3 - WRITE(iString, '(i2)') i - CALL config%get(object // '.velDist('// TRIM(iString) //')', fvType, found) - IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // iString // & + WRITE(istring, '(i2)') i + CALL config%get(object // '.velDist('// TRIM(istring) //')', fvType, found) + IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // istring // & " found for " // object, 'readVelDistr') SELECT CASE(fvType) CASE ("Maxwellian") - temperature = inj%temperature(i) - CALL initVelDistMaxwellian(inj%v(i)%obj, temperature, m) + T = inj%T(i) + CALL initVelDistMaxwellian(inj%v(i)%obj, t, m) CASE ("Half-Maxwellian") - temperature = inj%temperature(i) - CALL initVelDistHalfMaxwellian(inj%v(i)%obj, temperature, m) + T = inj%T(i) + CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m) CASE ("Delta") v = inj%vMod*inj%n(i) @@ -1459,37 +1363,5 @@ MODULE moduleInput END SUBROUTINE readParallel - SUBROUTINE initOutput(inputFile) - USE moduleRefParam - USE moduleMesh, ONLY: mesh, doubleMesh, pathMeshParticle, pathMeshColl - USE moduleOutput, ONLY: path, folder - IMPLICIT NONE - - CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile - INTEGER:: fileReference = 30 - !If everything is correct, creates the output folder - CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder ) - !Copies input file to output folder - CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder) - !Copies particle mesh - IF (mesh%dimen > 0) THEN - CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder) - IF (doubleMesh) THEN - CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder) - - END IF - - END IF - - ! Write commit of fpakc - CALL SYSTEM('git rev-parse HEAD > ' // path // folder // '/' // 'fpakc_commit.txt') - - ! Write file with reference values - OPEN (fileReference, file=path // folder // '/' // 'reference.txt') - WRITE(fileReference, "(7(1X,A20))") 'L_ref', 'v_ref', 'ti_ref', 'Vol_ref', 'EF_ref', 'Volt_ref', 'B_ref' - WRITE(fileReference, "(7(1X,ES20.6E3))") L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref - CLOSE(fileReference) - - END SUBROUTINE initOutput END MODULE moduleInput diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index a304c5e..269f157 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -122,8 +122,6 @@ MODULE moduleMesh1DCart self%x = r1(1) - self%surface = 1.D0 - self%normal = (/ 1.D0, 0.D0, 0.D0 /) !Boundary index diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index a46c6a2..d998267 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -122,8 +122,6 @@ MODULE moduleMesh1DRad self%r = r1(1) - self%surface = 1.D0 - self%normal = (/ 1.D0, 0.D0, 0.D0 /) !Boundary index diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index db1ef4f..4341cb0 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -163,7 +163,7 @@ MODULE moduleMesh2DCart r2 = self%n2%getCoordinates() self%x = (/r1(1), r2(1)/) self%y = (/r1(2), r2(2)/) - self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2) + self%weight = 1.D0 !Normal vector self%normal = (/ -(self%y(2)-self%y(1)), & self%x(2)-self%x(1) , & @@ -318,8 +318,6 @@ MODULE moduleMesh2DCart INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) - fPsi = 0.D0 - fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), & (1.D0 + Xi(1)) * (1.D0 - Xi(2)), & (1.D0 + Xi(1)) * (1.D0 + Xi(2)), & @@ -494,36 +492,34 @@ MODULE moduleMesh2DCart END FUNCTION insideQuad - !Transform physical coordinates to element coordinates with a Taylor series + !Transform physical coordinates to element coordinates PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2) + REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) REAL(8):: dPsi(1:3,1:4), fPsi(1:4) REAL(8):: pDer(1:3, 1:3) REAL(8):: conv !Iterative newton method to transform coordinates - conv = 1.D0 - Xi0 = 0.D0 - Xi(3) = 0.D0 + conv = 1.D0 + XiO = 0.D0 DO WHILE(conv > 1.D-4) - fPsi = self%fPsi(Xi0, 4) - x0(1) = dot_product(fPsi, self%x) - x0(2) = dot_product(fPsi, self%y) - deltaR = r(1:2) - x0 - dPsi = self%dPsi(Xi0, 4) - pDer = self%partialDer(4, dPsi) - detJ = self%detJac(pDer) - pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) - pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) - Xi(1:2) = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ - conv = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1) - Xi0(1:2) = Xi(1:2) + dPsi = self%dPsi(XiO, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + fPsi = self%fPsi(XiO, 4) + f = (/ DOT_PRODUCT(fPsi,self%x), & + DOT_PRODUCT(fPsi,self%y), & + 0.D0 /) - r + Xi = XiO - MATMUL(invJ, f)/detJ + conv = MAXVAL(DABS(Xi-XiO),1) + XiO = Xi END DO @@ -573,7 +569,6 @@ MODULE moduleMesh2DCart pDer = self%partialDer(4, dPsi) detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 4) - !Compute total volume of the cell self%volume = detJ*4.D0 !Compute volume per node @@ -680,8 +675,8 @@ MODULE moduleMesh2DCart dPsi = 0.D0 - dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /) - dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /) + dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /) + dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTria @@ -767,7 +762,6 @@ MODULE moduleMesh2DCart pDer = self%partialDer(3, dPsi) detJ = self%detJac(pDer) invJ = self%invJac(pDer) - localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ END DO @@ -826,19 +820,19 @@ MODULE moduleMesh2DCart CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2) - REAL(8):: dPsi(1:3,1:4) + REAL(8):: deltaR(1:3) + REAL(8):: dPsi(1:3, 1:3) REAL(8):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3, 1:3), detJ !Direct method to convert coordinates - Xi(3) = 0.D0 - deltaR = (/ r(1) - self%x(1), r(2) - self%y(1) /) - dPsi = self%dPsi(Xi, 3) - pDer = self%partialDer(3, dPsi) - detJ = self%detJac(pDer) - pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) - pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) - Xi(1:2) = MATMUL(pDerInv,deltaR)/detJ + Xi = 0.D0 + deltaR = (/ r(1) - self%x(1), r(2) - self%y(1), 0.D0 /) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + invJ = self%invJac(pDer) + detJ = self%detJac(pDer) + Xi = MATMUL(invJ,deltaR)/detJ END FUNCTION phy2logTria @@ -913,8 +907,8 @@ MODULE moduleMesh2DCart invJ = 0.D0 - invJ(1, 1:2) = (/ pDer(2,2), -pDer(2,1) /) - invJ(2, 1:2) = (/ -pDer(1,2), pDer(1,1) /) + invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) + invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) invJ(3, 3) = 1.D0 END FUNCTION invJ2DCart diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 1f5f33c..d4baedd 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -144,7 +144,6 @@ MODULE moduleMesh2DCyl USE moduleSpecies USE moduleBoundary USE moduleErrors - USE moduleConstParam, ONLY: PI IMPLICIT NONE CLASS(meshEdge2DCyl), INTENT(out):: self @@ -164,15 +163,7 @@ MODULE moduleMesh2DCyl r2 = self%n2%getCoordinates() self%z = (/r1(1), r2(1)/) self%r = (/r1(2), r2(2)/) - !Edge surface - IF (self%z(2) /= self%z(1)) THEN - self%surface = ABS(self%r(2) + self%r(1))*ABS(self%z(2) - self%z(1)) - - ELSE - self%surface = ABS(self%r(2)**2 - self%r(1)**2) - - END IF - self%surface = self%surface * PI + self%weight = r2(2)**2 - r1(2)**2 !Normal vector self%normal = (/ -(self%r(2)-self%r(1)), & self%z(2)-self%z(1) , & @@ -232,13 +223,21 @@ MODULE moduleMesh2DCyl CLASS(meshEdge2DCyl), INTENT(in):: self REAL(8):: rnd REAL(8):: r(1:3) - REAL(8):: p1(1:2), p2(1:2) + REAL(8):: dr, dz rnd = random() + dr = self%r(2) - self%r(1) + dz = self%z(2) - self%z(1) + IF (dr /= 0.D0) THEN + r(2) = dr * DSQRT(rnd) + self%r(1) + r(1) = dz * (r(2) - self%r(1))/dr + self%z(1) + + ELSE + r(2) = self%r(1) + r(1) = dz * rnd + self%z(1) + + END IF - p1 = (/self%z(1), self%r(1) /) - p2 = (/self%z(2), self%r(2) /) - r(1:2) = (1.D0 - rnd)*p1 + rnd*p2 r(3) = 0.D0 END FUNCTION randPosEdge @@ -247,6 +246,7 @@ MODULE moduleMesh2DCyl !QUAD FUNCTIONS !Init element SUBROUTINE initCellQuad2DCyl(self, n, p, nodes) + USE moduleRefParam IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(out):: self @@ -326,8 +326,6 @@ MODULE moduleMesh2DCyl INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) - fPsi = 0.D0 - fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), & (1.D0 + Xi(1)) * (1.D0 - Xi(2)), & (1.D0 + Xi(1)) * (1.D0 + Xi(2)), & @@ -498,7 +496,7 @@ MODULE moduleMesh2DCyl END FUNCTION elemFQuad - !Check if Xi is inside the element + !Checks if Xi is inside the element PURE FUNCTION insideQuad(Xi) RESULT(ins) IMPLICIT NONE @@ -510,36 +508,34 @@ MODULE moduleMesh2DCyl END FUNCTION insideQuad - !Transform physical coordinates to element coordinates with a Taylor series + !Transform physical coordinates to element coordinates PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2) + REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) REAL(8):: dPsi(1:3,1:4), fPsi(1:4) REAL(8):: pDer(1:3, 1:3) REAL(8):: conv !Iterative newton method to transform coordinates - conv = 1.D0 - Xi0 = 0.D0 - Xi(3) = 0.D0 + conv = 1.D0 + XiO = 0.D0 DO WHILE(conv > 1.D-4) - fPsi = self%fPsi(Xi0, 4) - x0(1) = dot_product(fPsi, self%z) - x0(2) = dot_product(fPsi, self%r) - deltaR = r(1:2) - x0 - dPsi = self%dPsi(Xi0, 4) - pDer = self%partialDer(4, dPsi) - detJ = self%detJac(pDer) - pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) - pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) - Xi(1:2) = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ - conv = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1) - Xi0(1:2) = Xi(1:2) + dPsi = self%dPsi(XiO, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + fPsi = self%fPsi(XiO, 4) + f = (/ DOT_PRODUCT(fPsi,self%z), & + DOT_PRODUCT(fPsi,self%r), & + 0.D0 /) - r + Xi = XiO - MATMUL(invJ, f)/detJ + conv = MAXVAL(DABS(Xi-XiO),1) + XiO = Xi END DO @@ -557,7 +553,7 @@ MODULE moduleMesh2DCyl XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /) nextInt = MAXLOC(XiArray,1) - !Select the higher value of directions and searches in that direction + !Selects the higher value of directions and searches in that direction NULLIFY(neighbourElement) SELECT CASE (nextInt) CASE (1) @@ -585,7 +581,6 @@ MODULE moduleMesh2DCyl REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3) self%volume = 0.D0 - !2D 1 point Gauss Quad Integral Xi = 0.D0 dPsi = self%dPsi(Xi, 4) @@ -594,18 +589,18 @@ MODULE moduleMesh2DCyl fPsi = self%fPsi(Xi, 4) r = DOT_PRODUCT(fPsi,self%r) !Computes total volume of the cell - self%volume = r*detJ*PI8 !2*pi * 4 (weight of 1 point 2D-Gaussian integral) - !Computes volume per node. Change the radius point to calculate the area to improve accuracy near the axis. - Xi = (/-5.D-1, -0.25D0, 0.D0/) + self%volume = r*detJ*PI8 !4*2*pi + !Computes volume per node + Xi = (/-5.D-1, -5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8 - Xi = (/ 5.D-1, -0.25D0, 0.D0/) + Xi = (/ 5.D-1, -5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8 - Xi = (/ 5.D-1, 0.75D0, 0.D0/) + Xi = (/ 5.D-1, 5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8 - Xi = (/-5.D-1, 0.75D0, 0.D0/) + Xi = (/-5.D-1, 5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8 @@ -707,8 +702,8 @@ MODULE moduleMesh2DCyl dPsi = 0.D0 - dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /) - dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /) + dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /) + dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTria @@ -860,19 +855,19 @@ MODULE moduleMesh2DCyl CLASS(meshCell2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2) - REAL(8):: dPsi(1:3,1:4) + REAL(8):: deltaR(1:3) + REAL(8):: dPsi(1:3, 1:3) REAL(8):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3, 1:3), detJ !Direct method to convert coordinates - Xi(3) = 0.D0 - deltaR = (/ r(1) - self%z(1), r(2) - self%r(1) /) - dPsi = self%dPsi(Xi, 3) - pDer = self%partialDer(3, dPsi) - detJ = self%detJac(pDer) - pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) - pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) - Xi(1:2) = MATMUL(pDerInv,deltaR)/detJ + Xi = 0.D0 + deltaR = (/ r(1) - self%z(1), r(2) - self%r(1), 0.D0 /) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + invJ = self%invJac(pDer) + detJ = self%detJac(pDer) + Xi = MATMUL(invJ,deltaR)/detJ END FUNCTION phy2logTria @@ -950,8 +945,8 @@ MODULE moduleMesh2DCyl invJ = 0.D0 - invJ(1, 1:2) = (/ pDer(2,2), -pDer(2,1) /) - invJ(2, 1:2) = (/ -pDer(1,2), pDer(1,1) /) + invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) + invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) invJ(3, 3) = 1.D0 END FUNCTION invJ2DCyl diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 8170df7..c451689 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -109,7 +109,6 @@ MODULE moduleMesh3DCart USE moduleBoundary USE moduleErrors USE moduleMath - USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge3DCartTria), INTENT(out):: self @@ -143,8 +142,6 @@ MODULE moduleMesh3DCart self%normal = crossProduct(vec1, vec2) self%normal = normalize(self%normal) - self%surface = 1.D0/L_ref**2 !TODO: FIX THIS WHEN MOVING TO 3D - !Boundary index self%boundary => boundary(bt) ALLOCATE(self%fBoundary(1:nSpecies)) diff --git a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 index c0dcfbb..97ec729 100644 --- a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 +++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 @@ -1,22 +1,22 @@ MODULE moduleMeshOutput0D CONTAINS - SUBROUTINE printOutput0D(self) + SUBROUTINE printOutput0D(self, t) USE moduleMesh USE moduleRefParam USE moduleSpecies USE moduleOutput - USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: t INTEGER:: i TYPE(outputFormat):: output CHARACTER(:), ALLOCATABLE:: fileName DO i = 1, nSpecies fileName='OUTPUT_' // species(i)%obj%name // '.dat' - IF (timeStep == 0) THEN + IF (t == 0) THEN OPEN(20, file = path // folder // '/' // fileName, action = 'write') WRITE(20, "(A1, 14X, A5, A20, 40X, A20, 2(A20))") "#","t (s)","density (m^-3)", "velocity (m/s)", & "pressure (Pa)", "temperature (K)" @@ -27,17 +27,14 @@ MODULE moduleMeshOutput0D OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') CALL calculateOutput(self%nodes(1)%obj%output(i), output, self%nodes(1)%obj%v, species(i)%obj) - WRITE(20, "(7(ES20.6E3))") REAL(timeStep)*tauMin*ti_ref, output%density, & - output%velocity, & - output%pressure, & - output%temperature + WRITE(20, "(7(ES20.6E3))") REAL(t)*tauMin*ti_ref, output%density, output%velocity, output%pressure, output%temperature CLOSE(20) END DO END SUBROUTINE printOutput0D - SUBROUTINE printColl0D(self) + SUBROUTINE printColl0D(self, t) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -46,11 +43,12 @@ MODULE moduleMeshOutput0D IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self + INTEGER, INTENT(in):: t CHARACTER(:), ALLOCATABLE:: fileName INTEGER:: k fileName='OUTPUT_Collisions.dat' - IF (timeStep == tInitial) THEN + IF (t == tInitial) THEN OPEN(20, file = path // folder // '/' // fileName, action = 'write') WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions" WRITE(*, "(6X,A15,A)") "Creating file: ", fileName @@ -59,12 +57,12 @@ MODULE moduleMeshOutput0D END IF OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - WRITE(20, "(ES20.6E3, 10I20)") REAL(timeStep)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs) + WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs) CLOSE(20) END SUBROUTINE printColl0D - SUBROUTINE printEM0D(self) + SUBROUTINE printEM0D(self, t) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -72,6 +70,7 @@ MODULE moduleMeshOutput0D IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: t END SUBROUTINE printEM0D diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index 1ad6a36..0df1289 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -108,7 +108,6 @@ MODULE moduleMeshInputGmsh2 READ(10, *) totalNumElem !Count edges and volume elements - numEdges = 0 SELECT TYPE(self) TYPE IS(meshParticles) self%numEdges = 0 @@ -329,7 +328,7 @@ MODULE moduleMeshInputGmsh2 DO i = 1, numNodes !Reads the density - READ(10, *) e, density(i) + READ(10, *), e, density(i) END DO @@ -340,7 +339,7 @@ MODULE moduleMeshInputGmsh2 DO i = 1, numNodes !Reads the velocity - READ(10, *) e, velocity(i, 1:3) + READ(10, *), e, velocity(i, 1:3) END DO diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index 8176bb5..ccdcf3d 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -80,50 +80,50 @@ MODULE moduleMeshOutputGmsh2 END SUBROUTINE writeGmsh2FooterElementData !Prints the scattered properties of particles into the nodes - SUBROUTINE printOutputGmsh2(self) + SUBROUTINE printOutputGmsh2(self, t) USE moduleMesh USE moduleRefParam USE moduleSpecies USE moduleOutput USE moduleMeshInoutCommon - USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: t INTEGER:: n, i TYPE(outputFormat):: output(1:self%numNodes) REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName - time = DBLE(timeStep)*tauMin*ti_ref + time = DBLE(t)*tauMin*ti_ref DO i = 1, nSpecies - fileName = formatFileName(prefix, species(i)%obj%name, 'msh', timeStep) + fileName = formatFileName(prefix, species(i)%obj%name, 'msh', t) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (60, file = path // folder // '/' // fileName) CALL writeGmsh2HeaderMesh(60) - CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', timeStep, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', t, time, 1, self%numNodes) DO n=1, self%numNodes CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj) WRITE(60, "(I6,ES20.6E3)") n, output(n)%density END DO CALL writeGmsh2FooterNodeData(60) - CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', timeStep, time, 3, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', t, time, 3, self%numNodes) DO n=1, self%numNodes WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity END DO CALL writeGmsh2FooterNodeData(60) - CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', timeStep, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', t, time, 1, self%numNodes) DO n=1, self%numNodes WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure END DO CALL writeGmsh2FooterNodeData(60) - CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', timeStep, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', t, time, 1, self%numNodes) DO n=1, self%numNodes WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature END DO @@ -135,7 +135,7 @@ MODULE moduleMeshOutputGmsh2 END SUBROUTINE printOutputGmsh2 !Prints the number of collisions into the volumes - SUBROUTINE printCollGmsh2(self) + SUBROUTINE printCollGmsh2(self, t) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -145,6 +145,7 @@ MODULE moduleMeshOutputGmsh2 IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self + INTEGER, INTENT(in):: t INTEGER:: numEdges INTEGER:: k, c INTEGER:: n @@ -166,9 +167,9 @@ MODULE moduleMeshOutputGmsh2 END SELECT IF (collOutput) THEN - time = DBLE(timeStep)*tauMin*ti_ref + time = DBLE(t)*tauMin*ti_ref - fileName = formatFileName(prefix, 'Collisions', 'msh', timeStep) + fileName = formatFileName(prefix, 'Collisions', 'msh', t) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (60, file = path // folder // '/' // fileName) @@ -178,7 +179,7 @@ MODULE moduleMeshOutputGmsh2 DO c = 1, interactionMatrix(k)%amount WRITE(cString, "(I2)") c title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString - CALL writeGmsh2HeaderElementData(60, title, timeStep, time, 1, self%numCells) + CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numCells) DO n=1, self%numCells WRITE(60, "(I6,I10)") n + numEdges, self%cells(n)%obj%tallyColl(k)%tally(c) END DO @@ -195,7 +196,7 @@ MODULE moduleMeshOutputGmsh2 END SUBROUTINE printCollGmsh2 !Prints the electrostatic EM properties into the nodes and volumes - SUBROUTINE printEMGmsh2(self) + SUBROUTINE printEMGmsh2(self, t) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -204,6 +205,7 @@ MODULE moduleMeshOutputGmsh2 IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: t INTEGER:: n, e REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName @@ -212,27 +214,27 @@ MODULE moduleMeshOutputGmsh2 Xi = (/ 0.D0, 0.D0, 0.D0 /) IF (emOutput) THEN - time = DBLE(timeStep)*tauMin*ti_ref + time = DBLE(t)*tauMin*ti_ref - fileName = formatFileName(prefix, 'EMField', 'msh', timeStep) + fileName = formatFileName(prefix, 'EMField', 'msh', t) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (20, file = path // folder // '/' // fileName) CALL writeGmsh2HeaderMesh(20) - CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', timeStep, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', t, time, 1, self%numNodes) DO n=1, self%numNodes WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref END DO CALL writeGmsh2FooterNodeData(20) - CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', timeStep, time, 3, self%numCells) + CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numCells) DO e=1, self%numCells WRITE(20, *) e+self%numEdges, self%cells(e)%obj%gatherElectricField(Xi)*EF_ref END DO CALL writeGmsh2FooterElementData(20) - CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', timeStep, time, 3, self%numNodes) + CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', t, time, 3, self%numNodes) DO n=1, self%numNodes WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref END DO diff --git a/src/modules/mesh/inout/makefile b/src/modules/mesh/inout/makefile index 1b8883d..a93161d 100644 --- a/src/modules/mesh/inout/makefile +++ b/src/modules/mesh/inout/makefile @@ -1,4 +1,4 @@ -all: vtu.o gmsh2.o 0D.o text.o +all: vtu.o gmsh2.o 0D.o vtu.o: moduleMeshInoutCommon.o $(MAKE) -C vtu all @@ -9,8 +9,5 @@ gmsh2.o: 0D.o: $(MAKE) -C 0D all -text.o: - $(MAKE) -C text all - %.o: %.f90 $(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/inout/moduleMeshInoutCommon.f90 b/src/modules/mesh/inout/moduleMeshInoutCommon.f90 index 7dc2e84..e4a6c72 100644 --- a/src/modules/mesh/inout/moduleMeshInoutCommon.f90 +++ b/src/modules/mesh/inout/moduleMeshInoutCommon.f90 @@ -3,17 +3,17 @@ MODULE moduleMeshInoutCommon CHARACTER(LEN=4):: prefix = 'Step' CONTAINS - PURE FUNCTION formatFileName(prefix, suffix, extension, timeStep) RESULT(fileName) + PURE FUNCTION formatFileName(prefix, suffix, extension, t) RESULT(fileName) USE moduleOutput IMPLICIT NONE CHARACTER(*), INTENT(in):: prefix, suffix, extension - INTEGER, INTENT(in), OPTIONAL:: timeStep + INTEGER, INTENT(in), OPTIONAL:: t CHARACTER (LEN=iterationDigits):: tString CHARACTER(:), ALLOCATABLE:: fileName - IF (PRESENT(timeStep)) THEN - WRITE(tString, iterationFormat) timeStep + IF (PRESENT(t)) THEN + WRITE(tString, iterationFormat) t fileName = prefix // '_' // tString // '_' // suffix // '.' // extension ELSE diff --git a/src/modules/mesh/inout/text/makefile b/src/modules/mesh/inout/text/makefile deleted file mode 100644 index 26d0f10..0000000 --- a/src/modules/mesh/inout/text/makefile +++ /dev/null @@ -1,7 +0,0 @@ -all: moduleMeshInputText.o moduleMeshOutputText.o - -moduleMeshInputText.o: moduleMeshOutputText.o moduleMeshInputText.f90 - $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ - -%.o: %.f90 - $(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/inout/text/moduleMeshInputText.f90 b/src/modules/mesh/inout/text/moduleMeshInputText.f90 deleted file mode 100644 index a41e6d3..0000000 --- a/src/modules/mesh/inout/text/moduleMeshInputText.f90 +++ /dev/null @@ -1,232 +0,0 @@ -module moduleMeshInputText - !The mesh is stored as a column-wise text file. - !Aimed for simple geometries in 1D - - contains - !Inits the text mesh - subroutine initText(self) - use moduleMesh - use moduleMeshOutputText - implicit none - - class(meshGeneric), intent(inout), target:: self - - if (associated(meshForMCC,self)) then - self%printColl => printCollText - - end if - - select type(self) - type is (meshParticles) - self%printOutput => printOutputText - self%printEM => printEMText - self%printAverage => printAverageText - - self%readInitial => readInitialText - - end select - - self%readMesh => readText - - end subroutine initText - - !Reads the text mesh - subroutine readText(self, filename) - use moduleMesh - use moduleMesh1DCart - use moduleMesh1DRad - use moduleErrors - implicit none - - class(meshGeneric), intent(inout):: self - character(:), allocatable, intent(in):: filename !Dummy file, not used - integer:: fileID, reason - character(len=256):: line - integer:: nNodes - real(8):: r(1:3) !dummy 3D coordinate - integer:: physicalID - integer:: n, c - integer, allocatable:: p(:) - integer:: bt - - fileID = 10 - - open(fileID, file=trim(filename)) - - !Skip header - read(fileID, *) - - !Get number of nodes - nNodes = 0 - do - read(fileID, *, iostat=reason) line - - if (reason > 0) then - call criticalError('Error reading mesh file', 'readText') - - else if (reason < 0) then - exit - - else if (len(line) > 0) then - nNodes = nNodes + 1 - - end if - - end do - - if (nNodes == 0) then - call criticalError('No nodes read in mesh file', 'readText') - - end if - - self%numNodes = nNodes - allocate(self%nodes(1:self%numNodes)) - - SELECT TYPE(self) - TYPE IS(meshParticles) - ALLOCATE(self%K(1:self%numNodes, 1:self%numNodes)) - ALLOCATE(self%IPIV(1:self%numNodes, 1:self%numNodes)) - self%K = 0.D0 - self%IPIV = 0 - - END SELECT - - self%numCells = nNodes - 1 - allocate(self%cells(1:self%numCells)) - - select type(self) - type is (meshParticles) - self%numEdges = 2 - - allocate(self%edges(1:self%numEdges)) - - end select - - !Read the mesh now - rewind(fileID) - - !Skip header - read(fileID, *) - - !Allocate nodes and edges - do n = 1, self%numNodes - r = 0.D0 - - read(fileID, *) r(1), physicalID - - select case(self%geometry) - case("Cart") - allocate(meshNode1DCart:: self%nodes(n)%obj) - - case("Rad") - allocate(meshNode1DRad:: self%nodes(n)%obj) - - end select - - !Init nodes - call self%nodes(n)%obj%init(n, r) - - !Allocate edges if required) - select type(self) - type is (meshParticles) - if ((physicalID == 1) .or. (physicalID == 2)) then - select case(self%geometry) - case("Cart") - allocate(meshEdge1DCart:: self%edges(physicalID)%obj) - - case("Rad") - allocate(meshEdge1DRad:: self%edges(physicalID)%obj) - - end select - - allocate(p(1)) - p(1) = n - bt = getBoundaryId(physicalID) - call self%edges(physicalID)%obj%init(physicalID, p, physicalID, physicalID) - deallocate(p) - - end if - - end select - - end do - - !Allocate cells - n = 1 - allocate(p(1:2)) - do c = 1, self%numCells - p(1) = n - n = n + 1 - p(2) = n - - select case(self%geometry) - case("Cart") - allocate(meshCell1DCartSegm:: self%cells(c)%obj) - - case("Rad") - allocate(meshCell1DRadSegm:: self%cells(c)%obj) - - end select - - call self%cells(c)%obj%init(c, p, self%nodes) - - - end do - deallocate(p) - - close(fileID) - - !Call mesh connectivity - CALL self%connectMesh - - end subroutine readText - - subroutine readInitialText(filename, density, velocity, temperature) - use moduleErrors - implicit none - - character(:), allocatable, intent(in):: filename - real(8), allocatable, intent(out), dimension(:):: density - real(8), allocatable, intent(out), dimension(:,:):: velocity - real(8), allocatable, intent(out), dimension(:):: temperature - integer:: fileID, reason - character(len=256):: line - integer:: nNodes - integer:: n - - fileID = 10 - - open(fileID, file=trim(filename)) - - do - read(fileID, *, iostat=reason) line - - if (reason > 0) then - call criticalError('Error reading mesh file', 'readText') - - else if (reason < 0) then - exit - - else if (len(line) > 0) then - nNodes = nNodes + 1 - - end if - - end do - - allocate(density(1:nNodes)) - allocate(velocity(1:nNodes, 1:3)) - allocate(temperature(1:nNodes)) - - rewind(fileID) - - do n = 1, nNodes - read(fileID, *) density(n), velocity(n, 1:3), temperature(n) - - end do - - close(fileID) - - end subroutine readInitialText - -end module moduleMeshInputText diff --git a/src/modules/mesh/inout/text/moduleMeshOutputText.f90 b/src/modules/mesh/inout/text/moduleMeshOutputText.f90 deleted file mode 100644 index d69ddda..0000000 --- a/src/modules/mesh/inout/text/moduleMeshOutputText.f90 +++ /dev/null @@ -1,265 +0,0 @@ -module moduleMeshOutputText - contains - - subroutine writeSpeciesOutput(self, fileID, speciesIndex) - use moduleMesh - use moduleOutput - use moduleRefParam, only: L_ref - implicit none - - class(meshParticles), INTENT(in):: self - integer, intent(in):: fileID - integer, intent(in):: speciesIndex - real(8):: r(1:3) - type(outputFormat):: output - integer:: n - - do n = 1, self%numNodes - r = self%nodes(n)%obj%getCoordinates() - call calculateOutput(self%nodes(n)%obj%output(speciesIndex), output, self%nodes(n)%obj%v, species(speciesIndex)%obj) - - write(fileID, '(5(ES0.6E3,","),ES0.6E3)') r(1)*L_ref, output%density, output%velocity, output%temperature - - end do - - end subroutine writeSpeciesOutput - - subroutine writeCollOutput(self, fileID) - use moduleMesh - use moduleCollisions - use moduleRefParam, only: L_ref - implicit none - - class(meshGeneric), intent(in):: self - integer, intent(in):: fileID - integer:: n, k, c - - do n = 1, self%numCells - write(fileID, '(I0)', advance='no') n - - do k = 1, nCollPairs - do c = 1, interactionMatrix(k)%amount - write(fileID, '(",",I0)', advance='no') self%cells(n)%obj%tallyColl(k)%tally(c) - - end do - - end do - write(fileID, *) - - end do - - end subroutine writeCollOutput - - subroutine writeEMOutput(self, fileID) - use moduleMesh - use moduleRefParam, only: L_ref, Volt_ref, B_ref, EF_ref - implicit none - - class(meshParticles), intent(in):: self - integer, intent(in):: fileID - integer:: n, c - real(8):: r(1:3), Xi(1:3) - - do n = 1, self%numNodes - r = self%nodes(n)%obj%getCoordinates() - if (n == self%numNodes) then - Xi = (/ 1.D0, 0.D0, 0.D0 /) - c = self%numNodes - 1 - - else - Xi = (/ 0.D0, 0.D0, 0.D0 /) - c = n - - end if - - associate(output => self%nodes(n)%obj%emData) - write(fileID, '(7(ES0.6E3,","),ES0.6E3)') r(1)*L_ref, & - output%phi*Volt_ref, & - self%cells(c)%obj%gatherElectricField(Xi)*EF_ref, & - output%B*B_ref - - end associate - - end do - - end subroutine writeEMOutput - - subroutine writeAverage(self, fileIDMean, & - fileIDDeviation, & - speciesIndex) - - use moduleMesh - use moduleOutput - use moduleAverage - use moduleRefParam, only: L_ref - implicit none - - class(meshParticles), intent(in):: self - integer, intent(in):: fileIDMean, fileIDDeviation - INTEGER, intent(in):: speciesIndex - real(8):: r(1:3) - type(outputFormat):: outputMean - type(outputFormat):: outputDeviation - integer:: n - - do n = 1, self%numNodes - r = self%nodes(n)%obj%getCoordinates() - - call calculateOutput(averageScheme(n)%mean%output(speciesIndex), outputMean, & - self%nodes(n)%obj%v, species(speciesIndex)%obj) - - write(fileIDMean, '(5(ES0.6E3,","),ES0.6E3)') r(1)*L_ref, outputMean%density, outputMean%velocity, outputMean%temperature - - call calculateOutput(averageScheme(n)%deviation%output(speciesIndex), outputDeviation, & - self%nodes(n)%obj%v, species(speciesIndex)%obj) - - write(fileIDDeviation, '(5(ES0.6E3,","),ES0.6E3)') r(1)*L_ref, outputDeviation%density, outputDeviation%velocity, outputDeviation%temperature - - end do - - end subroutine writeAverage - - subroutine printOutputText(self) - use moduleMesh - use moduleSpecies - use moduleMeshInoutCommon - use moduleCaseParam, ONLY: timeStep - implicit none - - class(meshParticles), intent(in):: self - - INTEGER:: s, fileID - character(:), allocatable:: fileName - - fileID = 60 - - do s = 1, nSpecies - fileName = formatFileName(prefix, species(s)%obj%name, 'csv', timeStep) - write(*, "(6X,A15,A)") "Creating file: ", fileName - open (fileID, file = path // folder // '/' // fileName) - - write(fileID, '(5(A,","),A)') 'Position (m)', & - 'Density (m^-3)', & - 'Velocity (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', & - 'Temperature (K)' - - call writeSpeciesOutput(self, fileID, s) - - close(fileID) - - end do - - end subroutine printOutputText - - subroutine printCollText(self) - use moduleMesh - use moduleOutput - use moduleMeshInoutCommon - use moduleCaseParam, only: timeStep - implicit none - - class(meshGeneric), intent(in):: self - integer:: fileID - character(:), allocatable:: fileName - integer:: k, c - character (len=2):: cString - - fileID = 62 - - if (collOutput) then - fileName = formatFileName(prefix, 'Collisions', 'csv', timeStep) - write(*, "(6X,A15,A)") "Creating file: ", fileName - open (fileID, file = path // folder // '/' // fileName) - - write(fileID, '(A)', advance='no') "Cell" - do k = 1, nCollPairs - do c = 1, interactionMatrix(k)%amount - write(cString, "(I2)") c - write(fileID, '(",",A)', advance='no') 'Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString - - end do - end do - - write(fileID, *) - - call writeCollOutput(self, fileID) - - close(fileID) - - end if - - end subroutine printCollText - - subroutine printEMText(self) - use moduleMesh - use moduleMeshInoutCommon - use moduleCaseParam, only: timeStep - implicit none - - class(meshParticles), intent(in):: self - integer:: fileID - character(:), allocatable:: fileName - - fileID = 64 - - if (emOutput) then - fileName = formatFileName(prefix, 'EMField', 'csv', timeStep) - write(*, "(6X,A15,A)") "Creating file: ", fileName - open (fileID, file = path // folder // '/' // fileName) - - write(fileID, '(8(A,","),A)') 'Position (m)', & - 'Potential (V)', & - 'Electric Field (V m^-1):0', 'Electric Field (V m^-1):1', 'Electric Field (V m^-1):2', & - 'Magnetic Field (T):0', 'Magnetic Field (T):1', 'Magnetic Field (T):2' - - call writeEMOutput(self, fileID) - - close(fileID) - - end if - - end subroutine printEMText - - subroutine printAverageText(self) - use moduleMesh - use moduleSpecies - use moduleMeshInoutCommon - implicit none - - class(meshParticles), intent(in):: self - integer:: s - integer:: fileIDMean, fileIDDeviation - character(:), allocatable:: fileNameMean, fileNameDeviation - - fileIDMean = 66 - fileIDDeviation = 67 - - do s = 1, nSpecies - fileNameMean = formatFileName('Average_mean', species(s)%obj%name, 'csv', timeStep) - write(*, "(6X,A15,A)") "Creating file: ", fileNameMean - open (fileIDMean, file = path // folder // '/' // fileNameMean) - - write(fileIDMean, '(5(A,","),A)') 'Position (m)', & - 'Density, mean (m^-3)', & - 'Velocity, mean (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', & - 'Temperature, mean (K)' - - fileNameDeviation = formatFileName('Average_deviation', species(s)%obj%name, 'csv', timeStep) - write(*, "(6X,A15,A)") "Creating file: ", fileNameDeviation - open (fileIDDeviation, file = path // folder // '/' // fileNameDeviation) - - write(fileIDDeviation, '(5(A,","),A)') 'Position (m)', & - 'Density, deviation (m^-3)', & - 'Velocity, deviation (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', & - 'Temperature, deviation (K)' - - call writeAverage(self, fileIDMean, fileIDDeviation, s) - - close(fileIDMean) - close(fileIDDeviation) - - end do - - end subroutine printAverageText - -end module moduleMeshOutputText diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index c7b89b5..763517f 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -167,7 +167,7 @@ MODULE moduleMeshInputVTU CLASS(meshGeneric), INTENT(inout):: self CHARACTER(:), ALLOCATABLE, INTENT(in):: filename REAL(8):: r(1:3) !3 generic coordinates - INTEGER:: fileID + INTEGER:: fileID, error, found CHARACTER(LEN=256):: line INTEGER:: numNodes, numElements, numEdges INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types @@ -275,7 +275,6 @@ MODULE moduleMeshInputVTU END DO !Count the number of edges - numEdges = 0 SELECT CASE(self%dimen) CASE(3) !Edges are triangles, type 5 in VTK @@ -496,7 +495,7 @@ MODULE moduleMeshInputVTU END SELECT END DO - + !Call mesh connectivity CALL self%connectMesh @@ -548,8 +547,6 @@ MODULE moduleMeshInputVTU CALL readDataBlock(fileID, numNodes, temperature) REWIND(fileID) - close(fileID) - END SUBROUTINE readInitialVTU END MODULE moduleMeshInputVTU diff --git a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 index 8e8f91a..6286cfc 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 @@ -11,7 +11,7 @@ MODULE moduleMeshOutputVTU WRITE(fileID,"(A)") '' WRITE(fileID,"(2X, A)") '' - WRITE(fileID,"(4X, A)") '' + WRITE(fileID,"(4X, A,ES20.6E3,A)") '' WRITE(fileID,"(6X, A, I10, A, I10, A)") '' END SUBROUTINE writeHeader @@ -209,22 +209,23 @@ MODULE moduleMeshOutputVTU WRITE(fileID,"(8X,A)") '' !Electric field WRITE(fileID,"(10X,A, A, A)") '' - WRITE(fileID,"(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells) + WRITE(fileID, "(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells) WRITE(fileID,"(10X,A)") '' WRITE(fileID,"(8X,A)") '' END SUBROUTINE writeEM - SUBROUTINE writeCollection(fileID, fileNameStep, fileNameCollection) + SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection) USE moduleCaseParam USE moduleOutput USE moduleRefParam IMPLICIT NONE INTEGER:: fileID + INTEGER, INTENT(in):: t CHARACTER(*):: fileNameStep, fileNameCollection - IF (timeStep == tInitial) THEN + IF (t == tInitial) THEN !Create collection file WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection) @@ -236,11 +237,10 @@ MODULE moduleMeshOutputVTU !Write iteration file in collection OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection, ACCESS='APPEND') - WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") & - '' + WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") '' !Close collection file - IF (timeStep == tFinal) THEN + IF (t == tFinal) THEN WRITE (fileID + 1, "(2X, A)") '' WRITE (fileID + 1, "(A)") '' @@ -307,21 +307,22 @@ MODULE moduleMeshOutputVTU END SUBROUTINE writeAverage - SUBROUTINE printOutputVTU(self) + SUBROUTINE printOutputVTU(self,t) USE moduleMesh USE moduleSpecies USE moduleMeshInoutCommon - USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER:: i, fileID + INTEGER, INTENT(in):: t + INTEGER:: n, i, fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection + TYPE(outputFormat):: output(1:self%numNodes) fileID = 60 DO i = 1, nSpecies - fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', timeStep) + fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) @@ -337,27 +338,29 @@ MODULE moduleMeshOutputVTU !Write collection file for time plotting fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd') - CALL writeCollection(fileID, fileName, filenameCollection) + CALL writeCollection(fileID, t, fileName, filenameCollection) END DO END SUBROUTINE printOutputVTU - SUBROUTINE printCollVTU(self) + SUBROUTINE printCollVTU(self,t) USE moduleMesh USE moduleOutput USE moduleMeshInoutCommon - USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self - INTEGER:: fileID + INTEGER, INTENT(in):: t + INTEGER:: n, i, fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection + CHARACTER (LEN=iterationDigits):: tstring + TYPE(outputFormat):: output(1:self%numNodes) fileID = 62 IF (collOutput) THEN - fileName = formatFileName(prefix, 'Collisions', 'vtu', timeStep) + fileName = formatFileName(prefix, 'Collisions', 'vtu', t) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) @@ -373,26 +376,26 @@ MODULE moduleMeshOutputVTU !Write collection file for time plotting fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd') - CALL writeCollection(fileID, fileName, filenameCollection) + CALL writeCollection(fileID, t, fileName, filenameCollection) END IF END SUBROUTINE printCollVTU - SUBROUTINE printEMVTU(self) + SUBROUTINE printEMVTU(self, t) USE moduleMesh USE moduleMeshInoutCommon - USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: t INTEGER:: fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection fileID = 64 IF (emOutput) THEN - fileName = formatFileName(prefix, 'EMField', 'vtu', timeStep) + fileName = formatFileName(prefix, 'EMField', 'vtu', t) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) @@ -408,7 +411,7 @@ MODULE moduleMeshOutputVTU !Write collection file for time plotting fileNameCollection = formatFileName('Collection', 'EMField', 'pvd') - CALL writeCollection(fileID, fileName, filenameCollection) + CALL writeCollection(fileID, t, fileName, filenameCollection) END IF @@ -421,8 +424,9 @@ MODULE moduleMeshOutputVTU IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER:: i, fileIDMean, fileIDDeviation + INTEGER:: n, i, fileIDMean, fileIDDeviation CHARACTER(:), ALLOCATABLE:: fileNameMean, fileNameDeviation + TYPE(outputFormat):: output(1:self%numNodes) fileIDMean = 66 fileIDDeviation = 67 diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 01b03c7..203047f 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -59,13 +59,6 @@ MODULE moduleMesh END TYPE meshNodeCont - ! Array of pointers to nodes. - TYPE:: meshNodePointer - CLASS(meshNode), POINTER:: obj - CONTAINS - - END TYPE meshNodePointer - !Type for array of boundary functions (one per species) TYPE, PUBLIC:: fBoundaryGeneric PROCEDURE(boundary_interface), POINTER, NOPASS:: apply => NULL() @@ -83,8 +76,8 @@ MODULE moduleMesh CLASS(meshCell), POINTER:: eColl => NULL() !Normal vector REAL(8):: normal(1:3) - ! Surface of edge - REAL(8):: surface = 0.D0 + !Weight for random injection of particles + REAL(8):: weight = 1.D0 !Pointer to boundary type TYPE(boundaryCont), POINTER:: boundary !Array of functions for boundary conditions @@ -344,10 +337,10 @@ MODULE moduleMesh !Array of cell elements TYPE(meshCellCont), ALLOCATABLE:: cells(:) !PROCEDURES SPECIFIC OF FILE TYPE - PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL() - PROCEDURE(readInitial_interface), POINTER, NOPASS:: readInitial => NULL() - PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL() - PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL() + PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL() + PROCEDURE(readInitial_interface), POINTER, NOPASS:: readInitial => NULL() + PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL() + PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL() CONTAINS !GENERIC PROCEDURES PROCEDURE, PASS:: doCollisions @@ -379,9 +372,10 @@ MODULE moduleMesh END SUBROUTINE connectMesh_interface !Prints number of collisions in each cell - SUBROUTINE printColl_interface(self) + SUBROUTINE printColl_interface(self, t) IMPORT meshGeneric CLASS(meshGeneric), INTENT(in):: self + INTEGER, INTENT(in):: t END SUBROUTINE printColl_interface @@ -409,16 +403,18 @@ MODULE moduleMesh ABSTRACT INTERFACE !Prints Species data - SUBROUTINE printOutput_interface(self) + SUBROUTINE printOutput_interface(self, t) IMPORT meshParticles CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: t END SUBROUTINE printOutput_interface !Prints EM info - SUBROUTINE printEM_interface(self) + SUBROUTINE printEM_interface(self, t) IMPORT meshParticles CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: t END SUBROUTINE printEM_interface @@ -499,29 +495,28 @@ MODULE moduleMesh IMPLICIT NONE CLASS(meshParticles), INTENT(inout):: self - INTEGER:: c + INTEGER:: e + INTEGER:: nNodes INTEGER, ALLOCATABLE:: n(:) REAL(8), ALLOCATABLE:: localK(:,:) INTEGER:: i, j - DO c = 1, self%numCells - associate(nNodes => self%cells(c)%obj%nNodes) - ALLOCATE(n(1:nNodes)) - ALLOCATE(localK(1:nNodes, 1:nNodes)) - n = self%cells(c)%obj%getNodes(nNodes) - localK = self%cells(c)%obj%elemK(nNodes) + DO e = 1, self%numCells + nNodes = self%cells(e)%obj%nNodes + ALLOCATE(n(1:nNodes)) + ALLOCATE(localK(1:nNodes, 1:nNodes)) + n = self%cells(e)%obj%getNodes(nNodes) + localK = self%cells(e)%obj%elemK(nNodes) - DO i = 1, nNodes - DO j = 1, nNodes - self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j) - - END DO + DO i = 1, nNodes + DO j = 1, nNodes + self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j) END DO - DEALLOCATE(n, localK) + END DO - end associate + DEALLOCATE(n, localK) END DO @@ -618,7 +613,6 @@ MODULE moduleMesh INTEGER:: sp INTEGER:: i CLASS(meshNode), POINTER:: node - REAL(8):: pFraction !Particle fraction cellNodes = self%getNodes(nNodes) fPsi = self%fPsi(part%Xi, nNodes) @@ -629,11 +623,10 @@ MODULE moduleMesh DO i = 1, nNodes node => mesh%nodes(cellNodes(i))%obj - pFraction = fPsi(i)*part%weight CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + pFraction - node%output(sp)%mom(:) = node%output(sp)%mom(:) + pFraction*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + pFraction*tensorS + node%output(sp)%den = node%output(sp)%den + part%weight*fPsi(i) + node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*fPsi(i)*part%v(:) + node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*fPsi(i)*tensorS CALL OMP_UNSET_LOCK(node%lock) END DO @@ -794,7 +787,7 @@ MODULE moduleMesh END FUNCTION findCellBrute !Computes collisions in element - SUBROUTINE doCollisions(self) + SUBROUTINE doCollisions(self, t) USE moduleCollisions USE moduleSpecies USE moduleList @@ -802,10 +795,10 @@ MODULE moduleMesh USE moduleRandom USE moduleOutput USE moduleMath - USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshGeneric), INTENT(inout), TARGET:: self + INTEGER, INTENT(in):: t INTEGER:: e CLASS(meshCell), POINTER:: cell INTEGER:: k, i, j @@ -821,11 +814,10 @@ MODULE moduleMesh REAL(8):: rnd_real !Random number for collision INTEGER:: rnd_int !Random number for collision - IF (MOD(timeStep, everyColl) == 0) THEN + IF (MOD(t, everyColl) == 0) THEN !Collisions need to be performed in this iteration - !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j) + !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j, cell) DO e=1, self%numCells - cell => self%cells(e)%obj !TODO: Simplify this, to many sublevels @@ -892,7 +884,7 @@ MODULE moduleMesh !Obtain the cross sections for the different processes !TODO: From here it might be a procedure in interactionMatrix vRel = NORM2(part_i%v-part_j%v) - rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m) + rMass = reducedMass(part_i%weight*part_i%species%mass, part_j%weight*part_j%species%mass) eRel = rMass*vRel**2 CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel) @@ -918,9 +910,7 @@ MODULE moduleMesh !Loop over collisions DO c = 1, interactionMatrix(k)%amount IF (rnd_real <= probabilityColl(c)) THEN - !$OMP CRITICAL CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel) - !$OMP END CRITICAL !If collisions are gonna be output, count the collision IF (collOutput) THEN @@ -1030,9 +1020,6 @@ MODULE moduleMesh ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3)) ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3)) ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount)) - deltaV_ij = 0.D0 - p_ij = 0.D0 - mass_ij = 0.D0 !Loop over particles of species_i partTemp => cell%listPart_in(i)%head p = 1 @@ -1092,7 +1079,7 @@ MODULE moduleMesh !Compute changes in velocity for each particle deltaV_ij(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v - mass_ij(p) = pair%sp_i%m*partTemp%part%weight + mass_ij(p) = pair%sp_i%mass*partTemp%part%weight p_ij(p,1:3) = mass_ij(p)*partTemp%part%v !Move to the next particle in the list @@ -1117,9 +1104,6 @@ MODULE moduleMesh ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3)) ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3)) ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount)) - deltaV_ji = 0.D0 - p_ji = 0.D0 - mass_ji = 0.D0 !Loop over particles of species_j partTemp => cell%listPart_in(j)%head p = 1 @@ -1178,7 +1162,7 @@ MODULE moduleMesh !Compute changes in velocity for each particle deltaV_ji(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v - mass_ji(p) = pair%sp_j%m*partTemp%part%weight + mass_ji(p) = pair%sp_j%mass*partTemp%part%weight p_ji(p,1:3) = mass_ji(p)*partTemp%part%v !Move to the next particle in the list diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index deb4459..e7e4d65 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -77,20 +77,6 @@ MODULE moduleMeshBoundary END SUBROUTINE transparent - !Symmetry axis. Reflects particles. - !Although this function should never be called, it is set as a reflective boundary - !to properly deal with possible particles reaching a corner and selecting this boundary. - SUBROUTINE symmetryAxis(edge, part) - USE moduleSpecies - IMPLICIT NONE - - CLASS(meshEdge), INTENT(inout):: edge - CLASS(particle), INTENT(inout):: part - - CALL reflection(edge, part) - - END SUBROUTINE symmetryAxis - !Wall with temperature SUBROUTINE wallTemperature(edge, part) USE moduleSpecies @@ -139,7 +125,7 @@ MODULE moduleMeshBoundary SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj) TYPE IS(boundaryIonization) - mRel = reducedMass(bound%m0, part%species%m) + mRel = reducedMass(bound%m0, part%species%mass) vRel = SUM(DABS(part%v-bound%v0)) eRel = mRel*vRel**2*5.D-1 @@ -161,13 +147,7 @@ MODULE moduleMeshBoundary ALLOCATE(newElectron) ALLOCATE(newIon) - IF (ASSOCIATED(bound%electronSecondary)) THEN - newElectron%species => bound%electronSecondary - - ELSE - newElectron%species => part%species - - END IF + newElectron%species => part%species newIon%species => bound%species newElectron%v = v0 + (1.D0 + bound%deltaV*v0/NORM2(v0)) @@ -213,32 +193,134 @@ MODULE moduleMeshBoundary END SELECT - !Removes ionizing electron regardless the number of pair created + !Removes ionizing electron regardless the number of pairs created part%n_in = .FALSE. END SUBROUTINE ionization - subroutine outflowAdaptive(edge, part) - use moduleRandom - implicit none + !Symmetry axis. Reflects particles. + !Although this function should never be called, it is set as a reflective boundary + !to properly deal with possible particles reaching a corner and selecting this boundary. + SUBROUTINE symmetryAxis(edge, part) + USE moduleSpecies + IMPLICIT NONE - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part - select type(bound => edge%boundary%bTypes(part%species%n)%obj) - type is(boundaryOutflowAdaptive) + CALL reflection(edge, part) - if (random() < 0.844d0) then - call reflection(edge, part) + END SUBROUTINE symmetryAxis - else - call transparent(edge, part) + !Secondary emission of electrons by impacting particle. + SUBROUTINE secondaryEmission(edge, part) + USE moduleSpecies + USE moduleRandom + USE moduleConstParam + IMPLICIT NONE - end if + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + REAL(8):: vRel, eRel + REAL(8), DIMENSION(1:3):: rElectron, XiElectron!Position of new electrons (impacting particle) + REAL(8), DIMENSION(1:3):: rIncident !Vector from imapcting particle position to particle position + REAL(8):: theta !incident angle + REAL(8):: yield + REAL(8):: energy !incident energy corrected by threshold and + INTEGER:: nNewElectrons + REAL(8), ALLOCATABLE:: weight(:) + INTEGER:: p + INTEGER:: cell + TYPE(particle), POINTER:: secondaryElectron - end select + SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj) + TYPE IS(boundarySEE) + !Get relative velocity + vRel = NORM2(part%v) + !Convert to relative energy + eRel = part%species%mass*vRel**2*5.D-1 - end subroutine outflowAdaptive + !If energy is abound threshold calculate secondary electrons + IF (eRel >= bound%energyThreshold) THEN + + !position of impacting ion + rElectron = edge%intersection(part%r) + XiElectron = mesh%cells(part%cell)%obj%phy2log(rElectron) + + !Calculate incident angle + rIncident = part%r - rElectron + theta = ACOS(DOT_PRODUCT(rIncident, edge%normal) / (NORM2(rIncident) * NORM2(edge%normal))) + + !Calculate incident energy + energy = (eRel - bound%energyThreshold) * PI2 / (PI2 + theta**2) + bound%energyThreshold + + !Get number of secondary electrons particles + yield = part%weight*bound%yield%get(eRel) * (1.D0 * theta**2 / PI2) !Check equation! + + !Convert the number to macro-particles + nNewElectrons = FLOOR(yield / bound%electron%weight) + + !If the weight of new macro-particles is below the yield, correct adding a particle + IF (REAL(nNewElectrons) * bound%electron%weight < yield) THEN + nNewElectrons = nNewElectrons + 1 + ALLOCATE(weight(1:nNewElectrons)) + weight(1:nNewElectrons-1) = bound%electron%weight + weight(nNewElectrons) = yield - SUM(weight(1:nNewElectrons-1)) + + ELSE + ALLOCATE(weight(1:nNewElectrons)) + weight(1:nNewElectrons) = bound%electron%weight + + END IF + + !New cell of origin + IF (ASSOCIATED(edge%e1)) THEN + cell = edge%e1%n + + ELSEIF (ASSOCIATED(edge%e2)) THEN + cell = edge%e2%n + + END IF + + !Create the new electron macro-particles + DO p = 1, nNewElectrons + !Create new macro-particle + ALLOCATE(secondaryElectron) + + !Assign species to electron + secondaryElectron%species => bound%electron + + !Assign position to particle + secondaryElectron%r = rElectron + secondaryElectron%Xi = XiElectron + + !Assign cell to electron + secondaryElectron%cell = cell + + !Assign weight + secondaryElectron%weight = weight(p) + + !Assume particle is inside the numerical domain + secondaryElectron%n_in = .TRUE. + + !Assign velocity + secondaryElectron%v = 2.D0 * edge%normal + 1.D-1 * (/ randomMaxwellian(), randomMaxwellian(), randomMaxwellian() /) + + !Add particle to list + CALL partSurfaces%setLock() + CALL partSurfaces%add(secondaryElectron) + CALL partSurfaces%unsetLock() + + END DO + + END IF + + !Absorb impacting particle + CALL absorption(edge, part) + + END SELECT + + END SUBROUTINE secondaryEmission !Points the boundary function to specific type SUBROUTINE pointBoundaryFunction(edge, s) @@ -258,17 +340,17 @@ MODULE moduleMeshBoundary TYPE IS(boundaryTransparent) edge%fBoundary(s)%apply => transparent - TYPE IS(boundaryAxis) - edge%fBoundary(s)%apply => symmetryAxis - TYPE IS(boundaryWallTemperature) edge%fBoundary(s)%apply => wallTemperature TYPE IS(boundaryIonization) edge%fBoundary(s)%apply => ionization - type is(boundaryOutflowAdaptive) - edge%fBoundary(s)%apply => outflowAdaptive + TYPE IS(boundaryAxis) + edge%fBoundary(s)%apply => symmetryAxis + + TYPE IS(boundarySEE) + edge%fBoundary(s)%apply => secondaryEmission CLASS DEFAULT CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction') diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index 67465c7..b3c1d12 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -26,12 +26,6 @@ MODULE moduleBoundary END TYPE boundaryTransparent - !Symmetry axis - TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis - CONTAINS - - END TYPE boundaryAxis - !Wall Temperature boundary TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature !Thermal velocity of the wall: square root(Wall temperature X specific heat) @@ -44,7 +38,6 @@ MODULE moduleBoundary TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals. CLASS(speciesGeneric), POINTER:: species !Ion species - CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron TYPE(table1D):: crossSection REAL(8):: effectiveTime REAL(8):: eThreshold @@ -53,13 +46,21 @@ MODULE moduleBoundary END TYPE boundaryIonization - !Boundary for quasi-neutral outflow adjusting reflection coefficient - type, public, extends(boundaryGeneric):: boundaryOutflowAdaptive - real(8):: outflowCurrent - real(8):: reflectionFraction - contains + !Secondary electron emission (by ion impact) + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundarySEE + !Yield as a function of ion energy + TYPE(table1D):: yield + CLASS(speciesGeneric), POINTER:: electron !Electron species for secondary emission + REAL(8):: energyThreshold + CONTAINS - end type boundaryOutflowAdaptive + END TYPE boundarySEE + + !Symmetry axis + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis + CONTAINS + + END TYPE boundaryAxis !Wrapper for boundary types (one per species) TYPE:: bTypesCont @@ -112,19 +113,17 @@ MODULE moduleBoundary END SUBROUTINE initWallTemperature - SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary) + SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, speciesID, effTime, crossSection, eThreshold) USE moduleRefParam USE moduleSpecies USE moduleCaseParam USE moduleConstParam - USE moduleErrors IMPLICIT NONE CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary REAL(8), INTENT(in):: me !Electron mass REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties - INTEGER, INTENT(in):: ion - INTEGER, OPTIONAL, INTENT(in):: electronSecondary + INTEGER:: speciesID REAL(8):: effTime CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection REAL(8), INTENT(in):: eThreshold @@ -137,22 +136,7 @@ MODULE moduleBoundary boundary%n0 = n0 * Vol_ref boundary%v0 = v0 / v_ref boundary%vTh = DSQRT(kb*T0/m0)/v_ref - boundary%species => species(ion)%obj - IF (PRESENT(electronSecondary)) THEN - SELECT TYPE(sp => species(electronSecondary)%obj) - TYPE IS(speciesCharged) - boundary%electronSecondary => sp - - CLASS DEFAULT - CALL criticalError("Species " // sp%name // " chosen for " // & - "secondary electron is not a charged species", 'initIonization') - - END SELECT - - ELSE - boundary%electronSecondary => NULL() - - END IF + boundary%species => species(speciesID)%obj boundary%effectiveTime = effTime / ti_ref CALL boundary%crossSection%init(crossSection) CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) @@ -163,4 +147,36 @@ MODULE moduleBoundary END SUBROUTINE initIonization + SUBROUTINE initSEE(boundary, tableFile, speciesID) + USE moduleRefParam + USE moduleConstParam + USE moduleSpecies + IMPLICIT NONE + + CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary + CHARACTER(:), ALLOCATABLE, INTENT(in):: tableFile + INTEGER:: speciesID + INTEGER:: i + + ALLOCATE(boundarySEE:: boundary) + SELECT TYPE(boundary) + TYPE IS(boundarySEE) + CALL boundary%yield%init(tableFile) + CALL boundary%yield%convert(eV2J/(m_ref*v_ref**2), 1.D0) + boundary%electron => species(speciesID)%obj + !Search for the threshold energy in the table + DO i = 1, SIZE(boundary%yield%f) + IF (boundary%yield%f(i) > 0.D0) THEN + boundary%energyThreshold = boundary%yield%x(i) + + EXIT + + END IF + + END DO + + END SELECT + + END SUBROUTINE initSEE + END MODULE moduleBoundary diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index d9a6ef2..f01c20b 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -43,8 +43,7 @@ MODULE moduleCollisions TYPE, EXTENDS(collisionBinary):: collisionBinaryIonization REAL(8):: eThreshold !Minimum energy (non-dimensional units) required for ionization REAL(8):: deltaV !Change in velocity due to exchange of eThreshold - CLASS(speciesCharged), POINTER:: electron !Pointer to species considerer as electrons - CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron + CLASS(speciesCharged), POINTER:: electron !Pointer to species considerer as electrons CONTAINS PROCEDURE, PASS:: collide => collideBinaryIonization @@ -165,8 +164,8 @@ MODULE moduleCollisions self%amount = amount - mass_i = species(i)%obj%m - mass_j = species(j)%obj%m + mass_i = species(i)%obj%mass + mass_j = species(j)%obj%mass ALLOCATE(self%collisions(1:self%amount)) @@ -228,8 +227,8 @@ MODULE moduleCollisions REAL(8):: m_i, m_j REAL(8), DIMENSION(1:3):: vCM, vp - m_i = part_i%species%m - m_j = part_j%species%m + m_i = part_i%species%mass + m_j = part_j%species%mass !Applies the collision vCM = velocityCM(part_i%weight*m_i, part_i%v, part_j%weight*m_j, part_j%v) vp = vRel*randomDirectionVHS() @@ -242,7 +241,7 @@ MODULE moduleCollisions !ELECTRON IMPACT IONIZATION !Inits electron impact ionization - SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, electron, electronSecondary) + SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, electron) USE moduleTable USE moduleRefParam USE moduleConstParam @@ -254,8 +253,7 @@ MODULE moduleCollisions CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename REAL(8), INTENT(in):: energyThreshold CHARACTER(:), ALLOCATABLE, INTENT(in):: electron - CHARACTER(:), ALLOCATABLE, OPTIONAL, INTENT(in):: electronSecondary - INTEGER:: electronIndex, electronSecondaryIndex + INTEGER:: electronIndex ALLOCATE(collisionBinaryIonization:: collision) @@ -280,29 +278,12 @@ MODULE moduleCollisions CLASS DEFAULT CALL criticalError("Species " // sp%name // " chosen for " // & - "impacting electron is not a charged species", 'initBinaryIonization') + "secondary electron is not a charged species", 'initBinaryIonization') END SELECT - IF (PRESENT(electronSecondary)) THEN - electronSecondaryIndex = speciesName2Index(electronSecondary) - SELECT TYPE(sp => species(electronSecondaryIndex)%obj) - TYPE IS(speciesCharged) - collision%electronSecondary => sp - - CLASS DEFAULT - CALL criticalError("Species " // sp%name // " chosen for " // & - "secondary electron is not a charged species", 'initBinaryIonization') - - END SELECT - - ELSE - collision%electronSecondary => NULL() - - END IF - !momentum change per ionization process - collision%deltaV = sqrt(collision%eThreshold / collision%electron%m) + collision%deltaV = sqrt(collision%eThreshold / collision%electron%mass) END SELECT @@ -326,7 +307,7 @@ MODULE moduleCollisions REAL(8), DIMENSION(1:3):: vChange TYPE(particle), POINTER:: newElectron => NULL(), remainingNeutral => NULL() - rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m) + rMass = reducedMass(part_i%weight*part_i%species%mass, part_j%weight*part_j%species%mass) eRel = rMass*vRel**2 !Relative energy must be higher than threshold IF (eRel > self%eThreshold) THEN @@ -355,12 +336,6 @@ MODULE moduleCollisions !Copy basic information from primary electron newElectron = electron - !If secondary electron species indicates, convert - IF (ASSOCIATED(self%electronSecondary)) THEN - newElectron%species => self%electronSecondary - - END IF - !Secondary electorn gains energy from ionization newElectron%v = vChange @@ -387,7 +362,7 @@ MODULE moduleCollisions CALL sp%ionize(neutral) CLASS DEFAULT - CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization') + ! CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization') RETURN END SELECT diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90 index 166b6b2..da03ad8 100644 --- a/src/modules/moduleCoulomb.f90 +++ b/src/modules/moduleCoulomb.f90 @@ -61,7 +61,7 @@ MODULE moduleCoulomb self%sp_i => species(i)%obj self%sp_j => species(j)%obj - self%one_plus_massRatio_ij = 1.D0 + self%sp_i%m/self%sp_j%m + self%one_plus_massRatio_ij = 1.D0 + self%sp_i%mass/self%sp_j%mass Z_i = 0.D0 Z_j = 0.D0 @@ -87,11 +87,11 @@ MODULE moduleCoulomb scaleFactor = (n_ref * qe**4 * ti_ref) / (eps_0**2 * m_ref**2 * v_ref**3) - self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%m**2) * scaleFactor !Missing density because it's cell dependent - self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%m**2) * scaleFactor !Missing density because it's cell dependent + self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%mass**2) * scaleFactor !Missing density because it's cell dependent + self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%mass**2) * scaleFactor !Missing density because it's cell dependent - self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent - self%l2_i = self%sp_i%m / 2.D0 !Missing temperature because it's cell dependent + self%l2_j = self%sp_j%mass / 2.D0 !Missing temperature because it's cell dependent + self%l2_i = self%sp_i%mass / 2.D0 !Missing temperature because it's cell dependent END SUBROUTINE initInteractionCoulomb diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index ff8a694..b5babfb 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -54,16 +54,15 @@ MODULE moduleInject INTEGER:: id CHARACTER(:), ALLOCATABLE:: name REAL(8):: vMod !Velocity (module) - REAL(8):: temperature(1:3) !Temperature + REAL(8):: T(1:3) !Temperature REAL(8):: n(1:3) !Direction of injection LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n INTEGER:: nParticles !Number of particles to introduce each time step CLASS(speciesGeneric), POINTER:: species !Species of injection INTEGER:: nEdges INTEGER, ALLOCATABLE:: edges(:) !Array with edges - INTEGER, ALLOCATABLE:: particlesPerEdge(:) ! Particles per edge - REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge - REAL(8):: surface ! Total surface of injection + REAL(8), ALLOCATABLE:: cumWeight(:) !Array of cummulative probability + REAL(8):: sumWeight TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction CONTAINS PROCEDURE, PASS:: init => initInject @@ -76,7 +75,7 @@ MODULE moduleInject CONTAINS !Initialize an injection of particles - SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, physicalSurface, particlesPerEdge) + SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface) USE moduleMesh USE moduleRefParam USE moduleConstParam @@ -87,29 +86,49 @@ MODULE moduleInject CLASS(injectGeneric), INTENT(inout):: self INTEGER, INTENT(in):: i - REAL(8), INTENT(in):: v, n(1:3), temperature(1:3) - INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge + REAL(8), INTENT(in):: v, n(1:3), T(1:3) + INTEGER, INTENT(in):: sp, physicalSurface REAL(8):: tauInject REAL(8), INTENT(in):: flow CHARACTER(:), ALLOCATABLE, INTENT(in):: units INTEGER:: e, et INTEGER:: phSurface(1:mesh%numEdges) INTEGER:: nVolColl - REAL(8):: fluxPerStep = 0.D0 - self%id = i - self%vMod = v / v_ref - self%n = n / NORM2(n) - self%temperature = temperature / T_ref + self%id = i + self%vMod = v / v_ref + self%n = n / NORM2(n) + self%T = T / T_ref + self%species => species(sp)%obj + tauInject = tau(self%species%n) + SELECT CASE(units) + CASE ("sccm") + !Standard cubic centimeter per minute + self%nParticles = INT(flow*sccm2atomPerS*tauInject*ti_ref/species(sp)%obj%weight) + + CASE ("A") + !Input current in Ampers + self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight)) + + CASE ("part/s") + !Input current in Ampers + self%nParticles = INT(flow*tauInject*ti_ref/species(sp)%obj%weight) + + CASE DEFAULT + CALL criticalError("No support for units: " // units, 'initInject') + + END SELECT + !Scale particles for different species steps + IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') + !Gets the edge elements from which particles are injected DO e = 1, mesh%numEdges phSurface(e) = mesh%edges(e)%obj%physicalSurface END DO + self%nEdges = COUNT(phSurface == physicalSurface) - ALLOCATE(self%edges(1:self%nEdges)) - ALLOCATE(self%particlesPerEdge(1:self%nEdges)) - ALLOCATE(self%weightPerEdge(1:self%nEdges)) + ALLOCATE(inject(i)%edges(1:self%nEdges)) et = 0 DO e=1, mesh%numEdges IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN @@ -141,82 +160,15 @@ MODULE moduleInject END DO - !Calculates total area - self%surface = 0.D0 - DO et = 1, self%nEdges - self%surface = self%surface + mesh%edges(self%edges(et))%obj%surface + !Calculates cumulative probability + ALLOCATE(self%cumWeight(1:self%nEdges)) + et = 1 + self%cumWeight(1) = mesh%edges(self%edges(et))%obj%weight + DO et = 2, self%nEdges + self%cumWeight(et) = mesh%edges(self%edges(et))%obj%weight + self%cumWeight(et-1) END DO - - ! Information about species and flux - self%species => species(sp)%obj - tauInject = tau(self%species%n) - ! Convert units - SELECT CASE(units) - CASE ("sccm") - !Standard cubic centimeter per minute - fluxPerStep = flow*sccm2atomPerS - - CASE ("A") - !Current in Ampers - SELECT TYPE(sp => self%species) - CLASS IS(speciesCharged) - fluxPerStep = flow/(qe*abs(sp%q)) - - CLASS DEFAULT - call criticalError('Attempted to assign a flux in "A" to a species without charge.', 'initInject') - - END SELECT - - CASE ("Am2") - !Input current in Ampers per square meter - SELECT TYPE(sp => self%species) - CLASS IS(speciesCharged) - fluxPerStep = flow*self%surface*L_ref**2/(qe*abs(sp%q)) - - CLASS DEFAULT - call criticalError('Attempted to assign a flux in "Am2" to a species without charge.', 'initInject') - - END SELECT - - - CASE ("part/s") - !Input current in Ampers - fluxPerStep = flow - - CASE DEFAULT - CALL criticalError("No support for units: " // units, 'initInject') - - END SELECT - fluxPerStep = fluxPerStep * tauInject * ti_ref / self%surface - - !Assign particles per edge - IF (particlesPerEdge > 0) THEN - ! Particles per edge defined by the user - self%particlesPerEdge = particlesPerEdge - DO et = 1, self%nEdges - self%weightPerEdge(et) = fluxPerStep*mesh%edges(self%edges(et))%obj%surface / REAL(particlesPerEdge) - - END DO - - self%nParticles = SUM(self%particlesPerEdge) - - ELSE - ! No particles assigned per edge, use the species weight - self%weightPerEdge = self%species%weight - DO et = 1, self%nEdges - self%particlesPerEdge(et) = max(1,FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface / self%species%weight)) - END DO - - self%nParticles = SUM(self%particlesPerEdge) - - !Rescale weight to match flux - self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles)) - - END IF - - !Scale particles for different species steps - IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') + self%sumWeight = self%cumWeight(self%nEdges) END SUBROUTINE initInject @@ -227,16 +179,21 @@ MODULE moduleInject IMPLICIT NONE INTEGER:: i + INTEGER, DIMENSION(1:nInject):: nMin, nMax !$OMP SINGLE nPartInj = 0 DO i = 1, nInject IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN + nMin(i) = nPartInj + 1 nPartInj = nPartInj + inject(i)%nParticles + nMax(i) = nPartInj END IF END DO + PRINT *, nMin + PRINT *, nMax IF (ALLOCATED(partInj)) DEALLOCATE(partInj) ALLOCATE(partInj(1:nPartInj)) @@ -244,30 +201,30 @@ MODULE moduleInject DO i=1, nInject IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN - CALL inject(i)%addParticles() + CALL inject(i)%addParticles(nMin(i), nMax(i)) END IF END DO END SUBROUTINE doInjects - SUBROUTINE initVelDistMaxwellian(velDist, temperature, m) + SUBROUTINE initVelDistMaxwellian(velDist, T, m) IMPLICIT NONE CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist - REAL(8), INTENT(in):: temperature, m + REAL(8), INTENT(in):: T, m - velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m)) + velDist = velDistMaxwellian(vTh = DSQRT(T/m)) END SUBROUTINE initVelDistMaxwellian - SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m) + SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m) IMPLICIT NONE CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist - REAL(8), INTENT(in):: temperature, m + REAL(8), INTENT(in):: T, m - velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m)) + velDist = velDistHalfMaxwellian(vTh = DSQRT(T/m)) END SUBROUTINE initVelDistHalfMaxwellian @@ -289,7 +246,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = self%vTh*randomMaxwellian()/sqrt(2.d0) + v = self%vTh*randomMaxwellian() END FUNCTION randomVelMaxwellian @@ -302,7 +259,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = self%vTh*randomHalfMaxwellian()/sqrt(2.d0) + v = self%vTh*randomHalfMaxwellian() END FUNCTION randomVelHalfMaxwellian @@ -327,82 +284,68 @@ MODULE moduleInject IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self - INTEGER, SAVE:: nMin - INTEGER:: i, e + INTEGER, INTENT(in):: nMin, nMax !Min and Max index in partInj array + INTEGER:: randomX + INTEGER:: i INTEGER:: n, sp CLASS(meshEdge), POINTER:: randomEdge REAL(8):: direction(1:3) !Insert particles !$OMP SINGLE - nMin = 0 - DO i = 1, self%id -1 - IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN - nMin = nMin + inject(i)%nParticles - - END IF - - END DO - nMin = nMin + 1 + !Assign weight to particle. + partInj(nMin:nMax)%weight = self%species%weight + !Particle is considered to be outside the domain + partInj(nMin:nMax)%n_in = .FALSE. !$OMP END SINGLE !$OMP DO - DO e = 1, self%nEdges - ! Select edge for injection - randomEdge => mesh%edges(self%edges(e))%obj - ! Inject particles in edge - DO i = 1, self%particlesPerEdge(e) - ! Index in the global partInj array - n = nMin - 1 + SUM(self%particlesPerEdge(1:e-1)) + i - !Particle is considered to be outside the domain - partInj(n)%n_in = .FALSE. - !Random position in edge - partInj(n)%r = randomEdge%randPos() - !Assign weight to particle. - partInj(n)%weight = self%weightPerEdge(e) - !Volume associated to the edge: - IF (ASSOCIATED(randomEdge%e1)) THEN - partInj(n)%cell = randomEdge%e1%n + DO n = nMin, nMax + randomX = randomWeighted(self%cumWeight, self%sumWeight) - ELSEIF (ASSOCIATED(randomEdge%e2)) THEN - partInj(n)%cell = randomEdge%e2%n + randomEdge => mesh%edges(self%edges(randomX))%obj + !Random position in edge + partInj(n)%r = randomEdge%randPos() + !Volume associated to the edge: + IF (ASSOCIATED(randomEdge%e1)) THEN + partInj(n)%cell = randomEdge%e1%n - ELSE - CALL criticalError("No Volume associated to edge", 'addParticles') + ELSEIF (ASSOCIATED(randomEdge%e2)) THEN + partInj(n)%cell = randomEdge%e2%n - END IF - partInj(n)%cellColl = randomEdge%eColl%n - sp = self%species%n + ELSE + CALL criticalError("No Volume associated to edge", 'addParticles') - !Assign particle type - partInj(n)%species => self%species + END IF + partInj(n)%cellColl = randomEdge%eColl%n + sp = self%species%n - if (all(self%n == 0.D0)) then - direction = randomEdge%normal + !Assign particle type + partInj(n)%species => self%species - else - direction = self%n + direction = self%n - end if + !Sample initial velocity + partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & + self%v(2)%obj%randomVel(), & + self%v(3)%obj%randomVel() /) - partInj(n)%v = 0.D0 + !For each direction, velocities have to agree with the direction of injection + DO i = 1, 3 + DO WHILE (partInj(n)%v(i)*direction(i) < 0) + partInj(n)%v(i) = self%vMod*direction(i) + self%v(i)%obj%randomVel() - do while(dot_product(partInj(n)%v, direction) <= 0.d0) - partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & - self%v(2)%obj%randomVel(), & - self%v(3)%obj%randomVel() /) - end do - - - !Obtain natural coordinates of particle in cell - partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) - !Push new particle with the minimum time step - CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp)) - !Assign cell to new particle - CALL solver%updateParticleCell(partInj(n)) + END DO END DO + !Obtain natural coordinates of particle in cell + partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) + !Push new particle with the minimum time step + CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp)) + !Assign cell to new particle + CALL solver%updateParticleCell(partInj(n)) + END DO !$OMP END DO diff --git a/src/modules/moduleProbe.f90 b/src/modules/moduleProbe.f90 index 754d56a..c7d3cf5 100644 --- a/src/modules/moduleProbe.f90 +++ b/src/modules/moduleProbe.f90 @@ -27,7 +27,7 @@ MODULE moduleProbe CONTAINS !Functions for probeDistFunc type - SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, everyTimeStep) + SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, timeStep) USE moduleCaseParam USE moduleRefParam USE moduleSpecies @@ -41,7 +41,7 @@ MODULE moduleProbe REAL(8), INTENT(in):: r(1:3) REAL(8), INTENT(in):: v1(1:2), v2(1:2), v3(1:2) INTEGER, INTENT(in):: points(1:3) - REAL(8), INTENT(in):: everyTimeStep + REAL(8), INTENT(in):: timeStep INTEGER:: sp, i REAL(8):: dv(1:3) @@ -91,17 +91,17 @@ MODULE moduleProbe 1:self%nv(3))) !Number of iterations between output - IF (everyTimeStep == 0.D0) THEN + IF (timeStep == 0.D0) THEN self%every = 1 ELSE - self%every = NINT(everyTimeStep/ tauMin / ti_ref) + self%every = NINT(timeStep/ tauMin / ti_ref) END IF !Maximum radius !TODO: Make this an input parameter - self%maxR = 1.D-2/L_ref + self%maxR = 1.D0 !Init the probe lock CALL OMP_INIT_LOCK(self%lock) @@ -148,7 +148,7 @@ MODULE moduleProbe deltaR = NORM2(self%r - part%r) !Only include particle if it is inside the maximum radius - ! IF (deltaR < self%maxR) THEN + IF (deltaR < self%maxR) THEN !find lower index for all dimensions CALL self%findLowerIndex(part%v, i, j, k, inside) @@ -162,40 +162,40 @@ MODULE moduleProbe fk = self%vk(k+1) - part%v(3) fk1 = part%v(3) - self%vk(k) - weight = part%weight * DEXP(-deltaR/self%maxR) - ! weight = part%weight + ! weight = part%weight * DEXP(deltaR/self%maxR) + weight = part%weight !Lock the probe CALL OMP_SET_LOCK(self%lock) !Assign particle weight to distribution function - self%f(i , j , k ) = self%f(i , j , k ) + fi * fj * fk * weight - self%f(i+1, j , k ) = self%f(i+1, j , k ) + fi1 * fj * fk * weight - self%f(i , j+1, k ) = self%f(i , j+1, k ) + fi * fj1 * fk * weight - self%f(i+1, j+1, k ) = self%f(i+1, j+1, k ) + fi1 * fj1 * fk * weight - self%f(i , j , k+1) = self%f(i , j , k+1) + fi * fj * fk1 * weight - self%f(i+1, j , k+1) = self%f(i+1, j , k+1) + fi1 * fj * fk1 * weight - self%f(i , j+1, k+1) = self%f(i , j+1, k+1) + fi * fj1 * fk1 * weight - self%f(i+1, j+1, k+1) = self%f(i+1, j+1, k+1) + fi1 * fj1 * fk1 * weight + self%f(i , j , k ) = fi * fj * fk * weight + self%f(i+1, j , k ) = fi1 * fj * fk * weight + self%f(i , j+1, k ) = fi * fj1 * fk * weight + self%f(i+1, j+1, k ) = fi1 * fj1 * fk * weight + self%f(i , j , k+1) = fi * fj * fk1 * weight + self%f(i+1, j , k+1) = fi1 * fj * fk1 * weight + self%f(i , j+1, k+1) = fi * fj1 * fk1 * weight + self%f(i+1, j+1, k+1) = fi1 * fj1 * fk1 * weight !Unlock the probe CALL OMP_UNSET_LOCK(self%lock) END IF - ! END IF + END IF END IF END SUBROUTINE calculate - SUBROUTINE output(self) + SUBROUTINE output(self, t) USE moduleOutput USE moduleRefParam - USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(probeDistFunc), INTENT(inout):: self + INTEGER, INTENT(in):: t CHARACTER (LEN=iterationDigits):: tstring CHARACTER (LEN=3):: pstring CHARACTER(:), ALLOCATABLE:: filename @@ -204,14 +204,14 @@ MODULE moduleProbe !Divide by the velocity cube volume self%f = self%f * self%dvInv - WRITE(tstring, iterationFormat) timeStep + WRITE(tstring, iterationFormat) t WRITE(pstring, "(I3.3)") self%id fileName='Probe_' // tstring// '_f_' // pstring // '.dat' WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (10, file = path // folder // '/' // fileName) WRITE(10, "(A1, 1X, A)") "# ", self%species%name WRITE(10, "(A6, 3(ES15.6E3), A2)") "# r = ", self%r(:)*L_ref, " m" - WRITE(10, "(A6, ES15.6E3, A2)") "# t = ", REAL(timeStep)*tauMin*ti_ref, " s" + WRITE(10, "(A6, ES15.6E3, A2)") "# t = ", REAL(t)*tauMin*ti_ref, " s" WRITE(10, "(A1, A19, 3(A20))") "#", "v1 (m s^-1)", "v2 (m s^-1)", "v3 (m s^-1)", "f" DO i = 1, self%nv(1) DO j = 1, self%nv(2) @@ -252,15 +252,15 @@ MODULE moduleProbe END SUBROUTINE doProbes - SUBROUTINE outputProbes() - USE moduleCaseParam, ONLY: timeStep + SUBROUTINE outputProbes(t) IMPLICIT NONE + INTEGER, INTENT(in):: t INTEGER:: i DO i = 1, nProbes IF (probe(i)%update) THEN - CALL probe(i)%output() + CALL probe(i)%output(t) END IF @@ -268,15 +268,15 @@ MODULE moduleProbe END SUBROUTINE outputProbes - SUBROUTINE resetProbes() - USE moduleCaseParam, ONLY: timeStep + SUBROUTINE resetProbes(t) IMPLICIT NONE + INTEGER, INTENT(in):: t INTEGER:: i DO i = 1, nProbes probe(i)%f = 0.D0 - probe(i)%update = timeStep == tFinal .OR. timeStep == tInitial .OR. MOD(timeStep, probe(i)%every) == 0 + probe(i)%update = t == tFinal .OR. t == tInitial .OR. MOD(t, probe(i)%every) == 0 END DO diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index ab08f08..d039247 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -4,12 +4,54 @@ MODULE moduleSpecies USE OMP_LIB IMPLICIT NONE + !Basic type that defines a macro-particle + TYPE:: particle + REAL(8):: r(1:3) !Position + REAL(8):: v(1:3) !Velocity + CLASS(speciesGeneric), POINTER:: species !Pointer to particle's species + INTEGER:: cell !Index of element in which the particle is located TODO: Make these pointers + INTEGER:: cellColl !Index of element in which the particle is located in the Collision Mesh + REAL(8):: Xi(1:3) !Logical coordinates of particle in cell + LOGICAL:: n_in !Flag that indicates if a particle is in the domain + REAL(8):: weight=0.D0 !weight of particle + + END TYPE particle + + !Wrapper to store the particles per species + TYPE:: particleArray + TYPE(particle), ALLOCATABLE, DIMENSION(:):: p + + END TYPE particleArray + + !Array of pointers for the particles to be pushed + TYPE:: particlePointer + TYPE(particle), POINTER:: p + + END TYPE particlePointer + + !Arrays that contain the particles + TYPE(particleArray), ALLOCATABLE, TARGET, DIMENSION(:):: particles !array of particles in the domain + TYPE(particle), ALLOCATABLE, TARGET, DIMENSION(:):: particlesInjection !array of inject particles + TYPE(particlePointer), ALLOCATABLE, DIMENSION(:):: particlesToPush !particles pushed in each iteration + !Integers to track number of particles in domain + INTEGER, ALLOCATABLE, DIMENSION(:):: nParticles + INTEGER:: nParticlesTotal + !Number of injected particles + INTEGER, ALLOCATABLE, DIMENSION(:):: nParticlesInject + INTEGER:: nPariclesInjectTotal + + !Generic species type TYPE, ABSTRACT:: speciesGeneric - CHARACTER(:), ALLOCATABLE:: name - REAL(8):: m=0.D0, weight=0.D0, qm=0.D0 - INTEGER:: n=0 + INTEGER:: n=0 !Index of species + CHARACTER(:), ALLOCATABLE:: name !Name of species + !Mass, default weight of species and charge over mass + REAL(8):: mass=0.D0, weight=0.D0, qm=0.D0 + INTEGER:: every !How many interations between advancing the species + LOGICAL:: pushSpecies !Boolean to indicate if the species is moved in the iteration + END TYPE speciesGeneric + !Neutral species TYPE, EXTENDS(speciesGeneric):: speciesNeutral CLASS(speciesGeneric), POINTER:: ion => NULL() CONTAINS @@ -17,6 +59,7 @@ MODULE moduleSpecies END TYPE speciesNeutral + !Charged species TYPE, EXTENDS(speciesGeneric):: speciesCharged REAL(8):: q=0.D0 CLASS(speciesGeneric), POINTER:: ion => NULL(), neutral => NULL() @@ -26,34 +69,17 @@ MODULE moduleSpecies END TYPE speciesCharged + !Wrapper for species TYPE:: speciesCont CLASS(speciesGeneric), ALLOCATABLE:: obj END TYPE + !Number of species INTEGER:: nSpecies + !Array for species TYPE(speciesCont), ALLOCATABLE, TARGET:: species(:) - TYPE particle - REAL(8):: r(1:3) !Position - REAL(8):: v(1:3) !Velocity - CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle - INTEGER:: cell !Index of element in which the particle is located - INTEGER:: cellColl !Index of element in which the particle is located in the Collision Mesh - REAL(8):: Xi(1:3) !Logical coordinates of particle in element e_p. - LOGICAL:: n_in !Flag that indicates if a particle is in the domain - REAL(8):: weight=0.D0 !weight of particle - - END TYPE particle - - !Number of old particles - INTEGER:: nPartOld - !Number of injected particles - INTEGER:: nPartInj - !Arrays that contain the particles - TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration - TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles - CONTAINS FUNCTION speciesName2Index(speciesName) RESULT(sp) USE moduleErrors diff --git a/src/modules/output/moduleOutput.f90 b/src/modules/output/moduleOutput.f90 index 2af3c70..16ed1b4 100644 --- a/src/modules/output/moduleOutput.f90 +++ b/src/modules/output/moduleOutput.f90 @@ -22,6 +22,7 @@ MODULE moduleOutput !Type for EM data in node TYPE emNode + CHARACTER(:), ALLOCATABLE:: type REAL(8):: phi REAL(8):: B(1:3) @@ -147,8 +148,8 @@ MODULE moduleOutput formatValues%density = rawValues%den*tempVol formatValues%velocity(:) = tempVel IF (tensorTrace(tensorTemp) > 0.D0) THEN - formatValues%pressure = speciesIn%m*tensorTrace(tensorTemp)*tempVol/3.D0 - formatValues%temperature = formatValues%pressure/(formatValues%density*kb) + formatValues%pressure = speciesIn%mass*tensorTrace(tensorTemp)*tempVol/3.D0 + formatValues%temperature = formatValues%pressure/(formatValues%density*kb) END IF END IF @@ -159,12 +160,12 @@ MODULE moduleOutput END SUBROUTINE calculateOutput - SUBROUTINE printTime(first) + SUBROUTINE printTime(t, first) USE moduleSpecies USE moduleCompTime - USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE + INTEGER, INTENT(in):: t LOGICAL, INTENT(in), OPTIONAL:: first CHARACTER(:), ALLOCATABLE:: fileName @@ -186,7 +187,7 @@ MODULE moduleOutput OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - WRITE (20, "(I10, I10, 7(ES20.6E3))") timeStep, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField + WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nParticlesTotal, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField CLOSE(20) diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index 7f6891c..bdf6b03 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -1,202 +1,55 @@ !Module to solve the electromagnetic field MODULE moduleEM - USE moduleMesh - USE moduleTable IMPLICIT NONE - ! Generic type for electromagnetic boundary conditions - TYPE, PUBLIC, ABSTRACT:: boundaryEMGeneric - INTEGER:: nNodes - TYPE(meshNodePointer), ALLOCATABLE:: nodes(:) - - CONTAINS - PROCEDURE(applyEM_interface), DEFERRED, PASS:: apply - - END TYPE boundaryEMGeneric - - ABSTRACT INTERFACE - ! Apply boundary condition to the load vector for the Poission equation - SUBROUTINE applyEM_interface(self, vectorF) - IMPORT boundaryEMGeneric - CLASS(boundaryEMGeneric), INTENT(in):: self - REAL(8), INTENT(inout):: vectorF(:) - - END SUBROUTINE applyEM_interface - - END INTERFACE - - TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichlet + TYPE:: boundaryEM + CHARACTER(:), ALLOCATABLE:: typeEM + INTEGER:: physicalSurface REAL(8):: potential CONTAINS - ! boundaryEMGeneric DEFERRED PROCEDURES - PROCEDURE, PASS:: apply => applyDirichlet + PROCEDURE, PASS:: apply - END TYPE boundaryEMDirichlet - - TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime - REAL(8):: potential - TYPE(table1D):: temporalProfile - - CONTAINS - ! boundaryEMGeneric DEFERRED PROCEDURES - PROCEDURE, PASS:: apply => applyDirichletTime - - END TYPE boundaryEMDirichletTime - - ! Container for boundary conditions - TYPE:: boundaryEMCont - CLASS(boundaryEMGeneric), ALLOCATABLE:: obj - - END TYPE boundaryEMCont + END TYPE boundaryEM INTEGER:: nBoundaryEM - TYPE(boundaryEMCont), ALLOCATABLE:: boundaryEM(:) + TYPE(boundaryEM), ALLOCATABLE:: boundEM(:) !Information of charge and reference parameters for rho vector REAL(8), ALLOCATABLE:: qSpecies(:) CONTAINS - SUBROUTINE findNodes(self, physicalSurface) + !Apply boundary conditions to the K matrix for Poisson's equation + SUBROUTINE apply(self, edge) USE moduleMesh IMPLICIT NONE - CLASS(boundaryEMGeneric), INTENT(inout):: self - INTEGER, INTENT(in):: physicalSurface - CLASS(meshEdge), POINTER:: edge - INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:) - INTEGER:: nNodes, nodesNew - INTEGER:: e, n + CLASS(boundaryEM), INTENT(in):: self + CLASS(meshEdge):: edge + INTEGER:: nNodes + INTEGER, ALLOCATABLE:: nodes(:) + INTEGER:: n - !Temporal array to hold nodes - ALLOCATE(nodes(0)) + nNodes = edge%nNodes + nodes = edge%getNodes(nNodes) - ! Loop thorugh the edges and identify those that are part of the boundary - DO e = 1, mesh%numEdges - edge => mesh%edges(e)%obj - IF (edge%physicalSurface == physicalSurface) THEN - ! Edge is of the right boundary index - ! Get nodes in the edge - nNodes = edge%nNodes - nodesEdge = edge%getNodes(nNodes) - ! Collect all nodes that are not already in the temporal array - DO n = 1, nNodes - IF (ANY(nodes == nodesEdge(n))) THEN - ! Node already in array, skip - CYCLE - - ELSE - ! If not, add element to array of nodes - nodes = [nodes, nodesEdge(n)] - - END IF - - END DO - - END IF - - END DO - - ! Point boundary to nodes - nNodes = SIZE(nodes) - ALLOCATE(self%nodes(nNodes)) - self%nNodes = nNodes DO n = 1, nNodes - self%nodes(n)%obj => mesh%nodes(nodes(n))%obj + SELECT CASE(self%typeEM) + CASE ("dirichlet") + mesh%K(nodes(n), :) = 0.D0 + mesh%K(nodes(n), nodes(n)) = 1.D0 + + mesh%nodes(nodes(n))%obj%emData%type = self%typeEM + mesh%nodes(nodes(n))%obj%emData%phi = self%potential + + END SELECT END DO - END SUBROUTINE findNodes - - ! Initialize Dirichlet boundary condition - SUBROUTINE initDirichlet(self, physicalSurface, potential) - USE moduleRefParam, ONLY: Volt_ref - IMPLICIT NONE - - CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self - INTEGER, INTENT(in):: physicalSurface - REAL(8), INTENT(in):: potential - - ! Allocate boundary edge - ALLOCATE(boundaryEMDirichlet:: self) - - SELECT TYPE(self) - TYPE IS(boundaryEMDirichlet) - self%potential = potential / Volt_ref - - CALL findNodes(self, physicalSurface) - - END SELECT - - END SUBROUTINE initDirichlet - - ! Initialize Dirichlet boundary condition - SUBROUTINE initDirichletTime(self, physicalSurface, potential, temporalProfile) - USE moduleRefParam, ONLY: Volt_ref, ti_ref - IMPLICIT NONE - - CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self - INTEGER, INTENT(in):: physicalSurface - REAL(8), INTENT(in):: potential - CHARACTER(:), ALLOCATABLE, INTENT(in):: temporalProfile - - ! Allocate boundary edge - ALLOCATE(boundaryEMDirichletTime:: self) - - SELECT TYPE(self) - TYPE IS(boundaryEMDirichletTime) - self%potential = potential / Volt_ref - - CALL findNodes(self, physicalSurface) - - CALL self%temporalProfile%init(temporalProfile) - - CALL self%temporalProfile%convert(1.D0/ti_ref, 1.D0) - - END SELECT - - END SUBROUTINE initDirichletTime - - !Apply Dirichlet boundary condition to the poisson equation - SUBROUTINE applyDirichlet(self, vectorF) - USE moduleMesh - IMPLICIT NONE - - CLASS(boundaryEMDirichlet), INTENT(in):: self - REAL(8), INTENT(inout):: vectorF(:) - INTEGER:: n, ni - - DO n = 1, self%nNodes - self%nodes(n)%obj%emData%phi = self%potential - vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi - - END DO - - END SUBROUTINE applyDirichlet - - !Apply Dirichlet boundary condition with time temporal profile - SUBROUTINE applyDirichletTime(self, vectorF) - USE moduleMesh - USE moduleCaseParam, ONLY: timeStep, tauMin - IMPLICIT NONE - - CLASS(boundaryEMDirichletTime), INTENT(in):: self - REAL(8), INTENT(inout):: vectorF(:) - REAL(8):: timeFactor - INTEGER:: n, ni - - timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin) - - DO n = 1, self%nNodes - self%nodes(n)%obj%emData%phi = self%potential * timeFactor - vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi - - END DO - - END SUBROUTINE applyDirichletTime + END SUBROUTINE !Assemble the source vector based on the charge density to solve Poisson's equation - SUBROUTINE assembleSourceVector(vectorF, n_e) + SUBROUTINE assembleSourceVector(vectorF) USE moduleMesh USE moduleRefParam IMPLICIT NONE @@ -205,9 +58,8 @@ MODULE moduleEM REAL(8), ALLOCATABLE:: localF(:) INTEGER, ALLOCATABLE:: nodes(:) REAL(8), ALLOCATABLE:: rho(:) - REAL(8), INTENT(in), OPTIONAL:: n_e(1:mesh%numNodes) INTEGER:: nNodes - INTEGER:: e, i, ni, b + INTEGER:: e, i, ni CLASS(meshNode), POINTER:: node !$OMP SINGLE @@ -225,10 +77,6 @@ MODULE moduleEM ni = nodes(i) node => mesh%nodes(ni)%obj rho(i) = DOT_PRODUCT(qSpecies(:), node%output(:)%den/(vol_ref*node%v*n_ref)) - IF (PRESENT(n_e)) THEN - rho(i) = rho(i) - n_e(i) - - END IF END DO @@ -249,12 +97,18 @@ MODULE moduleEM !$OMP END DO !Apply boundary conditions - !$OMP SINGLE - do b = 1, nBoundaryEM - call boundaryEM(b)%obj%apply(vectorF) + !$OMP DO + DO i = 1, mesh%numNodes + node => mesh%nodes(i)%obj - end do - !$OMP END SINGLE + SELECT CASE(node%emData%type) + CASE ("dirichlet") + vectorF(i) = node%emData%phi + + END SELECT + + END DO + !$OMP END DO END SUBROUTINE assembleSourceVector @@ -302,86 +156,4 @@ MODULE moduleEM END SUBROUTINE solveElecField - FUNCTION BoltzmannElectron(phi, n) RESULT(n_e) - USE moduleRefParam - USE moduleConstParam - IMPLICIT NONE - - INTEGER, INTENT(in):: n - REAL(8), INTENT(in):: phi(1:n) - REAL(8):: n_e(1:n) - REAL(8):: n_e0 = 1.0D16, phi_0 = -500.0D0, T_e = 11604.0 - INTEGER:: i - - n_e = n_e0 / n_ref * exp(qe * (phi*Volt_ref - phi_0) / (kb * T_e)) - - RETURN - - END FUNCTION BoltzmannElectron - - SUBROUTINE solveElecFieldBoltzmann() - USE moduleMesh - USE moduleErrors - IMPLICIT NONE - - INTEGER, SAVE:: INFO - INTEGER:: n - REAL(8), ALLOCATABLE, SAVE:: tempF(:) - REAL(8), ALLOCATABLE, SAVE:: n_e(:), phi_old(:), phi(:) - INTEGER:: k - EXTERNAL:: dgetrs - - !$OMP SINGLE - ALLOCATE(tempF(1:mesh%numNodes)) - ALLOCATE(n_e(1:mesh%numNodes)) - ALLOCATE(phi_old(1:mesh%numNodes)) - ALLOCATE(phi(1:mesh%numNodes)) - !$OMP END SINGLE - - !$OMP DO - DO n = 1, mesh%numNodes - phi_old(n) = mesh%nodes(n)%obj%emData%phi - - END DO - !$OMP END DO - - !$OMP SINGLE - DO k = 1, 100 - n_e = BoltzmannElectron(phi_old, mesh%numNodes) - CALL assembleSourceVector(tempF, n_e) - - CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, & - mesh%IPIV, tempF, mesh%numNodes, info) - phi = tempF - - PRINT *, MAXVAL(n_e), MINVAL(n_e) - PRINT *, MAXVAL(phi), MINVAL(phi) - PRINT*, k, "diff = ", MAXVAL(ABS(phi - phi_old)) - phi_old = phi - - END DO - !$OMP END SINGLE - - IF (info == 0) THEN - !Suscessful resolution of Poission equation - !$OMP DO - DO n = 1, mesh%numNodes - mesh%nodes(n)%obj%emData%phi = phi_old(n) - - END DO - !$OMP END DO - - ELSE - !$OMP SINGLE - CALL criticalError('Poisson equation failed', 'solveElecFieldBoltzmann') - !$OMP END SINGLE - - END IF - - !$OMP SINGLE - DEALLOCATE(tempF, n_e, phi_old, phi) - !$OMP END SINGLE - - END SUBROUTINE solveElecFieldBoltzmann - END MODULE moduleEM diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index df8d5e8..96e3f78 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -3,7 +3,7 @@ MODULE moduleSolver !Generic type for pusher of particles TYPE, PUBLIC:: pusherGeneric - PROCEDURE(push_interafece), POINTER, NOPASS:: pushParticle => NULL() + PROCEDURE(push_interface), POINTER, NOPASS:: pushParticle => NULL() !Boolean to indicate if the species is moved in the iteration LOGICAL:: pushSpecies !How many interations between advancing the species @@ -29,13 +29,13 @@ MODULE moduleSolver INTERFACE !Push a particle - PURE SUBROUTINE push_interafece(part, tauIn) + PURE SUBROUTINE push_interface(part, tauIn) USE moduleSpecies TYPE(particle), INTENT(inout):: part REAL(8), INTENT(in):: tauIn - END SUBROUTINE push_interafece + END SUBROUTINE push_interface !Solve the electromagnetic field SUBROUTINE solveEM_interface() @@ -123,7 +123,12 @@ MODULE moduleSolver END SELECT self%pushSpecies = .FALSE. - self%every = INT(tauSp/tau) + self%every = FLOOR(tauSp/tau) + !Correct value if not fulfilled + IF (tau*REAL(self%every) < tauSp) THEN + self%every = self%every + 1 + + END IF END SUBROUTINE initPusher @@ -138,9 +143,6 @@ MODULE moduleSolver CASE('Electrostatic','ConstantB') self%solveEM => solveElecField - CASE('ElectrostaticBoltzmann') - self%solveEM => solveElecFieldBoltzmann - END SELECT END SUBROUTINE initEM @@ -167,24 +169,75 @@ MODULE moduleSolver USE moduleMesh IMPLICIT NONE - INTEGER:: n - INTEGER:: sp + INTEGER:: p + INTEGER:: s, sp + INTEGER:: e - !$OMP DO - DO n = 1, nPartOld - !Select species type - sp = partOld(n)%species%n + !$OMP SINGLE + !Allocate the array of particles to push + nSpeciesToPush = COUNT(solver%pusher(:)%pushSpecies) + ALLOCATE(particlesToPush(1:nSpeciesToPush)) + ALLOCATE(nPartOldToPush(1:nSpeciesToPush)) + !Point the arrays to the location of the particles + sp = 0 + DO s = 1, nSpecies !Checks if the species sp is update this iteration IF (solver%pusher(sp)%pushSpecies) THEN - !Push particle - CALL solver%pusher(sp)%pushParticle(partOld(n), tau(sp)) - !Find cell in wich particle reside - CALL solver%updateParticleCell(partOld(n)) + sp = sp + 1 + particlesToPush(sp)%p = partOld(s)%p + nPartOldToPush(sp) = nPartOld(s) END IF END DO - !$OMP END DO + + !Delete list of particles in cells for collisions if the particle is pushed + IF (listInCells) THEN + DO e = 1, mesh%numCells + DO s = 1, nSpecies + IF (solver%pusher(s)%pushSpecies) THEN + CALL mesh%cells(e)%obj%listPart_in(s)%erase() + mesh%cells(e)%obj%totalWeight(s) = 0.D0 + + END IF + + END DO + + END DO + + END IF + + !Erase the list of particles inside the cell in coll mesh if the particle is pushed + IF (doubleMesh) THEN + DO e = 1, meshColl%numCells + DO s = 1, nSpecies + IF (solver%pusher(s)%pushSpecies) THEN + CALL meshColl%cells(e)%obj%listPart_in(s)%erase() + meshColl%cells(e)%obj%totalWeight(s) = 0.D0 + + END IF + + END DO + + END DO + + END IF + !$OMP END SINGLE + + !Now, push particles + !$OMP DO + DO sp = 1, nSpeciesToPush + DO p = 1, nPartOldToPush(sp) + !Push particle + CALL solver%pusher(sp)%pushParticle(particlesToPush(sp)%p(p), tau(sp)) + !Find cell in which particle reside + CALL solver%updateParticleCell(particlesToPush(sp)%p(p)) + + END DO + + END DO + !$END OMP DO + END SUBROUTINE doPushes @@ -245,7 +298,10 @@ MODULE moduleSolver !$OMP SECTION nOldIn = 0 IF (ALLOCATED(partOld)) THEN - nOldIn = COUNT(partOld%n_in) + DO s = 1, nSpecies + nOldIn = nOldin + COUNT(partOld(s)%p(:)%n_in) + + END DO END IF !$OMP SECTION @@ -322,40 +378,6 @@ MODULE moduleSolver END DO - !$OMP SECTION - !Erase the list of particles inside the cell if particles have been pushed - IF (listInCells) THEN - DO s = 1, nSpecies - DO e = 1, mesh%numCells - IF (solver%pusher(s)%pushSpecies) THEN - CALL mesh%cells(e)%obj%listPart_in(s)%erase() - mesh%cells(e)%obj%totalWeight(s) = 0.D0 - - END IF - - END DO - - END DO - - END IF - - !$OMP SECTION - !Erase the list of particles inside the cell in coll mesh - IF (doubleMesh) THEN - DO s = 1, nSpecies - DO e = 1, meshColl%numCells - IF (solver%pusher(s)%pushSpecies) THEN - CALL meshColl%cells(e)%obj%listPart_in(s)%erase() - meshColl%cells(e)%obj%totalWeight(s) = 0.D0 - - END IF - - END DO - - END DO - - END IF - !$OMP END SECTIONS !$OMP SINGLE @@ -370,14 +392,17 @@ MODULE moduleSolver USE moduleMesh IMPLICIT NONE - INTEGER:: n + INTEGER:: n, sp CLASS(meshCell), POINTER:: cell !Loops over the particles to scatter them !$OMP DO - DO n = 1, nPartOld - cell => mesh%cells(partOld(n)%cell)%obj - CALL cell%scatter(cell%nNodes, partOld(n)) + DO sp = 1, nSpeciesToPush + DO n = 1, nPartOldToPush(sp) + cell => mesh%cells(particlesToPush(sp)%p(n)%cell)%obj + CALL cell%scatter(cell%nNodes, particlesToPush(sp)%p(n)) + + END DO END DO !$OMP END DO @@ -494,46 +519,52 @@ MODULE moduleSolver END SUBROUTINE updateParticleCell !Update the information about if a species needs to be moved this iteration - SUBROUTINE updatePushSpecies(self) + SUBROUTINE updatePushSpecies(self, t) USE moduleSpecies - USE moduleCaseparam, ONLY: timeStep IMPLICIT NONE CLASS(solverGeneric), INTENT(inout):: self + INTEGER, INTENT(in):: t INTEGER:: s DO s=1, nSpecies - self%pusher(s)%pushSpecies = MOD(timeStep, self%pusher(s)%every) == 0 + self%pusher(s)%pushSpecies = MOD(t, self%pusher(s)%every) == 0 END DO END SUBROUTINE updatePushSpecies !Output the different data and information - SUBROUTINE doOutput() + SUBROUTINE doOutput(t) USE moduleMesh USE moduleOutput USE moduleSpecies USE moduleCompTime USE moduleProbe - USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - CALL outputProbes() + INTEGER, INTENT(in):: t + + IF (t == tInitial) THEN + CALL SYSTEM('git rev-parse HEAD > ' // path // folder // '/' // 'fpack_commit.txt') + + END IF + + CALL outputProbes(t) counterOutput = counterOutput + 1 IF (counterOutput >= triggerOutput .OR. & - timeStep == tFinal .OR. timeStep == tInitial) THEN + t == tFinal .OR. t == tInitial) THEN !Resets output counter counterOutput=0 - CALL mesh%printOutput() - IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl() - CALL mesh%printEM() - WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal + CALL mesh%printOutput(t) + IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t) + CALL mesh%printEM(t) + WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld - IF (timeStep == 0) THEN + IF (t == 0) THEN WRITE(*, "(5X,A21,F8.1,A2)") " init time: ", 1.D3*tStep, "ms" ELSE @@ -541,8 +572,8 @@ MODULE moduleSolver END IF - IF (nPartOld > 0) THEN - WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOld), "ns" + IF (nPartOldTotal > 0) THEN + WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOldTotal), "ns" END IF WRITE(*,*) @@ -551,32 +582,34 @@ MODULE moduleSolver counterCPUTime = counterCPUTime + 1 IF (counterCPUTime >= triggerCPUTime .OR. & - timeStep == tFinal .OR. timeStep == tInitial) THEN + t == tFinal .OR. t == tInitial) THEN !Reset CPU Time counter counterCPUTime = 0 - CALL printTime(timeStep == 0) + CALL printTime(t, t == 0) END IF !Output average values - IF (useAverage .AND. timeStep == tFinal) THEN + IF (useAverage .AND. t == tFinal) THEN CALL mesh%printAverage() END IF END SUBROUTINE doOutput - SUBROUTINE doAverage() + SUBROUTINE doAverage(t) USE moduleAverage USE moduleMesh IMPLICIT NONE + INTEGER, INTENT(in):: t INTEGER:: tAverage, n + IF (useAverage) THEN - tAverage = timeStep - tAverageStart + tAverage = t - tAverageStart IF (tAverage == 1) THEN !First iteration in which average scheme is used diff --git a/src/modules/solver/pusher/modulePusher.f90 b/src/modules/solver/pusher/modulePusher.f90 index f75c733..4c7e372 100644 --- a/src/modules/solver/pusher/modulePusher.f90 +++ b/src/modules/solver/pusher/modulePusher.f90 @@ -14,7 +14,7 @@ MODULE modulePusher END SUBROUTINE pushCartNeutral PURE SUBROUTINE pushCartElectrostatic(part, tauIn) - USE moduleSPecies + USE moduleSpecies USE moduleMesh IMPLICIT NONE