diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 0000000..18cd6d7 --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,50 @@ +# 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 new file mode 100644 index 0000000..532ca4d --- /dev/null +++ b/data/collisions/IO_e-Kr.dat @@ -0,0 +1,52 @@ +# 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 new file mode 100644 index 0000000..3067578 --- /dev/null +++ b/data/collisions/IO_e-Xe.dat @@ -0,0 +1,52 @@ +# 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/doc/logos/fpakc.pdf b/doc/logos/fpakc.pdf new file mode 100644 index 0000000..36bc5ab Binary files /dev/null and b/doc/logos/fpakc.pdf differ diff --git a/doc/user-manual/figures/logos/PPartiC.eps b/doc/user-manual/figures/logos/PPartiC.eps deleted file mode 100644 index 5517180..0000000 --- a/doc/user-manual/figures/logos/PPartiC.eps +++ /dev/null @@ -1,177 +0,0 @@ -%!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 deleted file mode 100644 index 9966134..0000000 --- a/doc/user-manual/figures/logos/PPartiC.svg +++ /dev/null @@ -1,159 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - image/svg+xml - - - - - - - PPartiC - - - - - - diff --git a/doc/user-manual/figures/logos/fpakc.eps b/doc/user-manual/figures/logos/fpakc.eps deleted file mode 100644 index 667a7f0..0000000 --- a/doc/user-manual/figures/logos/fpakc.eps +++ /dev/null @@ -1,294 +0,0 @@ -%!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 deleted file mode 100644 index ca8fd14..0000000 --- a/doc/user-manual/figures/scatteringQuad.eps +++ /dev/null @@ -1,303 +0,0 @@ -%!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 ef1494a..5321ba4 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 13be4cd..845cc9a 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} @@ -460,6 +460,10 @@ 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. @@ -585,12 +589,20 @@ make Type of boundary. Accepted values are: \begin{itemize} - \item \textbf{dirichlet}: Elastic reflection of particles. + \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}. \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} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -611,7 +623,7 @@ make \begin{itemize} \item \textbf{A}: Ampere. \item \textbf{Am2}: Ampere per square meter. - This value will be multiplied by the surface of injection. + This value will be multiplied by the area of injection. \item \textbf{sccm}: Standard cubic centimetre. \item \textbf{part/s}: Particles (real) per second. \end{itemize} @@ -717,7 +729,7 @@ make 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}. + File must be located in \textbf{output.path}. \item \textbf{particlesPerCell}: Integer. Optional. Initial number of particles per cell. diff --git a/src/fpakc.f90 b/src/fpakc.f90 index e90a5e0..ae6e7bb 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 @@ -32,10 +32,13 @@ PROGRAM fpakc CALL initOutput(inputFile) !Do '0' iteration - t = tInitial + timeStep = tInitial !$OMP PARALLEL DEFAULT(SHARED) !$OMP SINGLE + ! Initial reset of probes + CALL resetProbes() + CALL verboseError("Initial scatter of particles...") !$OMP END SINGLE CALL doScatter() @@ -49,19 +52,21 @@ PROGRAM fpakc tStep = omp_get_wtime() - tStep !Output initial state - CALL doOutput(t) + CALL doOutput() 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(t) + CALL solver%updatePushSpecies() !Checks if probes need to be calculated this iteration - CALL resetProbes(t) + CALL resetProbes() tPush = omp_get_wtime() !$OMP END SINGLE @@ -79,7 +84,7 @@ PROGRAM fpakc !$OMP END SINGLE IF (doMCCollisions) THEN - CALL meshForMCC%doCollisions(t) + CALL meshForMCC%doCollisions() END IF @@ -124,12 +129,12 @@ PROGRAM fpakc !$OMP SINGLE tEMField = omp_get_wtime() - tEMField - CALL doAverage(t) + CALL doAverage() tStep = omp_get_wtime() - tStep !Output data - CALL doOutput(t) + CALL doOutput() !$OMP END SINGLE END DO diff --git a/src/makefile b/src/makefile index 247c7d2..83cf8db 100644 --- a/src/makefile +++ b/src/makefile @@ -9,6 +9,7 @@ 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 8df3210..551d867 100644 --- a/src/modules/common/moduleCaseParam.f90 +++ b/src/modules/common/moduleCaseParam.f90 @@ -2,8 +2,13 @@ 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 ae5c548..1b8ec67 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -47,24 +47,41 @@ MODULE moduleRandom 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 randomMaxwellian() RESULT(rnd) - USE moduleConstParam, ONLY: PI + FUNCTION randomHalfMaxwellian() RESULT(rnd) IMPLICIT NONE REAL(8):: rnd - REAL(8):: x, y + REAL(8):: x rnd = 0.D0 x = 0.D0 DO WHILE (x == 0.D0) CALL RANDOM_NUMBER(x) END DO - CALL RANDOM_NUMBER(y) - rnd = DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + rnd = DSQRT(-DLOG(x)) - END FUNCTION randomMaxwellian + END FUNCTION randomHalfMaxwellian !Returns a random number weighted with the cumWeight array FUNCTION randomWeighted(cumWeight, sumWeight) RESULT(rnd) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 6f1d5bb..9f50947 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -259,13 +259,17 @@ 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') @@ -799,7 +803,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 @@ -815,8 +819,8 @@ MODULE moduleInput 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) @@ -825,71 +829,77 @@ 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 - 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) + 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) - CASE('absorption') - ALLOCATE(boundaryAbsorption:: boundary(i)%bTypes(s)%obj) + CASE('absorption') + ALLOCATE(boundaryAbsorption:: bound) - CASE('transparent') - ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj) + CASE('transparent') + ALLOCATE(boundaryTransparent:: 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('axis') + ALLOCATE(boundaryAxis:: bound) - CALL config%get(object // '.effectiveTime', effTime, found) - IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for 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 // '.energyThreshold', eThreshold, found) - IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary') + CALL initWallTemperature(bound, Tw, cw) - CALL config%get(object // '.crossSection', crossSection, found) - IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals 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 // '.electronSecondary', electronSecondary, found) - electronSecondaryID = speciesName2Index(electronSecondary) - IF (found) THEN - CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, & - speciesID, effTime, crossSection, eThreshold,electronSecondaryID) + CALL config%get(object // '.effectiveTime', effTime, found) + IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary') - ELSE - 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') - END IF + CALL config%get(object // '.crossSection', crossSection, found) + IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' 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 // '.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) - CALL initWallTemperature(boundary(i)%bTypes(s)%obj, Tw, cw) + ELSE + CALL initIonization(bound, species(s)%obj%m, m0, n0, v0, T0, & + speciesID, effTime, crossSection, eThreshold) - CASE('axis') - ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj) + END IF - CASE DEFAULT - CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary') + case('outflowAdaptive') + allocate(boundaryOutflowAdaptive:: bound) - END SELECT + CASE DEFAULT + CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary') + + END SELECT + + end associate END DO @@ -906,6 +916,7 @@ MODULE moduleInput USE moduleMeshInputGmsh2, ONLY: initGmsh2 USE moduleMeshInputVTU, ONLY: initVTU USE moduleMeshInput0D, ONLY: init0D + USE moduleMeshInputText, ONLY: initText USE moduleMesh3DCart USE moduleMesh2DCyl USE moduleMesh2DCart @@ -964,9 +975,9 @@ MODULE moduleInput !Read the 0D mesh CALL mesh%readMesh(pathMeshParticle) - !Get the volumne + !Get the volume CALL config%get(object // '.volume', volume, found) - !Rescale the volumne + !Rescale the volume 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 @@ -1054,6 +1065,20 @@ 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') @@ -1100,13 +1125,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):: timeStep + REAL(8):: everyTimeStep CALL config%info('output.probes', found, n_children = nProbes) @@ -1114,7 +1139,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) @@ -1122,16 +1147,14 @@ 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', timeStep, found) + CALL config%get(object // '.timeStep', everyTimeStep, 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, timeStep) + CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, everyTimeStep) END DO - CALL resetProbes(tInitial) - END SUBROUTINE readProbes SUBROUTINE readEMBoundary(config) @@ -1139,7 +1162,6 @@ MODULE moduleInput USE moduleOutput USE moduleErrors USE moduleEM - USE moduleRefParam USE moduleSpecies USE json_module IMPLICIT NONE @@ -1147,34 +1169,72 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config CHARACTER(:), ALLOCATABLE:: object LOGICAL:: found - CHARACTER(2):: istring - INTEGER:: i, e, s + CHARACTER(:), ALLOCATABLE:: typeEM + REAL(8):: potential + INTEGER:: physicalSurface + CHARACTER(:), ALLOCATABLE:: temporalProfile, temporalProfilePath + INTEGER:: b, s, n, ni + CHARACTER(2):: bString INTEGER:: info EXTERNAL:: dgetrf CALL config%info('boundaryEM', found, n_children = nBoundaryEM) - IF (found) ALLOCATE(boundEM(1:nBoundaryEM)) + IF (found) THEN + ALLOCATE(boundaryEM(1:nBoundaryEM)) + + END IF - DO i = 1, nBoundaryEM - WRITE(istring, '(I2)') i - object = 'boundaryEM(' // trim(istring) // ')' + DO b = 1, nBoundaryEM + WRITE(bString, '(I2)') b + object = 'boundaryEM(' // TRIM(bString) // ')' - CALL config%get(object // '.type', boundEM(i)%typeEM, found) + CALL config%get(object // '.type', typeEM, found) - SELECT CASE(boundEM(i)%typeEM) + SELECT CASE(typeEM) CASE ("dirichlet") - CALL config%get(object // '.potential', boundEM(i)%potential, found) - IF (.NOT. found) & + CALL config%get(object // '.potential', potential, found) + IF (.NOT. found) THEN CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary') - boundEM(i)%potential = boundEM(i)%potential/Volt_ref - CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found) - IF (.NOT. found) & - CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary') + 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) CASE DEFAULT - CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary') + CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary') END SELECT @@ -1193,18 +1253,28 @@ MODULE moduleInput END DO - 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) + ! 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 - END IF - END DO - END IF - END DO - END IF + 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 !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) @@ -1225,13 +1295,13 @@ 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:: T(:), normal(:) + REAL(8), ALLOCATABLE:: temperature(:), normal(:) REAL(8):: flow CHARACTER(:), ALLOCATABLE:: units INTEGER:: physicalSurface @@ -1242,8 +1312,8 @@ MODULE moduleInput 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) @@ -1251,7 +1321,7 @@ MODULE moduleInput CALL config%get(object // '.name', name, found) CALL config%get(object // '.v', v, found) - CALL config%get(object // '.T', T, found) + CALL config%get(object // '.T', temperature, found) CALL config%get(object // '.n', normal, found) IF (.NOT. found) THEN ALLOCATE(normal(1:3)) @@ -1263,7 +1333,7 @@ MODULE moduleInput particlesPerEdge = 0 CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found) - CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface, particlesPerEdge) + CALL inject(i)%init(i, v, normal, temperature, flow, units, sp, physicalSurface, particlesPerEdge) CALL readVelDistr(config, inject(i), object) @@ -1282,6 +1352,7 @@ 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 @@ -1295,8 +1366,11 @@ MODULE moduleInput CALL config%get('average.startTime', tStart, found) IF (found) THEN - tAverageStart = INT(tStart / tauMin) + tAverageStart = INT(tStart / ti_ref / tauMin) + ELSE + tAverageStart = 0 + END IF ALLOCATE(averageScheme(1:mesh%numNodes)) @@ -1323,28 +1397,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, T, m + REAL(8):: v, temperature, 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") - T = inj%T(i) - CALL initVelDistMaxwellian(inj%v(i)%obj, t, m) + temperature = inj%temperature(i) + CALL initVelDistMaxwellian(inj%v(i)%obj, temperature, m) CASE ("Half-Maxwellian") - T = inj%T(i) - CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m) + temperature = inj%temperature(i) + CALL initVelDistHalfMaxwellian(inj%v(i)%obj, temperature, m) CASE ("Delta") v = inj%vMod*inj%n(i) diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index f400ab0..a304c5e 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -104,7 +104,6 @@ MODULE moduleMesh1DCart USE moduleSpecies USE moduleBoundary USE moduleErrors - USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge1DCart), INTENT(out):: self @@ -123,7 +122,7 @@ MODULE moduleMesh1DCart self%x = r1(1) - self%surface = 1.D0 / L_ref**2 + self%surface = 1.D0 self%normal = (/ 1.D0, 0.D0, 0.D0 /) diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index fd617bd..a46c6a2 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -104,7 +104,6 @@ MODULE moduleMesh1DRad USE moduleSpecies USE moduleBoundary USE moduleErrors - USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge1DRad), INTENT(out):: self @@ -123,7 +122,7 @@ MODULE moduleMesh1DRad self%r = r1(1) - self%surface = 1.D0 / L_ref**2 + self%surface = 1.D0 self%normal = (/ 1.D0, 0.D0, 0.D0 /) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index dbc8b25..db1ef4f 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -144,7 +144,6 @@ MODULE moduleMesh2DCart USE moduleSpecies USE moduleBoundary USE moduleErrors - USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge2DCart), INTENT(out):: self @@ -164,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) / L_ref + self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2) !Normal vector self%normal = (/ -(self%y(2)-self%y(1)), & self%x(2)-self%x(1) , & @@ -495,34 +494,36 @@ MODULE moduleMesh2DCart END FUNCTION insideQuad - !Transform physical coordinates to element coordinates + !Transform physical coordinates to element coordinates with a Taylor series 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):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) + REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2) 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 - XiO = 0.D0 + conv = 1.D0 + Xi0 = 0.D0 + Xi(3) = 0.D0 - f(3) = 0.D0 DO WHILE(conv > 1.D-4) - dPsi = self%dPsi(XiO, 4) - pDer = self%partialDer(4, dPsi) - detJ = self%detJac(pDer) - invJ = self%invJac(pDer) - fPsi = self%fPsi(XiO, 4) - f(1:2) = (/ DOT_PRODUCT(fPsi,self%x), & - DOT_PRODUCT(fPsi,self%y) /) - r(1:2) - Xi = XiO - MATMUL(invJ, f)/detJ - conv = MAXVAL(DABS(Xi-XiO),1) - XiO = Xi + 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) END DO @@ -557,7 +558,6 @@ MODULE moduleMesh2DCart !Compute element volume PURE SUBROUTINE volumeQuad(self) - USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(inout):: self @@ -575,7 +575,7 @@ MODULE moduleMesh2DCart fPsi = self%fPsi(Xi, 4) !Compute total volume of the cell - self%volume = detJ*4.D0/L_ref + self%volume = detJ*4.D0 !Compute volume per node self%n1%v = self%n1%v + fPsi(1)*self%volume self%n2%v = self%n2%v + fPsi(2)*self%volume @@ -680,8 +680,8 @@ MODULE moduleMesh2DCart dPsi = 0.D0 - dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /) - dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /) + dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /) + dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTria @@ -826,19 +826,19 @@ MODULE moduleMesh2DCart CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: deltaR(1:3) - REAL(8):: dPsi(1:3, 1:3) + REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2) + REAL(8):: dPsi(1:3,1:4) REAL(8):: pDer(1:3, 1:3) - REAL(8):: invJ(1:3, 1:3), detJ !Direct method to convert coordinates - 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 + 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 END FUNCTION phy2logTria @@ -913,8 +913,8 @@ MODULE moduleMesh2DCart invJ = 0.D0 - invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) - invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) + invJ(1, 1:2) = (/ pDer(2,2), -pDer(2,1) /) + invJ(2, 1:2) = (/ -pDer(1,2), 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 ae1eb92..1f5f33c 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -510,34 +510,36 @@ MODULE moduleMesh2DCyl END FUNCTION insideQuad - !Transform physical coordinates to element coordinates + !Transform physical coordinates to element coordinates with a Taylor series 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):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) + REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2) 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 - XiO = 0.D0 + conv = 1.D0 + Xi0 = 0.D0 + Xi(3) = 0.D0 - f(3) = 0.D0 DO WHILE(conv > 1.D-4) - dPsi = self%dPsi(XiO, 4) - pDer = self%partialDer(4, dPsi) - detJ = self%detJac(pDer) - invJ = self%invJac(pDer) - fPsi = self%fPsi(XiO, 4) - f(1:2) = (/ DOT_PRODUCT(fPsi,self%z), & - DOT_PRODUCT(fPsi,self%r) /) - r(1:2) - Xi = XiO - MATMUL(invJ, f)/detJ - conv = MAXVAL(DABS(Xi-XiO),1) - XiO = Xi + 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) END DO @@ -705,8 +707,8 @@ MODULE moduleMesh2DCyl dPsi = 0.D0 - dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /) - dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /) + dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /) + dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTria @@ -858,19 +860,19 @@ MODULE moduleMesh2DCyl CLASS(meshCell2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: deltaR(1:3) - REAL(8):: dPsi(1:3, 1:3) + REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2) + REAL(8):: dPsi(1:3,1:4) REAL(8):: pDer(1:3, 1:3) - REAL(8):: invJ(1:3, 1:3), detJ !Direct method to convert coordinates - 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 + 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 END FUNCTION phy2logTria @@ -948,8 +950,8 @@ MODULE moduleMesh2DCyl invJ = 0.D0 - invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) - invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) + invJ(1, 1:2) = (/ pDer(2,2), -pDer(2,1) /) + invJ(2, 1:2) = (/ -pDer(1,2), pDer(1,1) /) invJ(3, 3) = 1.D0 END FUNCTION invJ2DCyl diff --git a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 index 97ec729..c0dcfbb 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, t) + SUBROUTINE printOutput0D(self) 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 (t == 0) THEN + IF (timeStep == 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,14 +27,17 @@ 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(t)*tauMin*ti_ref, output%density, output%velocity, output%pressure, output%temperature + WRITE(20, "(7(ES20.6E3))") REAL(timeStep)*tauMin*ti_ref, output%density, & + output%velocity, & + output%pressure, & + output%temperature CLOSE(20) END DO END SUBROUTINE printOutput0D - SUBROUTINE printColl0D(self, t) + SUBROUTINE printColl0D(self) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -43,12 +46,11 @@ MODULE moduleMeshOutput0D IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t CHARACTER(:), ALLOCATABLE:: fileName INTEGER:: k fileName='OUTPUT_Collisions.dat' - IF (t == tInitial) THEN + IF (timeStep == 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 @@ -57,12 +59,12 @@ MODULE moduleMeshOutput0D END IF OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs) + WRITE(20, "(ES20.6E3, 10I20)") REAL(timeStep)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs) CLOSE(20) END SUBROUTINE printColl0D - SUBROUTINE printEM0D(self, t) + SUBROUTINE printEM0D(self) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -70,7 +72,6 @@ MODULE moduleMeshOutput0D IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t END SUBROUTINE printEM0D diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index ccdcf3d..8176bb5 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, t) + SUBROUTINE printOutputGmsh2(self) 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(t)*tauMin*ti_ref + time = DBLE(timeStep)*tauMin*ti_ref DO i = 1, nSpecies - fileName = formatFileName(prefix, species(i)%obj%name, 'msh', t) + fileName = formatFileName(prefix, species(i)%obj%name, 'msh', timeStep) 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)', t, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', timeStep, 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)', t, time, 3, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', timeStep, 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)', t, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', timeStep, 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)', t, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', timeStep, 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, t) + SUBROUTINE printCollGmsh2(self) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -145,7 +145,6 @@ MODULE moduleMeshOutputGmsh2 IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: numEdges INTEGER:: k, c INTEGER:: n @@ -167,9 +166,9 @@ MODULE moduleMeshOutputGmsh2 END SELECT IF (collOutput) THEN - time = DBLE(t)*tauMin*ti_ref + time = DBLE(timeStep)*tauMin*ti_ref - fileName = formatFileName(prefix, 'Collisions', 'msh', t) + fileName = formatFileName(prefix, 'Collisions', 'msh', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (60, file = path // folder // '/' // fileName) @@ -179,7 +178,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, t, time, 1, self%numCells) + CALL writeGmsh2HeaderElementData(60, title, timeStep, 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 @@ -196,7 +195,7 @@ MODULE moduleMeshOutputGmsh2 END SUBROUTINE printCollGmsh2 !Prints the electrostatic EM properties into the nodes and volumes - SUBROUTINE printEMGmsh2(self, t) + SUBROUTINE printEMGmsh2(self) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -205,7 +204,6 @@ MODULE moduleMeshOutputGmsh2 IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: n, e REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName @@ -214,27 +212,27 @@ MODULE moduleMeshOutputGmsh2 Xi = (/ 0.D0, 0.D0, 0.D0 /) IF (emOutput) THEN - time = DBLE(t)*tauMin*ti_ref + time = DBLE(timeStep)*tauMin*ti_ref - fileName = formatFileName(prefix, 'EMField', 'msh', t) + fileName = formatFileName(prefix, 'EMField', 'msh', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (20, file = path // folder // '/' // fileName) CALL writeGmsh2HeaderMesh(20) - CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', t, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', timeStep, 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)', t, time, 3, self%numCells) + CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', timeStep, 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)', t, time, 3, self%numNodes) + CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', timeStep, 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 a93161d..1b8883d 100644 --- a/src/modules/mesh/inout/makefile +++ b/src/modules/mesh/inout/makefile @@ -1,4 +1,4 @@ -all: vtu.o gmsh2.o 0D.o +all: vtu.o gmsh2.o 0D.o text.o vtu.o: moduleMeshInoutCommon.o $(MAKE) -C vtu all @@ -9,5 +9,8 @@ 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 e4a6c72..7dc2e84 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, t) RESULT(fileName) + PURE FUNCTION formatFileName(prefix, suffix, extension, timeStep) RESULT(fileName) USE moduleOutput IMPLICIT NONE CHARACTER(*), INTENT(in):: prefix, suffix, extension - INTEGER, INTENT(in), OPTIONAL:: t + INTEGER, INTENT(in), OPTIONAL:: timeStep CHARACTER (LEN=iterationDigits):: tString CHARACTER(:), ALLOCATABLE:: fileName - IF (PRESENT(t)) THEN - WRITE(tString, iterationFormat) t + IF (PRESENT(timeStep)) THEN + WRITE(tString, iterationFormat) timeStep fileName = prefix // '_' // tString // '_' // suffix // '.' // extension ELSE diff --git a/src/modules/mesh/inout/text/makefile b/src/modules/mesh/inout/text/makefile new file mode 100644 index 0000000..26d0f10 --- /dev/null +++ b/src/modules/mesh/inout/text/makefile @@ -0,0 +1,7 @@ +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 new file mode 100644 index 0000000..a41e6d3 --- /dev/null +++ b/src/modules/mesh/inout/text/moduleMeshInputText.f90 @@ -0,0 +1,232 @@ +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 new file mode 100644 index 0000000..d69ddda --- /dev/null +++ b/src/modules/mesh/inout/text/moduleMeshOutputText.f90 @@ -0,0 +1,265 @@ +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 184d645..c7b89b5 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, error, found + INTEGER:: fileID CHARACTER(LEN=256):: line INTEGER:: numNodes, numElements, numEdges INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types @@ -548,6 +548,8 @@ 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 da90b6b..8e8f91a 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,ES20.6E3,A)") '' + WRITE(fileID,"(4X, A)") '' WRITE(fileID,"(6X, A, I10, A, I10, A)") '' END SUBROUTINE writeHeader @@ -215,17 +215,16 @@ MODULE moduleMeshOutputVTU END SUBROUTINE writeEM - SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection) + SUBROUTINE writeCollection(fileID, fileNameStep, fileNameCollection) USE moduleCaseParam USE moduleOutput USE moduleRefParam IMPLICIT NONE INTEGER:: fileID - INTEGER, INTENT(in):: t CHARACTER(*):: fileNameStep, fileNameCollection - IF (t == tInitial) THEN + IF (timeStep == tInitial) THEN !Create collection file WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection) @@ -237,10 +236,11 @@ 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 (t == tFinal) THEN + IF (timeStep == tFinal) THEN WRITE (fileID + 1, "(2X, A)") '' WRITE (fileID + 1, "(A)") '' @@ -307,21 +307,21 @@ MODULE moduleMeshOutputVTU END SUBROUTINE writeAverage - SUBROUTINE printOutputVTU(self,t) + SUBROUTINE printOutputVTU(self) USE moduleMesh USE moduleSpecies USE moduleMeshInoutCommon + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: i, fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection fileID = 60 DO i = 1, nSpecies - fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t) + fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) @@ -337,28 +337,27 @@ MODULE moduleMeshOutputVTU !Write collection file for time plotting fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd') - CALL writeCollection(fileID, t, fileName, filenameCollection) + CALL writeCollection(fileID, fileName, filenameCollection) END DO END SUBROUTINE printOutputVTU - SUBROUTINE printCollVTU(self,t) + SUBROUTINE printCollVTU(self) USE moduleMesh USE moduleOutput USE moduleMeshInoutCommon + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection - CHARACTER (LEN=iterationDigits):: tstring fileID = 62 IF (collOutput) THEN - fileName = formatFileName(prefix, 'Collisions', 'vtu', t) + fileName = formatFileName(prefix, 'Collisions', 'vtu', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) @@ -374,26 +373,26 @@ MODULE moduleMeshOutputVTU !Write collection file for time plotting fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd') - CALL writeCollection(fileID, t, fileName, filenameCollection) + CALL writeCollection(fileID, fileName, filenameCollection) END IF END SUBROUTINE printCollVTU - SUBROUTINE printEMVTU(self, t) + SUBROUTINE printEMVTU(self) 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', t) + fileName = formatFileName(prefix, 'EMField', 'vtu', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) @@ -409,7 +408,7 @@ MODULE moduleMeshOutputVTU !Write collection file for time plotting fileNameCollection = formatFileName('Collection', 'EMField', 'pvd') - CALL writeCollection(fileID, t, fileName, filenameCollection) + CALL writeCollection(fileID, fileName, filenameCollection) END IF diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index ae53aa0..01b03c7 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -59,6 +59,13 @@ 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() @@ -337,10 +344,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 @@ -372,10 +379,9 @@ MODULE moduleMesh END SUBROUTINE connectMesh_interface !Prints number of collisions in each cell - SUBROUTINE printColl_interface(self, t) + SUBROUTINE printColl_interface(self) IMPORT meshGeneric CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t END SUBROUTINE printColl_interface @@ -403,18 +409,16 @@ MODULE moduleMesh ABSTRACT INTERFACE !Prints Species data - SUBROUTINE printOutput_interface(self, t) + SUBROUTINE printOutput_interface(self) IMPORT meshParticles CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t END SUBROUTINE printOutput_interface !Prints EM info - SUBROUTINE printEM_interface(self, t) + SUBROUTINE printEM_interface(self) IMPORT meshParticles CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t END SUBROUTINE printEM_interface @@ -495,28 +499,29 @@ MODULE moduleMesh IMPLICIT NONE CLASS(meshParticles), INTENT(inout):: self - INTEGER:: e - INTEGER:: nNodes + INTEGER:: c INTEGER, ALLOCATABLE:: n(:) REAL(8), ALLOCATABLE:: localK(:,:) INTEGER:: i, j - 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 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 i = 1, nNodes - DO j = 1, nNodes - self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j) + 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 END DO - END DO + DEALLOCATE(n, localK) - DEALLOCATE(n, localK) + end associate END DO @@ -789,7 +794,7 @@ MODULE moduleMesh END FUNCTION findCellBrute !Computes collisions in element - SUBROUTINE doCollisions(self, t) + SUBROUTINE doCollisions(self) USE moduleCollisions USE moduleSpecies USE moduleList @@ -797,10 +802,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 @@ -816,7 +821,7 @@ MODULE moduleMesh REAL(8):: rnd_real !Random number for collision INTEGER:: rnd_int !Random number for collision - IF (MOD(t, everyColl) == 0) THEN + IF (MOD(timeStep, everyColl) == 0) THEN !Collisions need to be performed in this iteration !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j) DO e=1, self%numCells diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 091e52e..deb4459 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -77,6 +77,20 @@ 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 @@ -204,19 +218,27 @@ MODULE moduleMeshBoundary END SUBROUTINE ionization - !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 + subroutine outflowAdaptive(edge, part) + use moduleRandom + implicit none - CLASS(meshEdge), INTENT(inout):: edge - CLASS(particle), INTENT(inout):: part + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part - CALL reflection(edge, part) + select type(bound => edge%boundary%bTypes(part%species%n)%obj) + type is(boundaryOutflowAdaptive) - END SUBROUTINE symmetryAxis + if (random() < 0.844d0) then + call reflection(edge, part) + + else + call transparent(edge, part) + + end if + + end select + + end subroutine outflowAdaptive !Points the boundary function to specific type SUBROUTINE pointBoundaryFunction(edge, s) @@ -236,14 +258,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(boundaryAxis) - edge%fBoundary(s)%apply => symmetryAxis + type is(boundaryOutflowAdaptive) + edge%fBoundary(s)%apply => outflowAdaptive 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 0b76105..67465c7 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -26,6 +26,12 @@ 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) @@ -47,11 +53,13 @@ MODULE moduleBoundary END TYPE boundaryIonization - !Symmetry axis - TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis - CONTAINS + !Boundary for quasi-neutral outflow adjusting reflection coefficient + type, public, extends(boundaryGeneric):: boundaryOutflowAdaptive + real(8):: outflowCurrent + real(8):: reflectionFraction + contains - END TYPE boundaryAxis + end type boundaryOutflowAdaptive !Wrapper for boundary types (one per species) TYPE:: bTypesCont diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 18c7fbb..ff8a694 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -54,7 +54,7 @@ MODULE moduleInject INTEGER:: id CHARACTER(:), ALLOCATABLE:: name REAL(8):: vMod !Velocity (module) - REAL(8):: T(1:3) !Temperature + REAL(8):: temperature(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 @@ -76,7 +76,7 @@ MODULE moduleInject CONTAINS !Initialize an injection of particles - SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface, particlesPerEdge) + SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, physicalSurface, particlesPerEdge) USE moduleMesh USE moduleRefParam USE moduleConstParam @@ -87,7 +87,7 @@ MODULE moduleInject CLASS(injectGeneric), INTENT(inout):: self INTEGER, INTENT(in):: i - REAL(8), INTENT(in):: v, n(1:3), T(1:3) + REAL(8), INTENT(in):: v, n(1:3), temperature(1:3) INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge REAL(8):: tauInject REAL(8), INTENT(in):: flow @@ -97,10 +97,10 @@ MODULE moduleInject INTEGER:: nVolColl REAL(8):: fluxPerStep = 0.D0 - self%id = i - self%vMod = v / v_ref - self%n = n / NORM2(n) - self%T = T / T_ref + self%id = i + self%vMod = v / v_ref + self%n = n / NORM2(n) + self%temperature = temperature / T_ref !Gets the edge elements from which particles are injected DO e = 1, mesh%numEdges phSurface(e) = mesh%edges(e)%obj%physicalSurface @@ -159,11 +159,26 @@ MODULE moduleInject CASE ("A") !Current in Ampers - fluxPerStep = flow/qe + 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 - fluxPerStep = flow*self%surface*L_ref**2/qe + 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 @@ -184,17 +199,21 @@ MODULE moduleInject 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) = FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface /self%species%weight) - + self%particlesPerEdge(et) = max(1,FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface / self%species%weight)) END DO - END IF + self%nParticles = SUM(self%particlesPerEdge) - 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') @@ -232,23 +251,23 @@ MODULE moduleInject END SUBROUTINE doInjects - SUBROUTINE initVelDistMaxwellian(velDist, T, m) + SUBROUTINE initVelDistMaxwellian(velDist, temperature, m) IMPLICIT NONE CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist - REAL(8), INTENT(in):: T, m + REAL(8), INTENT(in):: temperature, m - velDist = velDistMaxwellian(vTh = DSQRT(T/m)) + velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m)) END SUBROUTINE initVelDistMaxwellian - SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m) + SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m) IMPLICIT NONE CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist - REAL(8), INTENT(in):: T, m + REAL(8), INTENT(in):: temperature, m - velDist = velDistHalfMaxwellian(vTh = DSQRT(T/m)) + velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m)) END SUBROUTINE initVelDistHalfMaxwellian @@ -270,7 +289,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = self%vTh*randomMaxwellian() + v = self%vTh*randomMaxwellian()/sqrt(2.d0) END FUNCTION randomVelMaxwellian @@ -283,10 +302,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - DO WHILE (v <= 0.D0) - v = self%vTh*randomMaxwellian() - - END DO + v = self%vTh*randomHalfMaxwellian()/sqrt(2.d0) END FUNCTION randomVelHalfMaxwellian @@ -361,19 +377,22 @@ MODULE moduleInject !Assign particle type partInj(n)%species => self%species - direction = self%n + if (all(self%n == 0.D0)) then + direction = randomEdge%normal + + else + direction = self%n + + end if partInj(n)%v = 0.D0 - partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & - self%v(2)%obj%randomVel(), & - self%v(3)%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 - !If velocity is not in the right direction, invert it - IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN - partInj(n)%v = - partInj(n)%v - - END IF !Obtain natural coordinates of particle in cell partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) diff --git a/src/modules/moduleProbe.f90 b/src/modules/moduleProbe.f90 index c29eeda..754d56a 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, timeStep) + SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, everyTimeStep) 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):: timeStep + REAL(8), INTENT(in):: everyTimeStep INTEGER:: sp, i REAL(8):: dv(1:3) @@ -91,11 +91,11 @@ MODULE moduleProbe 1:self%nv(3))) !Number of iterations between output - IF (timeStep == 0.D0) THEN + IF (everyTimeStep == 0.D0) THEN self%every = 1 ELSE - self%every = NINT(timeStep/ tauMin / ti_ref) + self%every = NINT(everyTimeStep/ tauMin / ti_ref) END IF @@ -189,13 +189,13 @@ MODULE moduleProbe END SUBROUTINE calculate - SUBROUTINE output(self, t) + SUBROUTINE output(self) 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) t + WRITE(tstring, iterationFormat) timeStep 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(t)*tauMin*ti_ref, " s" + WRITE(10, "(A6, ES15.6E3, A2)") "# t = ", REAL(timeStep)*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(t) + SUBROUTINE outputProbes() + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - INTEGER, INTENT(in):: t INTEGER:: i DO i = 1, nProbes IF (probe(i)%update) THEN - CALL probe(i)%output(t) + CALL probe(i)%output() END IF @@ -268,15 +268,15 @@ MODULE moduleProbe END SUBROUTINE outputProbes - SUBROUTINE resetProbes(t) + SUBROUTINE resetProbes() + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - INTEGER, INTENT(in):: t INTEGER:: i DO i = 1, nProbes probe(i)%f = 0.D0 - probe(i)%update = t == tFinal .OR. t == tInitial .OR. MOD(t, probe(i)%every) == 0 + probe(i)%update = timeStep == tFinal .OR. timeStep == tInitial .OR. MOD(timeStep, probe(i)%every) == 0 END DO diff --git a/src/modules/output/moduleOutput.f90 b/src/modules/output/moduleOutput.f90 index 18fbb7f..2af3c70 100644 --- a/src/modules/output/moduleOutput.f90 +++ b/src/modules/output/moduleOutput.f90 @@ -22,7 +22,6 @@ MODULE moduleOutput !Type for EM data in node TYPE emNode - CHARACTER(:), ALLOCATABLE:: type REAL(8):: phi REAL(8):: B(1:3) @@ -160,12 +159,12 @@ MODULE moduleOutput END SUBROUTINE calculateOutput - SUBROUTINE printTime(t, first) + SUBROUTINE printTime(first) USE moduleSpecies USE moduleCompTime + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - INTEGER, INTENT(in):: t LOGICAL, INTENT(in), OPTIONAL:: first CHARACTER(:), ALLOCATABLE:: fileName @@ -187,7 +186,7 @@ MODULE moduleOutput OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField + WRITE (20, "(I10, I10, 7(ES20.6E3))") timeStep, nPartOld, 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 d5d0793..7f6891c 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -1,56 +1,202 @@ !Module to solve the electromagnetic field MODULE moduleEM + USE moduleMesh + USE moduleTable IMPLICIT NONE - TYPE:: boundaryEM - CHARACTER(:), ALLOCATABLE:: typeEM - INTEGER:: physicalSurface + ! 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 REAL(8):: potential CONTAINS - PROCEDURE, PASS:: apply + ! boundaryEMGeneric DEFERRED PROCEDURES + PROCEDURE, PASS:: apply => applyDirichlet - END TYPE boundaryEM + 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 INTEGER:: nBoundaryEM - TYPE(boundaryEM), ALLOCATABLE:: boundEM(:) + TYPE(boundaryEMCont), ALLOCATABLE:: boundaryEM(:) !Information of charge and reference parameters for rho vector REAL(8), ALLOCATABLE:: qSpecies(:) CONTAINS - !Apply boundary conditions to the K matrix for Poisson's equation - SUBROUTINE apply(self, edge) + SUBROUTINE findNodes(self, physicalSurface) USE moduleMesh IMPLICIT NONE - CLASS(boundaryEM), INTENT(in):: self - CLASS(meshEdge):: edge - INTEGER:: nNodes - INTEGER, ALLOCATABLE:: nodes(:) - INTEGER:: n + CLASS(boundaryEMGeneric), INTENT(inout):: self + INTEGER, INTENT(in):: physicalSurface + CLASS(meshEdge), POINTER:: edge + INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:) + INTEGER:: nNodes, nodesNew + INTEGER:: e, n - nNodes = 1 - nNodes = edge%nNodes - nodes = edge%getNodes(nNodes) + !Temporal array to hold nodes + ALLOCATE(nodes(0)) - DO n = 1, nNodes - 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 + ! 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 - END SELECT + ELSE + ! If not, add element to array of nodes + nodes = [nodes, nodesEdge(n)] + + END IF + + END DO + + END IF END DO - END SUBROUTINE + ! 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 + + 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 !Assemble the source vector based on the charge density to solve Poisson's equation - SUBROUTINE assembleSourceVector(vectorF) + SUBROUTINE assembleSourceVector(vectorF, n_e) USE moduleMesh USE moduleRefParam IMPLICIT NONE @@ -59,8 +205,9 @@ 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 + INTEGER:: e, i, ni, b CLASS(meshNode), POINTER:: node !$OMP SINGLE @@ -78,6 +225,10 @@ 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 @@ -98,18 +249,12 @@ MODULE moduleEM !$OMP END DO !Apply boundary conditions - !$OMP DO - DO i = 1, mesh%numNodes - node => mesh%nodes(i)%obj + !$OMP SINGLE + do b = 1, nBoundaryEM + call boundaryEM(b)%obj%apply(vectorF) - SELECT CASE(node%emData%type) - CASE ("dirichlet") - vectorF(i) = node%emData%phi - - END SELECT - - END DO - !$OMP END DO + end do + !$OMP END SINGLE END SUBROUTINE assembleSourceVector @@ -157,4 +302,86 @@ 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 0b2837f..df8d5e8 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -138,6 +138,9 @@ MODULE moduleSolver CASE('Electrostatic','ConstantB') self%solveEM => solveElecField + CASE('ElectrostaticBoltzmann') + self%solveEM => solveElecFieldBoltzmann + END SELECT END SUBROUTINE initEM @@ -491,47 +494,46 @@ MODULE moduleSolver END SUBROUTINE updateParticleCell !Update the information about if a species needs to be moved this iteration - SUBROUTINE updatePushSpecies(self, t) + SUBROUTINE updatePushSpecies(self) 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(t, self%pusher(s)%every) == 0 + self%pusher(s)%pushSpecies = MOD(timeStep, self%pusher(s)%every) == 0 END DO END SUBROUTINE updatePushSpecies !Output the different data and information - SUBROUTINE doOutput(t) + SUBROUTINE doOutput() USE moduleMesh USE moduleOutput USE moduleSpecies USE moduleCompTime USE moduleProbe + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - INTEGER, INTENT(in):: t - - CALL outputProbes(t) + CALL outputProbes() counterOutput = counterOutput + 1 IF (counterOutput >= triggerOutput .OR. & - t == tFinal .OR. t == tInitial) THEN + timeStep == tFinal .OR. timeStep == tInitial) THEN !Resets output counter counterOutput=0 - 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 + CALL mesh%printOutput() + IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl() + CALL mesh%printEM() + WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld - IF (t == 0) THEN + IF (timeStep == 0) THEN WRITE(*, "(5X,A21,F8.1,A2)") " init time: ", 1.D3*tStep, "ms" ELSE @@ -549,34 +551,32 @@ MODULE moduleSolver counterCPUTime = counterCPUTime + 1 IF (counterCPUTime >= triggerCPUTime .OR. & - t == tFinal .OR. t == tInitial) THEN + timeStep == tFinal .OR. timeStep == tInitial) THEN !Reset CPU Time counter counterCPUTime = 0 - CALL printTime(t, t == 0) + CALL printTime(timeStep == 0) END IF !Output average values - IF (useAverage .AND. t == tFinal) THEN + IF (useAverage .AND. timeStep == tFinal) THEN CALL mesh%printAverage() END IF END SUBROUTINE doOutput - SUBROUTINE doAverage(t) + SUBROUTINE doAverage() USE moduleAverage USE moduleMesh IMPLICIT NONE - INTEGER, INTENT(in):: t INTEGER:: tAverage, n - IF (useAverage) THEN - tAverage = t - tAverageStart + tAverage = timeStep - tAverageStart IF (tAverage == 1) THEN !First iteration in which average scheme is used