Merge branch 'feature/electromagnetic' into 'development'

Restructuring pushsers and addition to first electromagnetic pusher

See merge request JorgeGonz/fpakc!27
This commit is contained in:
Jorge Gonzalez 2022-12-12 18:20:34 +00:00
commit 239552f715
66 changed files with 42726 additions and 3010 deletions

1
.gitignore vendored
View file

@ -4,4 +4,5 @@ obj/
doc/user_manual/
doc/coding_style/
json-fortran-8.2.0/
json-fortran/
runs/

204
data/collisions/EL_e-Ar.dat Normal file
View file

@ -0,0 +1,204 @@
# 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)
0.000000e+0 7.100000e-20
1.000000e-3 6.298400e-20
2.000000e-3 5.835000e-20
5.000000e-3 4.977500e-20
1.000000e-2 4.113000e-20
2.000000e-2 3.071300e-20
5.000000e-2 1.570900e-20
1.000000e-1 6.062200e-21
1.092000e-1 5.117000e-21
1.482000e-1 2.519100e-21
1.885000e-1 1.330700e-21
2.303000e-1 9.762900e-22
2.735000e-1 1.139000e-21
3.183000e-1 1.628900e-21
3.646000e-1 2.326700e-21
4.125000e-1 3.155200e-21
4.622000e-1 4.064500e-21
5.136000e-1 5.023100e-21
5.668000e-1 6.012300e-21
6.218000e-1 7.023300e-21
6.788000e-1 8.054600e-21
7.378000e-1 9.110500e-21
7.989000e-1 1.020000e-20
8.621000e-1 1.133700e-20
9.275000e-1 1.253700e-20
9.953000e-1 1.382100e-20
1.065400e+0 1.479100e-20
1.138000e+0 1.576700e-20
1.213100e+0 1.677000e-20
1.290900e+0 1.778100e-20
1.371400e+0 1.882800e-20
1.454700e+0 1.991100e-20
1.541000e+0 2.107400e-20
1.630300e+0 2.232400e-20
1.722700e+0 2.358000e-20
1.818400e+0 2.476000e-20
1.917400e+0 2.598200e-20
2.020000e+0 2.729100e-20
2.126100e+0 2.884100e-20
2.235900e+0 3.044500e-20
2.349700e+0 3.210500e-20
2.467400e+0 3.382400e-20
2.589200e+0 3.558500e-20
2.715400e+0 3.740100e-20
2.845900e+0 3.928100e-20
2.981100e+0 4.122700e-20
3.121000e+0 4.331500e-20
3.265800e+0 4.548700e-20
3.415700e+0 4.773600e-20
3.570900e+0 5.006300e-20
3.731500e+0 5.247300e-20
3.897800e+0 5.496700e-20
4.069900e+0 5.775100e-20
4.248100e+0 6.093800e-20
4.432500e+0 6.423700e-20
4.623400e+0 6.765200e-20
4.821000e+0 7.118700e-20
5.025600e+0 7.507600e-20
5.237300e+0 7.901500e-20
5.456500e+0 8.309200e-20
5.683400e+0 8.731200e-20
5.918300e+0 9.168100e-20
6.161400e+0 9.628400e-20
6.413100e+0 1.010900e-19
6.673600e+0 1.060800e-19
6.943300e+0 1.118000e-19
7.222400e+0 1.170000e-19
7.511400e+0 1.222000e-19
7.810500e+0 1.275900e-19
8.120100e+0 1.326900e-19
8.440600e+0 1.372100e-19
8.772400e+0 1.416500e-19
9.115800e+0 1.451600e-19
9.471300e+0 1.487100e-19
9.839300e+0 1.523900e-19
1.022020e+1 1.548800e-19
1.061450e+1 1.564600e-19
1.100000e+1 1.580000e-19
1.160000e+1 1.580000e-19
1.220000e+1 1.571000e-19
1.280380e+1 1.547800e-19
1.328890e+1 1.525600e-19
1.379110e+1 1.495700e-19
1.431090e+1 1.458200e-19
1.484890e+1 1.420600e-19
1.540590e+1 1.373500e-19
1.598240e+1 1.321600e-19
1.657920e+1 1.291500e-19
1.719700e+1 1.225700e-19
1.783650e+1 1.157400e-19
1.849840e+1 1.110100e-19
1.918370e+1 1.069000e-19
1.989300e+1 1.026400e-19
2.062720e+1 9.899000e-20
2.138720e+1 9.534100e-20
2.217390e+1 9.156500e-20
2.298830e+1 8.765600e-20
2.383130e+1 8.361000e-20
2.470400e+1 7.942100e-20
2.560730e+1 7.611800e-20
2.654230e+1 7.321900e-20
2.751020e+1 7.021800e-20
2.851210e+1 6.711300e-20
2.954920e+1 6.389700e-20
3.062280e+1 6.137900e-20
3.173410e+1 5.937900e-20
3.288440e+1 5.730800e-20
3.407520e+1 5.516500e-20
3.530780e+1 5.294600e-20
3.658370e+1 5.064900e-20
3.790450e+1 4.827200e-20
3.927170e+1 4.581100e-20
4.068690e+1 4.384700e-20
4.215190e+1 4.245600e-20
4.366840e+1 4.101500e-20
4.523810e+1 3.952400e-20
4.686300e+1 3.798000e-20
4.854500e+1 3.638200e-20
5.028610e+1 3.480000e-20
5.208840e+1 3.353800e-20
5.395410e+1 3.223200e-20
5.588530e+1 3.088000e-20
5.788440e+1 2.948100e-20
5.995370e+1 2.803200e-20
6.209570e+1 2.674300e-20
6.431310e+1 2.541200e-20
6.660830e+1 2.403500e-20
6.898420e+1 2.260900e-20
7.144360e+1 2.171100e-20
7.398940e+1 2.120200e-20
7.662470e+1 2.067500e-20
7.935260e+1 2.012900e-20
8.217640e+1 1.940100e-20
8.509940e+1 1.859800e-20
8.812510e+1 1.776600e-20
9.125710e+1 1.690400e-20
9.449930e+1 1.601300e-20
9.785530e+1 1.509000e-20
1.013293e+2 1.435400e-20
1.049254e+2 1.395800e-20
1.086478e+2 1.354900e-20
1.125011e+2 1.312500e-20
1.164898e+2 1.268600e-20
1.206186e+2 1.223200e-20
1.248925e+2 1.176200e-20
1.293167e+2 1.127500e-20
1.338963e+2 1.077100e-20
1.386368e+2 1.025000e-20
1.435440e+2 9.710200e-21
1.486236e+2 9.151400e-21
1.538817e+2 8.790400e-21
1.593245e+2 8.496500e-21
1.649587e+2 8.192200e-21
1.707908e+2 7.877300e-21
1.768279e+2 7.551300e-21
1.830772e+2 7.213800e-21
1.895461e+2 6.864500e-21
1.962423e+2 6.502900e-21
2.031738e+2 6.244500e-21
2.103489e+2 6.118900e-21
2.177762e+2 5.988900e-21
2.254644e+2 5.854400e-21
2.334229e+2 5.715100e-21
2.416610e+2 5.570900e-21
2.501886e+2 5.421700e-21
2.590160e+2 5.267200e-21
2.681535e+2 5.107300e-21
2.776121e+2 4.941800e-21
2.874032e+2 4.770400e-21
2.975383e+2 4.593100e-21
3.080295e+2 4.409500e-21
3.188895e+2 4.219400e-21
3.301311e+2 4.022700e-21
3.417678e+2 3.819100e-21
3.538134e+2 3.608300e-21
3.662823e+2 3.390100e-21
3.791894e+2 3.164200e-21
3.925501e+2 2.930400e-21
4.063803e+2 2.789400e-21
4.355158e+2 2.740800e-21
4.667351e+2 2.688800e-21
4.831724e+2 2.661400e-21
5.001872e+2 2.633000e-21
5.178000e+2 2.603700e-21
5.360318e+2 2.573300e-21
5.549043e+2 2.541800e-21
5.744399e+2 2.509300e-21
5.946621e+2 2.475600e-21
6.155950e+2 2.440700e-21
6.372635e+2 2.404600e-21
6.596934e+2 2.367200e-21
6.829116e+2 2.328500e-21
7.069458e+2 2.288400e-21
7.318245e+2 2.247000e-21
7.575776e+2 2.204000e-21
7.842356e+2 2.159600e-21
8.118305e+2 2.113600e-21
8.403951e+2 2.066000e-21
8.699636e+2 2.016700e-21
9.005711e+2 1.965700e-21
9.322543e+2 1.912900e-21
9.650509e+2 1.858200e-21

View file

@ -1,34 +1,203 @@
# H. C. Straub et. al, Physical Review A, 55,2(1995)
# 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)
17 1.700E-22
20 4.600E-21
25 1.240E-20
30 1.840E-20
35 2.260E-20
40 2.550E-20
45 2.660E-20
50 2.700E-20
55 2.690E-20
60 2.670E-20
65 2.670E-20
70 2.670E-20
75 2.660E-20
80 2.690E-20
85 2.700E-20
90 2.690E-20
95 2.670E-20
100 2.640E-20
110 2.610E-20
120 2.550E-20
140 2.450E-20
160 2.350E-20
180 2.270E-20
200 2.180E-20
225 2.100E-20
250 1.990E-20
275 1.870E-20
300 1.790E-20
350 1.630E-20
400 1.510E-20
450 1.390E-20
500 1.310E-20
1.570000e+1 0.000000e+0
1.571000e+1 1.033000e-25
1.574000e+1 3.631000e-23
1.577000e+1 7.390000e-23
1.581000e+1 1.128000e-22
1.585000e+1 1.531000e-22
1.589000e+1 1.948000e-22
1.593000e+1 2.379000e-22
1.597000e+1 2.826000e-22
1.602000e+1 3.330000e-22
1.606000e+1 3.914000e-22
1.611000e+1 4.518000e-22
1.616000e+1 5.143000e-22
1.621000e+1 5.791000e-22
1.627000e+1 6.461000e-22
1.632000e+1 7.155000e-22
1.638000e+1 7.873000e-22
1.644000e+1 8.616000e-22
1.650000e+1 9.386000e-22
1.656000e+1 1.026000e-21
1.663000e+1 1.116000e-21
1.670000e+1 1.209000e-21
1.677000e+1 1.306000e-21
1.684000e+1 1.406000e-21
1.691000e+1 1.510000e-21
1.699000e+1 1.617000e-21
1.707000e+1 1.733000e-21
1.715000e+1 1.853000e-21
1.724000e+1 1.977000e-21
1.733000e+1 2.106000e-21
1.742000e+1 2.239000e-21
1.752000e+1 2.378000e-21
1.762000e+1 2.526000e-21
1.772000e+1 2.680000e-21
1.783000e+1 2.839000e-21
1.794000e+1 3.004000e-21
1.805000e+1 3.175000e-21
1.817000e+1 3.354000e-21
1.829000e+1 3.540000e-21
1.842000e+1 3.731000e-21
1.855000e+1 3.933000e-21
1.868000e+1 4.146000e-21
1.882000e+1 4.367000e-21
1.897000e+1 4.596000e-21
1.912000e+1 4.837000e-21
1.927000e+1 5.089000e-21
1.943000e+1 5.349000e-21
1.960000e+1 5.618000e-21
1.977000e+1 5.897000e-21
1.995000e+1 6.186000e-21
2.013000e+1 6.498000e-21
2.032000e+1 6.826000e-21
2.052000e+1 7.161000e-21
2.073000e+1 7.464000e-21
2.094000e+1 7.777000e-21
2.116000e+1 8.092000e-21
2.138000e+1 8.414000e-21
2.162000e+1 8.757000e-21
2.186000e+1 9.122000e-21
2.211000e+1 9.468000e-21
2.237000e+1 9.786000e-21
2.264000e+1 1.013000e-20
2.292000e+1 1.050000e-20
2.321000e+1 1.085000e-20
2.351000e+1 1.121000e-20
2.382000e+1 1.158000e-20
2.414000e+1 1.197000e-20
2.447000e+1 1.237000e-20
2.482000e+1 1.278000e-20
2.517000e+1 1.317000e-20
2.554000e+1 1.355000e-20
2.592000e+1 1.400000e-20
2.631000e+1 1.440000e-20
2.672000e+1 1.479000e-20
2.715000e+1 1.519000e-20
2.758000e+1 1.560000e-20
2.804000e+1 1.604000e-20
2.850000e+1 1.650000e-20
2.899000e+1 1.699000e-20
2.949000e+1 1.749000e-20
3.001000e+1 1.801000e-20
3.055000e+1 1.844000e-20
3.111000e+1 1.888000e-20
3.168000e+1 1.935000e-20
3.228000e+1 1.981000e-20
3.290000e+1 2.027000e-20
3.354000e+1 2.075000e-20
3.420000e+1 2.123000e-20
3.488000e+1 2.167000e-20
3.559000e+1 2.214000e-20
3.633000e+1 2.255000e-20
3.709000e+1 2.289000e-20
3.787000e+1 2.324000e-20
3.869000e+1 2.351000e-20
3.953000e+1 2.376000e-20
4.040000e+1 2.398000e-20
4.131000e+1 2.416000e-20
4.224000e+1 2.435000e-20
4.321000e+1 2.454000e-20
4.421000e+1 2.474000e-20
4.525000e+1 2.492000e-20
4.632000e+1 2.501000e-20
4.743000e+1 2.509000e-20
4.858000e+1 2.519000e-20
4.978000e+1 2.528000e-20
5.101000e+1 2.544000e-20
5.228000e+1 2.562000e-20
5.360000e+1 2.580000e-20
5.497000e+1 2.600000e-20
5.639000e+1 2.617000e-20
5.785000e+1 2.634000e-20
5.937000e+1 2.652000e-20
6.094000e+1 2.673000e-20
6.256000e+1 2.696000e-20
6.425000e+1 2.719000e-20
6.599000e+1 2.738000e-20
6.779000e+1 2.752000e-20
6.965000e+1 2.767000e-20
7.159000e+1 2.786000e-20
7.358000e+1 2.806000e-20
7.565000e+1 2.823000e-20
7.780000e+1 2.831000e-20
8.001000e+1 2.840000e-20
8.231000e+1 2.845000e-20
8.468000e+1 2.849000e-20
8.714000e+1 2.854000e-20
8.969000e+1 2.859000e-20
9.232000e+1 2.860000e-20
9.505000e+1 2.860000e-20
9.788000e+1 2.854000e-20
1.008000e+2 2.848000e-20
1.038000e+2 2.842000e-20
1.070000e+2 2.836000e-20
1.102000e+2 2.830000e-20
1.136000e+2 2.823000e-20
1.170000e+2 2.816000e-20
1.206000e+2 2.807000e-20
1.243000e+2 2.788000e-20
1.282000e+2 2.769000e-20
1.322000e+2 2.753000e-20
1.363000e+2 2.741000e-20
1.406000e+2 2.727000e-20
1.450000e+2 2.705000e-20
1.496000e+2 2.682000e-20
1.543000e+2 2.654000e-20
1.592000e+2 2.625000e-20
1.643000e+2 2.598000e-20
1.696000e+2 2.572000e-20
1.750000e+2 2.545000e-20
1.807000e+2 2.516000e-20
1.865000e+2 2.478000e-20
1.925000e+2 2.439000e-20
1.988000e+2 2.398000e-20
2.052000e+2 2.367000e-20
2.119000e+2 2.337000e-20
2.189000e+2 2.307000e-20
2.260000e+2 2.275000e-20
2.335000e+2 2.243000e-20
2.412000e+2 2.209000e-20
2.491000e+2 2.174000e-20
2.574000e+2 2.142000e-20
2.659000e+2 2.110000e-20
2.747000e+2 2.076000e-20
2.839000e+2 2.041000e-20
2.933000e+2 2.005000e-20
3.031000e+2 1.969000e-20
3.132000e+2 1.935000e-20
3.237000e+2 1.899000e-20
3.346000e+2 1.862000e-20
3.458000e+2 1.824000e-20
3.575000e+2 1.791000e-20
3.695000e+2 1.759000e-20
3.820000e+2 1.727000e-20
3.949000e+2 1.693000e-20
4.083000e+2 1.659000e-20
4.221000e+2 1.623000e-20
4.364000e+2 1.585000e-20
4.512000e+2 1.548000e-20
4.666000e+2 1.520000e-20
4.824000e+2 1.492000e-20
4.989000e+2 1.462000e-20
5.159000e+2 1.435000e-20
5.335000e+2 1.406000e-20
5.517000e+2 1.377000e-20
5.706000e+2 1.347000e-20
5.901000e+2 1.316000e-20
6.104000e+2 1.285000e-20
6.313000e+2 1.256000e-20
6.530000e+2 1.226000e-20
6.754000e+2 1.194000e-20
6.986000e+2 1.162000e-20
7.226000e+2 1.137000e-20
7.475000e+2 1.112000e-20
7.733000e+2 1.087000e-20
7.999000e+2 1.060000e-20
8.275000e+2 1.039000e-20
8.561000e+2 1.018000e-20
8.857000e+2 9.958000e-21
9.163000e+2 9.736000e-21
9.480000e+2 9.514000e-21
9.808000e+2 9.285000e-21

148
data/collisions/IO_e-Li.dat Normal file
View file

@ -0,0 +1,148 @@
# Rusudan I., et al., Atoms 9(4):90 2021
# Relative energy (eV) cross section (m^2)
5.393 1.061E-23
6.000 5.885E-21
7.000 1.352E-20
8.000 1.927E-20
9.000 2.364E-20
10.000 2.698E-20
11.000 2.956E-20
12.000 3.156E-20
13.000 3.309E-20
14.000 3.427E-20
15.000 3.517E-20
16.000 3.585E-20
17.000 3.635E-20
18.000 3.670E-20
19.000 3.694E-20
20.000 3.708E-20
21.000 3.714E-20
22.000 3.714E-20
23.000 3.709E-20
24.000 3.700E-20
25.000 3.686E-20
26.000 3.671E-20
27.000 3.652E-20
28.000 3.632E-20
29.000 3.610E-20
30.000 3.588E-20
31.000 3.564E-20
32.000 3.539E-20
33.000 3.514E-20
34.000 3.488E-20
35.000 3.462E-20
36.000 3.435E-20
37.000 3.409E-20
38.000 3.383E-20
39.000 3.356E-20
40.000 3.330E-20
41.000 3.304E-20
42.000 3.277E-20
43.000 3.252E-20
44.000 3.226E-20
45.000 3.201E-20
46.000 3.175E-20
47.000 3.151E-20
48.000 3.126E-20
49.000 3.102E-20
50.000 3.078E-20
51.000 3.054E-20
52.000 3.031E-20
53.000 3.008E-20
54.000 2.985E-20
55.000 2.963E-20
56.000 2.941E-20
57.000 2.919E-20
58.000 2.898E-20
59.000 2.877E-20
60.000 2.856E-20
61.000 2.835E-20
62.000 2.815E-20
63.000 2.795E-20
64.000 2.776E-20
65.000 2.756E-20
66.000 2.737E-20
67.000 2.719E-20
68.000 2.700E-20
69.000 2.682E-20
70.000 2.664E-20
71.000 2.646E-20
72.000 2.629E-20
73.000 2.612E-20
74.000 2.595E-20
75.000 2.578E-20
76.000 2.562E-20
77.000 2.546E-20
78.000 2.530E-20
79.000 2.514E-20
80.000 2.499E-20
81.000 2.483E-20
82.000 2.468E-20
83.000 2.453E-20
84.000 2.439E-20
85.000 2.424E-20
86.000 2.410E-20
87.000 2.396E-20
88.000 2.382E-20
89.000 2.369E-20
90.000 2.355E-20
91.000 2.342E-20
92.000 2.329E-20
93.000 2.316E-20
94.000 2.303E-20
95.000 2.290E-20
96.000 2.278E-20
97.000 2.266E-20
98.000 2.253E-20
99.000 2.241E-20
100.000 2.230E-20
101.000 2.218E-20
102.000 2.206E-20
103.000 2.195E-20
104.000 2.184E-20
105.000 2.173E-20
106.000 2.162E-20
107.000 2.151E-20
108.000 2.140E-20
109.000 2.129E-20
110.000 2.119E-20
111.000 2.109E-20
112.000 2.098E-20
113.000 2.088E-20
114.000 2.078E-20
115.000 2.068E-20
116.000 2.059E-20
117.000 2.049E-20
118.000 2.039E-20
119.000 2.030E-20
120.000 2.021E-20
121.000 2.011E-20
122.000 2.002E-20
123.000 1.993E-20
124.000 1.984E-20
125.000 1.976E-20
126.000 1.967E-20
127.000 1.958E-20
128.000 1.950E-20
129.000 1.941E-20
130.000 1.933E-20
131.000 1.924E-20
132.000 1.916E-20
133.000 1.908E-20
134.000 1.900E-20
135.000 1.892E-20
136.000 1.884E-20
137.000 1.877E-20
138.000 1.869E-20
139.000 1.861E-20
140.000 1.854E-20
141.000 1.846E-20
142.000 1.839E-20
143.000 1.831E-20
144.000 1.824E-20
145.000 1.817E-20
146.000 1.810E-20
147.000 1.803E-20
148.000 1.796E-20
149.000 1.789E-20
150.000 1.782E-20

Binary file not shown.

View file

@ -47,6 +47,7 @@
\newglossaryentry{git}{name={Git},description={Git is a distributed version-control system for tracking changes in a set of files}}
\newglossaryentry{gitlab}{name={GitLab},description={GitLab is a web-based lifecycle tool that provides a \Gls{git}-repository manager}}
\newglossaryentry{gmsh}{name={Gmsh},description={A three-dimensional finite element mesh generator with built-in pre- and post-processing facilities.}}
\newglossaryentry{gnuplot}{name={Gnuplot},description={A portable command-line driven graphing utility for Linux, OS/2, MS Windows, OSX, VMS, and many other platforms.}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{bibliography}
@ -287,7 +288,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Gmsh}
Although \Gls{gmsh}\cite{gmshURL} is not required to compile and run \Gls{fpakc}, it is the default tool to generate finite element meshes and post-processing.
Right now, the only \acrshort{io} format available in \Gls{fpakc} is the v2.0 .msh format.
Right now, the only \acrshort{io} format available in \Gls{fpakc} is the v2.0 \lstinline|.msh| format.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Installation steps}
@ -340,7 +341,7 @@ make
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Mesh file}
\Gls{fpakc} accepts right now the version 2.0 of \Gls{gmsh} mesh format .msh in ASCII format.
\Gls{fpakc} accepts right now the version 2.0 of \Gls{gmsh} mesh format \lstinline|.msh| in ASCII format.
This file contains information about the nodes, edges and volumes that define the finite element mesh used by \Gls{fpakc} to scatter particle properties and compute the self-consistent electromagnetic field.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -405,28 +406,30 @@ make
The object \textbf{geometry} contains information about the type of geometry, the mesh file format and the mesh filename.
The accepted parameters are:
\begin{itemize}
\item \textbf{type}: Character.
\item \textbf{dimension}: Integer.
Number of spatial dimensions of the geometry.
Current values are: $0$, $1$, $2$ or $3$.
Zero dimension is a fictitious volume.
Geometry used mostly to test collisional effects.
No boundary or EM field is solved.
No injection can be implemented.
Initial state must be read from file.
No mesh file is required.
The optional argument \textbf{geometry.volume} can be used to set a value for the fictitious volume.
Otherwise, the volume is set to 1 in non-dimensional units.
\item \textbf{type}: Character.
Type of geometry.
Current accepted vaules are
\begin{itemize}
\item \textbf{3DCart}: Three-dimensional grid ($x \hyphen y \hyphen z$) in Cartesian coordinates..
For \Gls{gmsh} mesh format, the coordinates $x$, $y$ and $z$ correspond to $x$, $y$ and $z$ respectively.
\item \textbf{2DCyl}: Two-dimensional grid ($z \hyphen r$) with symmetry axis at $r = 0$.
For \Gls{gmsh} mesh format, the coordinates $x$ and $y$ correspond to $z$ and $r$ respectively.
\item \textbf{2DCart}: Two-dimensional grid ($x \hyphen y$) in Cartesian coordinates..
For \Gls{gmsh} mesh format, the coordinates $x$ and $y$ correspond to $x$ and $y$ respectively.
\item \textbf{1DRad}: One-dimensional grid ($r$) in radial coordinates
For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $r$.
\item \textbf{1DCart}: One-dimensional grid ($x$) in Cartesian coordinates
For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $x$.
\item \textbf{0D}: Zero dimension ficticius volume.
Geometry used mostly to test collisional effects.
No boundary or EM field is solved.
No injection can be implemented.
Initial state must be read from file.
No mesh file is required.
The optional argument \textbf{geometry.volume} can be used to set a ficticius volume.
Otherwise, the volume is set to 1 in non-dimensional units.
\item \textbf{Cart}: Cartesian coordinates.
Available for \textbf{geometry.dimension} $1$, $2$ and $3$.
For \Gls{gmsh} mesh format, the coordinates $x$, $y$ and $z$ correspond to $x$, $y$ and $z$ respectively.
\item \textbf{Cyl}: Cylindrical coordinates ($z \hyphen r$) with symmetry axis at $r = 0$.
Only available for \textbf{geometry.dimension} $2$.
For \Gls{gmsh} mesh format, the coordinates $x$ and $y$ correspond to $z$ and $r$ respectively.
\item \textbf{Rad}: One-dimensional radial space ($r$).
Only available for \textbf{geometry.dimension} $1$.
For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $r$.
\end{itemize}
\item \textbf{meshType}: Character.
Format of mesh file.
@ -437,9 +440,9 @@ make
\item \textbf{meshFile}: Character.
Mesh filename.
This file is searched in the path \textbf{output.path} and must contain the file extension.
\item \textbf{volume}: Real
\item \textbf{volume}: Real.
Units of $\unit{m^-3}$.
Used to set a ficticius volume for the \textbf{0D} geometry.
Used to set a fictitious volume for the $0$ dimension.
Ignored in the other cases.
\end{itemize}
@ -624,11 +627,12 @@ make
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{case}
This object determines the simulation time, time step, pushers, weighting scheme and solver for the electromagnetic field.
\subsection{solver}
This object determines the input parameters for the solvers used in the case, both for particle pushers and electromagnetic field.
Accepted variables are:
\begin{itemize}
\item \textbf{tau}: Real.
Units of $\unit{s}$.
Array dimension 'number of species'.
Defines the different time steps for each species.
Even if all time steps are equal, they need to be defined as an array.
@ -643,18 +647,9 @@ make
Array dimension 'number of species'.
Indicates the type of pusher used for each species:
\begin{itemize}
\item \textbf{3DCartNeutral}: Pushes particles in a 3D Cartesian space ($x \hyphen y \hyphen z$) without any external force.
\item \textbf{3DCartCharged}: Pushes particles in a 3D Cartesian space ($x \hyphen y \hyphen z$) including the effect of the electrostatic field.
\item \textbf{2DCylNeutral}: Pushes particles in a 2D cylindrical space ($z \hyphen r$) without any external force.
\item \textbf{2DCylCharged}: Pushes particles in a 2D cylindrical space ($z \hyphen r$) including the effect of the electrostatic field.
\item \textbf{2DCartNeutral}: Pushes particles in a 2D Cartesian space ($x \hyphen y$) without any external force.
\item \textbf{2DCartCharged}: Pushes particles in a 2D Cartesian space ($x \hyphen y$) including the effect of the electrostatic field.
\item \textbf{1DRadNeutral}: Pushes particles in a 1D cylindrical space ($r$) without any external force.
\item \textbf{1DRadCharged}: Pushes particles in a 1D cylindrical space ($r$) accounting the the electrostatic field.
\item \textbf{1DCartNeutral}: Pushes particles in a 1D Cartesian space ($x$) without any external force.
\item \textbf{1DCartCharged}: Pushes particles in a 1D Cartesian space ($x$) accounting the the electrostatic field.
\item \textbf{0D}: Dummy pusher for 0D geometry.
No pushing is actually done.
\item \textbf{Neutral}: Pushes a particle without any external force.
\item \textbf{Electrostatic}: Pushes a particle including the effect of the electrostatic field.
\item \textbf{Electromagnetic}: Pushes particles accounting for the electromagnetic field.
\end{itemize}
\item \textbf{WeightingScheme}: Character.
Indicates the variable weighting scheme to be used in the simulation.
@ -669,7 +664,13 @@ make
If no value is supplied, no field is solved.
\begin{itemize}
\item \textbf{Electrostatic}: Solves the Poison equation to obtain the self-consistent electrostatic potential.
\item \textbf{ConstantB}: Assumes a constant magnetic field in all the domain.
It solves the Poisson equation as in the \textbf{solver.EMSolver} option.
\end{itemize}
\item \textbf{B}: Real.
Units of $\unit{T}$.
Array of dimension $3$.
Provides the value of constant magnetic field for the option \textbf{solver.EMSolver} \textbf{ConstantB}.
\item \textbf{initial}: Array of objects.
Determines initial values for the species.
Required values are:
@ -743,6 +744,8 @@ make
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Example runs}\label{ch:exampleRuns}
This chapter presents a description of the different example files distributed with \acrshort{fpakc}.
All examples in the repository have a \lstinline|README.txt| file and a reference output.
Plotting of the output is done with \Gls{gnuplot} or \Gls{gmsh}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{1D Emissive Cathode (1D\_Cathode)}
@ -754,7 +757,7 @@ make
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{0D \ce{Ar}-\ce{Ar+} Elastic Collision (0D\_Argon)}
Test case to check the 0D geometry that includes the elastic collision between \ce{Ar} and \ce{Ar+}.
Initial states are read from the Argon\_Initial.dat and Argon+\_Initial.dat
Initial states are read from the \lstinline|Argon_Initial.dat| and \lstinline|Argon+_Initial.dat|.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{ALPHIE Grid system (ALPHIE\_Grid)}
@ -765,7 +768,7 @@ make
\section{Flow around cylinder (cylFlow)}
Simple case of neutral Argon flow around a cylinder in a 2D axial-symmetry geometry.
Elastic collisions between argon particles are included.
Two cases are presented here: one in which the same mesh (meshSingle.msh) for scattering and collisions is used (input.json) and a second one (inputDualMesh.json) in which a mesh is used for scattering (mesh.msh) and a second one is used only for collisions (meshColl.msh).
Two cases are presented here: one in which the same mesh (\lstinline|meshSingle.msh|) for scattering and collisions is used (\lstinline|input.json|) and a second one (\lstinline|inputDualMesh.json|) in which a mesh is used for scattering (\lstinline|mesh.msh|) and a second one is used only for collisions (\lstinline|meshColl.msh|).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\printglossaries

View file

@ -6,7 +6,7 @@ SRCDIR := $(TOPDIR)/src# source folder
# compiler
# gfortran:
FC := gfortran
JSONDIR := $(TOPDIR)/json-fortran-8.2.0/build-gfortran
JSONDIR := $(TOPDIR)/json-fortran/build-gfortran
# ifort:
# FC := ifort
# JSONDIR := $(TOPDIR)/json-fortran-8.2.0/build-ifort

11
runs/0D_Argon/README.txt Normal file
View file

@ -0,0 +1,11 @@
Example of 0D geometry.
Mostly used to test collision processes.
This example includes Ar and Ar+ with self-collisions for Ar and elastic collisions for AR-Ar+.
Different starting temperatures for each species that end up in equilibrium.
The gnuplpot script 'plot_Temperature.gp' generates a plot of the species temperatures.
The result is provided in the output folder to compare if modifications to the code are made.

View file

@ -1,7 +1,7 @@
{
"output": {
"path": "./runs/0D_test/",
"path": "./runs/0D_Argon/",
"triggerOutput": 1,
"numColl": true,
"folder": "test"
@ -13,21 +13,19 @@
"radius": 1.88e-10
},
"geometry": {
"type": "1DCart",
"meshType": "0D",
"dimension": 0,
"volume": 1e-11
},
"species": [
{"name": "Argon+", "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e0},
{"name": "Argon", "type": "neutral", "mass": 6.633e-26, "weight": 1.0e0}
],
"case": {
"tau": [1.0e-6, 1.0e-6],
"time": 1.0e-3,
"pusher": ["0D", "0D"],
"solver": {
"tau": [1.0e-3, 1.0e-3],
"finalTime": 1.0e0,
"initial": [
{"speciesName": "Argon+", "initialState": "Argon+_Initial.dat"},
{"speciesName": "Argon", "initialState": "Argon_Initial.dat"}
{"species": "Argon+", "file": "Argon+_Initial.dat"},
{"species": "Argon", "file": "Argon_Initial.dat"}
]
},
"interactions": {

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

View file

@ -43,6 +43,8 @@ set output "comp_temp.eps"
filename1 = "OUTPUT_Argon.dat"
filename2 = "OUTPUT_Argon+.dat"
folder = 'output/'
set lmargin at screen lmar
set rmargin at screen rmar
@ -58,7 +60,6 @@ set ylabel "Temperature (K)" offset 1.4,0.0
set key width 0.5 height 0.1 spacing 1.3 samplen 0.2 box opaque font ",16"
set key at graph 0.95, graph 0.9 right top
plot filename1 u ($1):($7) t "Ar" ls 1 with lines, \
filename2 u ($1):($7) t "Ar^{+}" ls 2 with lines, \
'< paste OUTPUT_Argon.dat OUTPUT_Argon+.dat' u ($1):($7 + $14) t "Sum" ls 3 with lines
plot folder.filename1 u ($1):($7) t "Ar" ls 1 with lines, \
folder.filename2 u ($1):($7) t "Ar^{+}" ls 2 with lines

View file

@ -1,51 +0,0 @@
#Element Density(m^-3) Velocity (m/2) Temperature (K)
1 1.8672E+15 -4.2433E+04 0.0000E+00 0.0000E+00 1.7536E+03
2 1.7491E+15 -4.1621E+04 0.0000E+00 0.0000E+00 1.8298E+03
3 1.7130E+15 -4.0665E+04 0.0000E+00 0.0000E+00 1.6838E+03
4 1.7346E+15 -3.9729E+04 0.0000E+00 0.0000E+00 1.5787E+03
5 1.7565E+15 -3.8817E+04 0.0000E+00 0.0000E+00 1.5641E+03
6 1.7107E+15 -3.7924E+04 0.0000E+00 0.0000E+00 1.4574E+03
7 1.6611E+15 -3.7036E+04 0.0000E+00 0.0000E+00 1.3561E+03
8 1.6543E+15 -3.6166E+04 0.0000E+00 0.0000E+00 1.2648E+03
9 1.6850E+15 -3.5302E+04 0.0000E+00 0.0000E+00 1.1812E+03
10 1.7615E+15 -3.4449E+04 0.0000E+00 0.0000E+00 1.1669E+03
11 1.7123E+15 -3.3627E+04 0.0000E+00 0.0000E+00 1.1458E+03
12 1.6138E+15 -3.2791E+04 0.0000E+00 0.0000E+00 1.0411E+03
13 1.6308E+15 -3.1965E+04 0.0000E+00 0.0000E+00 9.6919E+02
14 1.6581E+15 -3.1161E+04 0.0000E+00 0.0000E+00 9.8432E+02
15 1.6576E+15 -3.0359E+04 0.0000E+00 0.0000E+00 9.8741E+02
16 1.6780E+15 -2.9560E+04 0.0000E+00 0.0000E+00 9.1285E+02
17 1.6906E+15 -2.8768E+04 0.0000E+00 0.0000E+00 8.5525E+02
18 1.6773E+15 -2.7987E+04 0.0000E+00 0.0000E+00 8.5787E+02
19 1.6571E+15 -2.7208E+04 0.0000E+00 0.0000E+00 8.3556E+02
20 1.6547E+15 -2.6431E+04 0.0000E+00 0.0000E+00 7.9246E+02
21 1.6878E+15 -2.5659E+04 0.0000E+00 0.0000E+00 7.6169E+02
22 1.7458E+15 -2.4887E+04 0.0000E+00 0.0000E+00 7.4698E+02
23 1.7634E+15 -2.4124E+04 0.0000E+00 0.0000E+00 7.5560E+02
24 1.8070E+15 -2.3350E+04 0.0000E+00 0.0000E+00 7.3385E+02
25 1.8616E+15 -2.2579E+04 0.0000E+00 0.0000E+00 7.1475E+02
26 1.8601E+15 -2.1822E+04 0.0000E+00 0.0000E+00 7.1208E+02
27 1.8291E+15 -2.1053E+04 0.0000E+00 0.0000E+00 6.8693E+02
28 1.8542E+15 -2.0280E+04 0.0000E+00 0.0000E+00 6.5924E+02
29 1.9410E+15 -1.9518E+04 0.0000E+00 0.0000E+00 6.4794E+02
30 1.9685E+15 -1.8749E+04 0.0000E+00 0.0000E+00 6.4988E+02
31 2.0048E+15 -1.7973E+04 0.0000E+00 0.0000E+00 6.2752E+02
32 2.1060E+15 -1.7194E+04 0.0000E+00 0.0000E+00 6.1567E+02
33 2.1705E+15 -1.6418E+04 0.0000E+00 0.0000E+00 6.1351E+02
34 2.1947E+15 -1.5629E+04 0.0000E+00 0.0000E+00 6.0244E+02
35 2.2625E+15 -1.4824E+04 0.0000E+00 0.0000E+00 6.0286E+02
36 2.4096E+15 -1.4023E+04 0.0000E+00 0.0000E+00 5.9199E+02
37 2.5632E+15 -1.3209E+04 0.0000E+00 0.0000E+00 5.8388E+02
38 2.6774E+15 -1.2384E+04 0.0000E+00 0.0000E+00 5.8875E+02
39 2.8107E+15 -1.1544E+04 0.0000E+00 0.0000E+00 5.8973E+02
40 2.9409E+15 -1.0692E+04 0.0000E+00 0.0000E+00 5.9631E+02
41 3.1254E+15 -9.8219E+03 0.0000E+00 0.0000E+00 5.9943E+02
42 3.4330E+15 -8.9344E+03 0.0000E+00 0.0000E+00 5.9491E+02
43 3.7218E+15 -8.0260E+03 0.0000E+00 0.0000E+00 5.9949E+02
44 4.0857E+15 -7.0817E+03 0.0000E+00 0.0000E+00 6.4591E+02
45 4.5046E+15 -6.1015E+03 0.0000E+00 0.0000E+00 5.9669E+02
46 5.4123E+15 -5.0997E+03 0.0000E+00 0.0000E+00 5.7835E+02
47 6.7144E+15 -4.2183E+03 0.0000E+00 0.0000E+00 4.7454E+02
48 9.1344E+15 -3.1141E+03 0.0000E+00 0.0000E+00 3.0000E+02
49 1.4157E+16 -1.8530E+03 0.0000E+00 0.0000E+00 3.0000E+02
50 2.3064E+16 -9.7547E+02 0.0000E+00 0.0000E+00 2.7186E+02

View file

@ -1,51 +0,0 @@
#Element Density(m^-3) Velocity (m/2) Temperature (K)
1 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
2 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
3 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
4 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
5 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
6 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
7 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
8 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
9 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
10 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
11 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
12 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
13 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
14 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
15 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
16 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
17 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
18 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
19 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
20 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
21 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
22 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
23 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
24 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
25 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
26 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
27 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
28 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
29 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
30 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
31 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
32 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
33 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
34 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
35 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
36 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
37 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
38 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
39 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
40 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
41 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
42 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
43 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
44 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
45 7.5916E+12 1.8413E+05 0.0000E+00 0.0000E+00 5.7527E+03
46 7.5916E+12 1.8413E+05 0.0000E+00 0.0000E+00 5.7527E+03
47 4.0201E+14 -2.3775E+04 0.0000E+00 0.0000E+00 1.1359E+04
48 4.0201E+14 -2.3775E+04 0.0000E+00 0.0000E+00 1.1359E+04
49 6.9489E+15 -8.8002E+03 0.0000E+00 0.0000E+00 9.5125E+03
50 6.9489E+15 -8.8002E+03 0.0000E+00 0.0000E+00 9.5125E+03

View file

@ -0,0 +1,8 @@
This case is useful to explain how fpakc can work with different geometries (Cartisian and Radial in this case) using very similar input files and the same mesh.
From the position to the left, electrons are emitted at a constant rate.
The electric potential at the left is fixed at 0.
Depending on the geometry (which affects the symmetry) different distributions of the electron density and electric potential are achieved.
The last iteration obtained for these cases is in the output folder.

View file

@ -1,22 +0,0 @@
#Export the potential curves from Gmsh into a txt file.
reset
set terminal qt persist enhanced 1 size 600, 400
set style line 1 pt 1 lc rgb "red"
set style line 2 pt 2 lc rgb "blue"
set style line 3 pt 4 lc rgb "#006400"
set style line 4 pt 6 lc rgb "orange"
set style line 5 pt 8 lc rgb "black"
set style line 6 pt 10 lc rgb "#ADD8E6"
set xlabel "x/r (mm)"
set autoscale y
set ylabel "Potential (V)"
set key box at 90,-10
set title "Effect of geometry"
plot "Cartesian.dat" u ($5*1e3):($8) every 100 ls 1 t "Cartesian", \
"Radial.dat" u ($5*1e3):($8) every 100 ls 2 t "Radial"

View file

@ -14,21 +14,19 @@
"temperature": 11604.0
},
"geometry": {
"type": "1DCart",
"dimension": 1,
"type": "Cart",
"meshType": "gmsh2",
"meshFile": "mesh.msh"
},
"species": [
{"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1},
{"name": "Argon+", "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e1}
{"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1}
],
"boundary": [
{"name": "Cathode", "physicalSurface": 1, "bTypes": [
{"type": "absorption"},
{"type": "absorption"}
]},
{"name": "Infinite", "physicalSurface": 2, "bTypes": [
{"type": "transparent"},
{"type": "transparent"}
]}
],
@ -36,22 +34,14 @@
{"name": "Cathode", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1}
],
"inject": [
{"name": "Plasma Inf Ar+", "species": "Argon+", "flow": 1.00e-6, "units": "A", "v": 300.0, "T": [ 500.0, 500.0, 500.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2},
{"name": "Plasma Inf e", "species": "Electron", "flow": 2.64e-4, "units": "A", "v": 0.0, "T": [11604.0, 11604.0, 11604.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2},
{"name": "Plasma Cat e", "species": "Electron", "flow": 2.64e-5, "units": "A", "v": 180000.0, "T": [ 2300.0, 2300.0, 2300.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1}
],
"case": {
"tau": [1.0e-11, 1.0e-11],
"time": 4.0e-6,
"pusher": ["1DCartCharged", "1DCartCharged"],
"EMSolver": "Electrostatic",
"initial": [
{"speciesName": "Argon+", "initialState": "Argon+_Background.dat"},
{"speciesName": "Electron", "initialState": "Electron_Background.dat"}
]
"solver": {
"tau": [1.0e-11],
"finalTime": 1.0e-8,
"pusher": ["Electrostatic"],
"EMSolver": "Electrostatic"
},
"parallel": {
"OpenMP":{

View file

@ -14,21 +14,19 @@
"temperature": 11604.0
},
"geometry": {
"type": "1DRad",
"dimension": 1,
"type": "Rad",
"meshType": "gmsh2",
"meshFile": "mesh.msh"
},
"species": [
{"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1},
{"name": "Argon+", "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e1}
{"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1}
],
"boundary": [
{"name": "Cathode", "physicalSurface": 1, "bTypes": [
{"type": "absorption"},
{"type": "absorption"}
]},
{"name": "Infinite", "physicalSurface": 2, "bTypes": [
{"type": "transparent"},
{"type": "transparent"}
]}
],
@ -36,22 +34,14 @@
{"name": "Cathode", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1}
],
"inject": [
{"name": "Plasma Inf Ar+", "species": "Argon+", "flow": 1.00e-6, "units": "A", "v": 323.0, "T": [ 500.0, 500.0, 500.0],
"velDist": ["Delta", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2},
{"name": "Plasma Inf e", "species": "Electron", "flow": 2.64e-4, "units": "A", "v": 0.0, "T": [11604.0, 11604.0, 11604.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2},
{"name": "Plasma Cat e", "species": "Electron", "flow": 2.64e-2, "units": "A", "v": 180000.0, "T": [ 2300.0, 2300.0, 2300.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1}
],
"case": {
"tau": [1.0e-11, 1.0e-11],
"time": 3.0e-7,
"pusher": ["1DRadCharged", "1DRadCharged"],
"EMSolver": "Electrostatic",
"initial": [
{"speciesName": "Argon+", "initialState": "Argon+_Background.dat"},
{"speciesName": "Electron", "initialState": "Electron_Background.dat"}
]
"solver": {
"tau": [1.0e-11],
"finalTime": 1.0e-8,
"pusher": ["Electrostatic"],
"EMSolver": "Electrostatic"
},
"parallel": {
"OpenMP":{

View file

@ -0,0 +1,185 @@
$MeshFormat
2.2 0 8
$EndMeshFormat
$NodeData
1
"Potential (V)"
1
1.0000000000000000E-008
3
1000
1
51
1 -7.2180834934653183E-012
2 -1.6164415181357834
3 -0.36980461527280883
4 -0.57492573345192088
5 -0.71403100967357258
6 -0.82100545766928923
7 -0.90752842743341133
8 -0.97972921034487059
9 -1.0411791818607907
10 -1.0939318279574932
11 -1.1399404313046071
12 -1.1807139724083482
13 -1.2174163304531369
14 -1.2504990809463861
15 -1.2803540949915220
16 -1.3073500616819811
17 -1.3316282169084139
18 -1.3537385275137412
19 -1.3739696197563589
20 -1.3927592194128697
21 -1.4103076231740377
22 -1.4266496533259061
23 -1.4418485519615323
24 -1.4562354448986332
25 -1.4697852468324024
26 -1.4825221356246889
27 -1.4946048104025964
28 -1.5059524644639768
29 -1.5166545901603037
30 -1.5269837565446036
31 -1.5366931422515575
32 -1.5457467079858906
33 -1.5541298965687140
34 -1.5618551829370724
35 -1.5689172611267244
36 -1.5752150587183313
37 -1.5808462624732060
38 -1.5860612571248669
39 -1.5910404498436017
40 -1.5957279636246224
41 -1.5998828959091791
42 -1.6032394878414955
43 -1.6059267946667977
44 -1.6082423003097928
45 -1.6102897646308323
46 -1.6120678175170864
47 -1.6134998424675298
48 -1.6146588662864141
49 -1.6155317189248377
50 -1.6160774066120194
51 -1.6163242439958503
$EndNodeData
$ElementData
1
"Electric Field (V m^-1)"
1
1.0000000000000000E-008
3
1000
3
50
3 9245.1153816612405 0.0000000000000000 0.0000000000000000
4 5128.0279544843424 0.0000000000000000 0.0000000000000000
5 3477.6319055457411 0.0000000000000000 0.0000000000000000
6 2674.3611998963270 0.0000000000000000 0.0000000000000000
7 2163.0742441058046 0.0000000000000000 0.0000000000000000
8 1805.0195727920382 0.0000000000000000 0.0000000000000000
9 1536.2492879082888 0.0000000000000000 0.0000000000000000
10 1318.8161524257080 0.0000000000000000 0.0000000000000000
11 1150.2150836799428 0.0000000000000000 0.0000000000000000
12 1019.3385275948289 0.0000000000000000 0.0000000000000000
13 917.55895112089058 0.0000000000000000 0.0000000000000000
14 827.06876233228445 0.0000000000000000 0.0000000000000000
15 746.37535112673595 0.0000000000000000 0.0000000000000000
16 674.89916726183060 0.0000000000000000 0.0000000000000000
17 606.95388066159535 0.0000000000000000 0.0000000000000000
18 552.75776513387154 0.0000000000000000 0.0000000000000000
19 505.77730606825133 0.0000000000000000 0.0000000000000000
20 469.73999141591514 0.0000000000000000 0.0000000000000000
21 438.71009403214083 0.0000000000000000 0.0000000000000000
22 408.55075379944572 0.0000000000000000 0.0000000000000000
23 379.97246589319053 0.0000000000000000 0.0000000000000000
24 359.67232342854714 0.0000000000000000 0.0000000000000000
25 338.74504834466700 0.0000000000000000 0.0000000000000000
26 318.42221980693205 0.0000000000000000 0.0000000000000000
27 302.06686944643235 0.0000000000000000 0.0000000000000000
28 283.69135153333235 0.0000000000000000 0.0000000000000000
29 267.55314240706213 0.0000000000000000 0.0000000000000000
30 258.22915960642655 0.0000000000000000 0.0000000000000000
31 242.73464267283171 0.0000000000000000 0.0000000000000000
32 226.33914335739505 0.0000000000000000 0.0000000000000000
33 209.57971456971421 0.0000000000000000 0.0000000000000000
34 193.13215920815529 0.0000000000000000 0.0000000000000000
35 176.55195474056683 0.0000000000000000 0.0000000000000000
36 157.44493979054459 0.0000000000000000 0.0000000000000000
37 140.78009387281193 0.0000000000000000 0.0000000000000000
38 130.37486629147924 0.0000000000000000 0.0000000000000000
39 124.47981796785177 0.0000000000000000 0.0000000000000000
40 117.18784452503036 0.0000000000000000 0.0000000000000000
41 103.87330711348883 0.0000000000000000 0.0000000000000000
42 83.914798307876424 0.0000000000000000 0.0000000000000000
43 67.182670633004349 0.0000000000000000 0.0000000000000000
44 57.887641075013939 0.0000000000000000 0.0000000000000000
45 51.186608025777616 0.0000000000000000 0.0000000000000000
46 44.451322156171258 0.0000000000000000 0.0000000000000000
47 35.800623760933746 0.0000000000000000 0.0000000000000000
48 28.975595471990374 0.0000000000000000 0.0000000000000000
49 21.821315960497699 0.0000000000000000 0.0000000000000000
50 13.642192179482622 0.0000000000000000 0.0000000000000000
51 6.1709345957464548 0.0000000000000000 0.0000000000000000
52 2.9318534983186093 0.0000000000000000 0.0000000000000000
$EndElementData
$NodeData
1
"Magnetic Field (T)"
1
1.0000000000000000E-008
3
1000
3
51
1 0.0000000000000000 0.0000000000000000 0.0000000000000000
2 0.0000000000000000 0.0000000000000000 0.0000000000000000
3 0.0000000000000000 0.0000000000000000 0.0000000000000000
4 0.0000000000000000 0.0000000000000000 0.0000000000000000
5 0.0000000000000000 0.0000000000000000 0.0000000000000000
6 0.0000000000000000 0.0000000000000000 0.0000000000000000
7 0.0000000000000000 0.0000000000000000 0.0000000000000000
8 0.0000000000000000 0.0000000000000000 0.0000000000000000
9 0.0000000000000000 0.0000000000000000 0.0000000000000000
10 0.0000000000000000 0.0000000000000000 0.0000000000000000
11 0.0000000000000000 0.0000000000000000 0.0000000000000000
12 0.0000000000000000 0.0000000000000000 0.0000000000000000
13 0.0000000000000000 0.0000000000000000 0.0000000000000000
14 0.0000000000000000 0.0000000000000000 0.0000000000000000
15 0.0000000000000000 0.0000000000000000 0.0000000000000000
16 0.0000000000000000 0.0000000000000000 0.0000000000000000
17 0.0000000000000000 0.0000000000000000 0.0000000000000000
18 0.0000000000000000 0.0000000000000000 0.0000000000000000
19 0.0000000000000000 0.0000000000000000 0.0000000000000000
20 0.0000000000000000 0.0000000000000000 0.0000000000000000
21 0.0000000000000000 0.0000000000000000 0.0000000000000000
22 0.0000000000000000 0.0000000000000000 0.0000000000000000
23 0.0000000000000000 0.0000000000000000 0.0000000000000000
24 0.0000000000000000 0.0000000000000000 0.0000000000000000
25 0.0000000000000000 0.0000000000000000 0.0000000000000000
26 0.0000000000000000 0.0000000000000000 0.0000000000000000
27 0.0000000000000000 0.0000000000000000 0.0000000000000000
28 0.0000000000000000 0.0000000000000000 0.0000000000000000
29 0.0000000000000000 0.0000000000000000 0.0000000000000000
30 0.0000000000000000 0.0000000000000000 0.0000000000000000
31 0.0000000000000000 0.0000000000000000 0.0000000000000000
32 0.0000000000000000 0.0000000000000000 0.0000000000000000
33 0.0000000000000000 0.0000000000000000 0.0000000000000000
34 0.0000000000000000 0.0000000000000000 0.0000000000000000
35 0.0000000000000000 0.0000000000000000 0.0000000000000000
36 0.0000000000000000 0.0000000000000000 0.0000000000000000
37 0.0000000000000000 0.0000000000000000 0.0000000000000000
38 0.0000000000000000 0.0000000000000000 0.0000000000000000
39 0.0000000000000000 0.0000000000000000 0.0000000000000000
40 0.0000000000000000 0.0000000000000000 0.0000000000000000
41 0.0000000000000000 0.0000000000000000 0.0000000000000000
42 0.0000000000000000 0.0000000000000000 0.0000000000000000
43 0.0000000000000000 0.0000000000000000 0.0000000000000000
44 0.0000000000000000 0.0000000000000000 0.0000000000000000
45 0.0000000000000000 0.0000000000000000 0.0000000000000000
46 0.0000000000000000 0.0000000000000000 0.0000000000000000
47 0.0000000000000000 0.0000000000000000 0.0000000000000000
48 0.0000000000000000 0.0000000000000000 0.0000000000000000
49 0.0000000000000000 0.0000000000000000 0.0000000000000000
50 0.0000000000000000 0.0000000000000000 0.0000000000000000
51 0.0000000000000000 0.0000000000000000 0.0000000000000000
$EndNodeData

View file

@ -0,0 +1,247 @@
$MeshFormat
2.2 0 8
$EndMeshFormat
$NodeData
1
"Electron density (m^-3)"
1
1.0000000000000000E-008
3
1000
1
51
1 1.663355E+017
2 1.635218E+014
3 6.971311E+016
4 2.630990E+016
5 1.405426E+016
6 9.465391E+015
7 6.637877E+015
8 4.988303E+015
9 4.189069E+015
10 3.156337E+015
11 2.494143E+015
12 1.855358E+015
13 1.740443E+015
14 1.546003E+015
15 1.316731E+015
16 1.372680E+015
17 9.737899E+014
18 9.388176E+014
19 6.512445E+014
20 5.832764E+014
21 5.692614E+014
22 5.935914E+014
23 3.292208E+014
24 4.143419E+014
25 4.100487E+014
26 2.728780E+014
27 3.714897E+014
28 3.455677E+014
29 9.442239E+013
30 3.445456E+014
31 3.018612E+014
32 3.256579E+014
33 3.148337E+014
34 2.986174E+014
35 3.894979E+014
36 3.315681E+014
37 1.927240E+014
38 8.916509E+013
39 1.257294E+014
40 2.430101E+014
41 4.270402E+014
42 3.345209E+014
43 1.510759E+014
44 1.256626E+014
45 1.136905E+014
46 1.909151E+014
47 1.133460E+014
48 1.373180E+014
49 1.567058E+014
50 1.725501E+014
51 8.718697E+012
$EndNodeData
$NodeData
1
"Electron velocity (m s^-1)"
1
1.0000000000000000E-008
3
1000
3
51
1 6.822573E+002 7.341184E+002 -1.294679E+003
2 2.275900E+005 -1.156323E+005 -3.994264E+004
3 6.011992E+003 -2.085702E+003 -6.955356E+003
4 2.674256E+003 -1.413909E+003 1.028663E+003
5 -6.262168E+003 -5.034055E+003 1.381887E+004
6 1.057957E+003 8.436700E+002 -1.017065E+004
7 1.182090E+004 1.797951E+004 6.390618E+003
8 1.868273E+004 1.131158E+004 1.107675E+004
9 1.821590E+004 5.287286E+003 9.958486E+003
10 5.545844E+003 -1.134700E+003 -1.747390E+004
11 4.595683E+004 1.909957E+004 4.055354E+004
12 5.242670E+004 -1.209486E+003 1.046959E+004
13 -6.508610E+003 4.362341E+004 -8.868030E+003
14 -5.779804E+003 7.888582E+003 -1.190795E+004
15 -3.256102E+003 -2.423327E+004 1.680404E+004
16 2.026515E+004 5.900866E+003 4.573129E+004
17 1.881938E+004 -5.907004E+004 4.787712E+003
18 2.471440E+004 3.126628E+004 3.276285E+004
19 1.927866E+004 -1.774386E+004 1.728677E+004
20 4.357578E+004 -1.754746E+004 2.441953E+004
21 -2.293113E+004 5.501146E+003 3.536741E+004
22 2.289073E+004 2.308340E+004 5.214795E+004
23 1.380393E+005 8.134296E+004 6.839466E+004
24 1.059159E+005 7.461802E+003 -9.633646E+003
25 1.093199E+004 -1.991657E+004 -1.139848E+005
26 2.642203E+004 -2.481740E+004 3.143777E+004
27 9.005040E+004 8.070628E+003 -3.359234E+004
28 1.241325E+005 5.096109E+004 7.419358E+002
29 1.346489E+005 -1.573873E+005 6.203455E+003
30 5.610495E+004 -3.368553E+004 1.507073E+005
31 3.535163E+004 1.988632E+004 6.165400E+004
32 9.376597E+004 1.073849E+004 1.451017E+004
33 1.084345E+005 3.421118E+004 2.083726E+004
34 7.746595E+004 3.203130E+004 6.852780E+003
35 1.133766E+005 2.337750E+003 -4.867282E+004
36 1.025276E+005 1.014552E+005 -8.405875E+004
37 7.241202E+004 2.972254E+002 1.876103E+004
38 4.204461E+004 9.161183E+004 -8.625228E+004
39 2.000094E+004 1.616608E+005 -4.124777E+004
40 1.065082E+005 1.060744E+005 -9.625557E+004
41 2.353957E+005 -3.660256E+003 -6.119245E+003
42 2.043649E+005 1.172239E+003 -2.283045E+004
43 1.632956E+005 -3.766183E+004 -5.714495E+004
44 2.139853E+005 -1.091968E+005 -7.473445E+004
45 6.523701E+004 -5.437454E+004 2.092413E+005
46 6.804248E+004 2.190917E+004 -3.051095E+004
47 1.247063E+005 9.312863E+004 -2.110944E+004
48 1.488446E+005 4.703577E+003 1.162251E+005
49 1.444633E+005 2.511461E+004 7.188944E+004
50 1.343558E+005 -1.246722E+004 2.197029E+005
51 1.789699E+005 -1.338478E+005 -5.783810E+004
$EndNodeData
$NodeData
1
"Electron pressure (Pa)"
1
1.0000000000000000E-008
3
1000
1
51
1 5.385580E-003
2 1.002571E-006
3 2.356496E-003
4 8.962727E-004
5 4.794124E-004
6 3.503273E-004
7 2.442948E-004
8 1.546660E-004
9 1.322598E-004
10 9.957574E-005
11 8.063749E-005
12 5.468381E-005
13 5.120468E-005
14 4.099773E-005
15 4.147335E-005
16 4.439409E-005
17 3.060211E-005
18 3.239301E-005
19 1.928871E-005
20 1.720661E-005
21 1.615443E-005
22 2.301910E-005
23 9.019928E-006
24 8.246059E-006
25 9.022171E-006
26 5.090623E-006
27 1.430550E-005
28 1.078930E-005
29 2.237441E-006
30 8.720790E-006
31 7.185960E-006
32 1.040168E-005
33 8.947370E-006
34 8.505472E-006
35 1.011847E-005
36 1.009116E-005
37 6.487404E-006
38 9.794017E-007
39 1.137051E-006
40 4.115733E-006
41 1.110521E-005
42 5.195716E-006
43 2.249886E-006
44 2.702675E-006
45 1.909427E-006
46 4.939937E-006
47 1.983386E-006
48 2.345377E-006
49 4.284590E-006
50 4.071858E-006
51 4.700991E-008
$EndNodeData
$NodeData
1
"Electron temperature (K)"
1
1.0000000000000000E-008
3
1000
1
51
1 2.345117E+003
2 4.440750E+002
3 2.448325E+003
4 2.467391E+003
5 2.470690E+003
6 2.680725E+003
7 2.665642E+003
8 2.245737E+003
9 2.286795E+003
10 2.285004E+003
11 2.341707E+003
12 2.134754E+003
13 2.130918E+003
14 1.920730E+003
15 2.281334E+003
16 2.342463E+003
17 2.276161E+003
18 2.499120E+003
19 2.145241E+003
20 2.136672E+003
21 2.055403E+003
22 2.808779E+003
23 1.984416E+003
24 1.441466E+003
25 1.593648E+003
26 1.351199E+003
27 2.789158E+003
28 2.261399E+003
29 1.716301E+003
30 1.833268E+003
31 1.724227E+003
32 2.313443E+003
33 2.058406E+003
34 2.063004E+003
35 1.881598E+003
36 2.204374E+003
37 2.438102E+003
38 7.955782E+002
39 6.550284E+002
40 1.226704E+003
41 1.883539E+003
42 1.124965E+003
43 1.078654E+003
44 1.557774E+003
45 1.216454E+003
46 1.874123E+003
47 1.267412E+003
48 1.237093E+003
49 1.980347E+003
50 1.709206E+003
51 3.905303E+002
$EndNodeData

View file

@ -0,0 +1,185 @@
$MeshFormat
2.2 0 8
$EndMeshFormat
$NodeData
1
"Potential (V)"
1
1.0000000000000000E-008
3
1000
1
51
1 -2.0125602360500233E-012
2 -2.1476540511245403
3 -0.92254600934286035
4 -1.1693722452494442
5 -1.3055892978798838
6 -1.4090596211031350
7 -1.4933143808451328
8 -1.5638362561671308
9 -1.6231896643635058
10 -1.6740175110100952
11 -1.7185366931231973
12 -1.7583784348355216
13 -1.7946512154913346
14 -1.8278278087327231
15 -1.8580164245432433
16 -1.8852831251773992
17 -1.9099218127680115
18 -1.9323133093677307
19 -1.9527025003199079
20 -1.9712632206620335
21 -1.9881621423695754
22 -2.0035781018067564
23 -2.0176479044718039
24 -2.0305087603061467
25 -2.0422748378993707
26 -2.0530348175074842
27 -2.0628746022811000
28 -2.0718874967016756
29 -2.0801381491382527
30 -2.0876833658364853
31 -2.0945891263065528
32 -2.1008985562257543
33 -2.1066381420356524
34 -2.1118556546612779
35 -2.1165897969031175
36 -2.1208699819402956
37 -2.1247303076336199
38 -2.1281989125689238
39 -2.1312983972494055
40 -2.1340576795392843
41 -2.1364955581097296
42 -2.1386339190704353
43 -2.1405011196666606
44 -2.1421213004925956
45 -2.1435109859394803
46 -2.1446785599633511
47 -2.1456400594108116
48 -2.1464001968748745
49 -2.1469649987104309
50 -2.1473539773488697
51 -2.1475805166422961
$EndNodeData
$ElementData
1
"Electric Field (V m^-1)"
1
1.0000000000000000E-008
3
1000
3
50
3 23063.650233574754 0.0000000000000000 0.0000000000000000
4 6170.6558976724700 0.0000000000000000 0.0000000000000000
5 3405.4263157653404 0.0000000000000000 0.0000000000000000
6 2586.7580805845814 0.0000000000000000 0.0000000000000000
7 2106.3689935526277 0.0000000000000000 0.0000000000000000
8 1763.0468830553750 0.0000000000000000 0.0000000000000000
9 1483.8352049193140 0.0000000000000000 0.0000000000000000
10 1270.6961661725816 0.0000000000000000 0.0000000000000000
11 1112.9795528295804 0.0000000000000000 0.0000000000000000
12 996.04354280937548 0.0000000000000000 0.0000000000000000
13 906.81951639648162 0.0000000000000000 0.0000000000000000
14 829.41483103577525 0.0000000000000000 0.0000000000000000
15 754.71539526132392 0.0000000000000000 0.0000000000000000
16 681.66751585425504 0.0000000000000000 0.0000000000000000
17 615.96718976609407 0.0000000000000000 0.0000000000000000
18 559.78741499368221 0.0000000000000000 0.0000000000000000
19 509.72977380725575 0.0000000000000000 0.0000000000000000
20 464.01800855624742 0.0000000000000000 0.0000000000000000
21 422.47304269137754 0.0000000000000000 0.0000000000000000
22 385.39898593210717 0.0000000000000000 0.0000000000000000
23 351.74506662853850 0.0000000000000000 0.0000000000000000
24 321.52139585948811 0.0000000000000000 0.0000000000000000
25 294.15193983097186 0.0000000000000000 0.0000000000000000
26 268.99949020263892 0.0000000000000000 0.0000000000000000
27 245.99461933938235 0.0000000000000000 0.0000000000000000
28 225.32236051344728 0.0000000000000000 0.0000000000000000
29 206.26631091357709 0.0000000000000000 0.0000000000000000
30 188.63041745502420 0.0000000000000000 0.0000000000000000
31 172.64401175097743 0.0000000000000000 0.0000000000000000
32 157.73574797938116 0.0000000000000000 0.0000000000000000
33 143.48964524685633 0.0000000000000000 0.0000000000000000
34 130.43781564009393 0.0000000000000000 0.0000000000000000
35 118.35355604550293 0.0000000000000000 0.0000000000000000
36 107.00462592970442 0.0000000000000000 0.0000000000000000
37 96.508142333755885 0.0000000000000000 0.0000000000000000
38 86.715123382567057 0.0000000000000000 0.0000000000000000
39 77.487117011718453 0.0000000000000000 0.0000000000000000
40 68.982057246684732 0.0000000000000000 0.0000000000000000
41 60.946964260881018 0.0000000000000000 0.0000000000000000
42 53.459024017616187 0.0000000000000000 0.0000000000000000
43 46.680014905952270 0.0000000000000000 0.0000000000000000
44 40.504520648464592 0.0000000000000000 0.0000000000000000
45 34.742136171983248 0.0000000000000000 0.0000000000000000
46 29.189350596637517 0.0000000000000000 0.0000000000000000
47 24.037486186421319 0.0000000000000000 0.0000000000000000
48 19.003436601487543 0.0000000000000000 0.0000000000000000
49 14.120045888854746 0.0000000000000000 0.0000000000000000
50 9.7244659609280415 0.0000000000000000 0.0000000000000000
51 5.6634823356386530 0.0000000000000000 0.0000000000000000
52 1.8383620561016323 0.0000000000000000 0.0000000000000000
$EndElementData
$NodeData
1
"Magnetic Field (T)"
1
1.0000000000000000E-008
3
1000
3
51
1 0.0000000000000000 0.0000000000000000 0.0000000000000000
2 0.0000000000000000 0.0000000000000000 0.0000000000000000
3 0.0000000000000000 0.0000000000000000 0.0000000000000000
4 0.0000000000000000 0.0000000000000000 0.0000000000000000
5 0.0000000000000000 0.0000000000000000 0.0000000000000000
6 0.0000000000000000 0.0000000000000000 0.0000000000000000
7 0.0000000000000000 0.0000000000000000 0.0000000000000000
8 0.0000000000000000 0.0000000000000000 0.0000000000000000
9 0.0000000000000000 0.0000000000000000 0.0000000000000000
10 0.0000000000000000 0.0000000000000000 0.0000000000000000
11 0.0000000000000000 0.0000000000000000 0.0000000000000000
12 0.0000000000000000 0.0000000000000000 0.0000000000000000
13 0.0000000000000000 0.0000000000000000 0.0000000000000000
14 0.0000000000000000 0.0000000000000000 0.0000000000000000
15 0.0000000000000000 0.0000000000000000 0.0000000000000000
16 0.0000000000000000 0.0000000000000000 0.0000000000000000
17 0.0000000000000000 0.0000000000000000 0.0000000000000000
18 0.0000000000000000 0.0000000000000000 0.0000000000000000
19 0.0000000000000000 0.0000000000000000 0.0000000000000000
20 0.0000000000000000 0.0000000000000000 0.0000000000000000
21 0.0000000000000000 0.0000000000000000 0.0000000000000000
22 0.0000000000000000 0.0000000000000000 0.0000000000000000
23 0.0000000000000000 0.0000000000000000 0.0000000000000000
24 0.0000000000000000 0.0000000000000000 0.0000000000000000
25 0.0000000000000000 0.0000000000000000 0.0000000000000000
26 0.0000000000000000 0.0000000000000000 0.0000000000000000
27 0.0000000000000000 0.0000000000000000 0.0000000000000000
28 0.0000000000000000 0.0000000000000000 0.0000000000000000
29 0.0000000000000000 0.0000000000000000 0.0000000000000000
30 0.0000000000000000 0.0000000000000000 0.0000000000000000
31 0.0000000000000000 0.0000000000000000 0.0000000000000000
32 0.0000000000000000 0.0000000000000000 0.0000000000000000
33 0.0000000000000000 0.0000000000000000 0.0000000000000000
34 0.0000000000000000 0.0000000000000000 0.0000000000000000
35 0.0000000000000000 0.0000000000000000 0.0000000000000000
36 0.0000000000000000 0.0000000000000000 0.0000000000000000
37 0.0000000000000000 0.0000000000000000 0.0000000000000000
38 0.0000000000000000 0.0000000000000000 0.0000000000000000
39 0.0000000000000000 0.0000000000000000 0.0000000000000000
40 0.0000000000000000 0.0000000000000000 0.0000000000000000
41 0.0000000000000000 0.0000000000000000 0.0000000000000000
42 0.0000000000000000 0.0000000000000000 0.0000000000000000
43 0.0000000000000000 0.0000000000000000 0.0000000000000000
44 0.0000000000000000 0.0000000000000000 0.0000000000000000
45 0.0000000000000000 0.0000000000000000 0.0000000000000000
46 0.0000000000000000 0.0000000000000000 0.0000000000000000
47 0.0000000000000000 0.0000000000000000 0.0000000000000000
48 0.0000000000000000 0.0000000000000000 0.0000000000000000
49 0.0000000000000000 0.0000000000000000 0.0000000000000000
50 0.0000000000000000 0.0000000000000000 0.0000000000000000
51 0.0000000000000000 0.0000000000000000 0.0000000000000000
$EndNodeData

View file

@ -0,0 +1,247 @@
$MeshFormat
2.2 0 8
$EndMeshFormat
$NodeData
1
"Electron density (m^-3)"
1
1.0000000000000000E-008
3
1000
1
51
1 6.926728E+018
2 4.338487E+014
3 1.201033E+018
4 1.505055E+017
5 7.986344E+016
6 4.403063E+016
7 3.228092E+016
8 2.775434E+016
9 2.012556E+016
10 1.435034E+016
11 1.012383E+016
12 7.086783E+015
13 6.226656E+015
14 6.362452E+015
15 6.664819E+015
16 5.946772E+015
17 4.875712E+015
18 4.419041E+015
19 4.092162E+015
20 3.771210E+015
21 3.309077E+015
22 3.071732E+015
23 2.732966E+015
24 2.484798E+015
25 2.313440E+015
26 2.142727E+015
27 1.892093E+015
28 1.776155E+015
29 1.672823E+015
30 1.493186E+015
31 1.407724E+015
32 1.403006E+015
33 1.262097E+015
34 1.182309E+015
35 1.134819E+015
36 1.047928E+015
37 9.865986E+014
38 9.544875E+014
39 8.680888E+014
40 8.440321E+014
41 7.961164E+014
42 7.174005E+014
43 6.589875E+014
44 6.158189E+014
45 6.203645E+014
46 5.604593E+014
47 5.677392E+014
48 5.696204E+014
49 5.003183E+014
50 4.735248E+014
51 4.522101E+014
$EndNodeData
$NodeData
1
"Electron velocity (m s^-1)"
1
1.0000000000000000E-008
3
1000
3
51
1 6.755624E+003 2.474633E+002 2.390054E+002
2 2.054354E+005 -5.163376E+003 -9.846850E+003
3 1.851665E+004 5.899050E+001 2.941097E+002
4 5.874981E+003 -2.544408E+001 3.338329E+002
5 5.745724E+003 -2.992103E+002 -1.562428E+002
6 9.042765E+003 5.755962E+002 -7.509840E+002
7 1.189973E+004 1.251817E+003 -1.278355E+003
8 1.171915E+004 2.483533E+002 -1.798212E+003
9 1.539467E+004 -4.794788E+002 -1.578985E+002
10 1.749366E+004 -2.133315E+003 2.884641E+003
11 3.073569E+004 -1.813877E+003 2.015157E+002
12 4.338021E+004 4.080113E+003 -1.903833E+002
13 4.108799E+004 1.245391E+003 -1.163594E+003
14 3.541737E+004 2.090477E+003 2.419701E+002
15 3.427041E+004 2.013895E+003 -2.234277E+003
16 4.230007E+004 -2.398689E+003 1.371723E+002
17 4.749322E+004 1.597246E+003 -3.628288E+003
18 4.997666E+004 1.160104E+003 -8.401792E+003
19 5.195058E+004 -8.364161E+002 -1.802101E+003
20 5.591312E+004 -3.939060E+003 2.388800E+003
21 6.333098E+004 -6.945454E+002 2.258530E+002
22 6.525072E+004 -3.210740E+003 1.801044E+003
23 7.399242E+004 -1.986068E+003 2.887907E+002
24 7.895156E+004 3.349907E+003 -2.339613E+003
25 8.100325E+004 3.214988E+003 9.916669E+002
26 8.973728E+004 7.594115E+002 -3.305504E+003
27 9.309609E+004 1.417677E+003 -1.436153E+002
28 9.726401E+004 -3.383752E+003 2.393699E+003
29 1.004495E+005 -2.281003E+003 3.209291E+003
30 1.107262E+005 -1.404186E+003 7.470032E+002
31 1.202574E+005 -2.703508E+003 4.930279E+003
32 1.224420E+005 -3.973287E+002 -7.416719E+003
33 1.312544E+005 -5.905453E+003 -9.093921E+002
34 1.347764E+005 -3.491930E+003 -4.010455E+003
35 1.376741E+005 1.357565E+003 2.096097E+003
36 1.419675E+005 1.982650E+003 6.535850E+002
37 1.437177E+005 9.484341E+002 -1.187650E+004
38 1.468643E+005 2.894937E+003 -5.271059E+002
39 1.559764E+005 2.200720E+003 -6.344991E+003
40 1.599147E+005 -2.663692E+003 5.719357E+003
41 1.567097E+005 -7.359025E+000 -4.033384E+003
42 1.690415E+005 -6.956095E+002 3.022908E+003
43 1.814588E+005 -9.615306E+001 6.917210E+003
44 1.863152E+005 -7.990127E+002 1.523384E+004
45 1.886311E+005 4.934329E+003 4.438134E+002
46 1.938832E+005 -3.121504E+002 -3.067774E+003
47 1.900862E+005 -1.160531E+003 -4.620352E+003
48 1.923784E+005 -3.210640E+003 2.151672E+003
49 2.043583E+005 2.368308E+003 1.749902E+002
50 2.054704E+005 4.364283E+003 1.393513E+003
51 2.125792E+005 -4.749345E+003 -2.018864E+003
$EndNodeData
$NodeData
1
"Electron pressure (Pa)"
1
1.0000000000000000E-008
3
1000
1
51
1 2.252567E-001
2 7.546705E-006
3 4.045878E-002
4 5.008147E-003
5 2.551765E-003
6 1.488666E-003
7 1.021884E-003
8 7.937038E-004
9 5.788516E-004
10 4.290930E-004
11 3.263161E-004
12 2.445814E-004
13 2.004129E-004
14 1.799644E-004
15 1.711394E-004
16 1.505493E-004
17 1.238843E-004
18 1.124456E-004
19 9.851639E-005
20 8.889251E-005
21 7.626551E-005
22 6.799265E-005
23 5.942787E-005
24 5.395244E-005
25 4.992406E-005
26 4.474483E-005
27 4.103992E-005
28 3.826699E-005
29 3.404848E-005
30 3.038250E-005
31 2.882086E-005
32 2.844685E-005
33 2.486692E-005
34 2.305064E-005
35 2.204186E-005
36 2.045655E-005
37 1.912248E-005
38 1.817270E-005
39 1.663693E-005
40 1.616886E-005
41 1.470869E-005
42 1.249547E-005
43 1.189200E-005
44 1.112158E-005
45 1.148139E-005
46 1.015354E-005
47 9.973794E-006
48 9.467041E-006
49 8.605269E-006
50 8.068109E-006
51 8.336811E-006
$EndNodeData
$NodeData
1
"Electron temperature (K)"
1
1.0000000000000000E-008
3
1000
1
51
1 2.355409E+003
2 1.259899E+003
3 2.439915E+003
4 2.410136E+003
5 2.314246E+003
6 2.448833E+003
7 2.292833E+003
8 2.071307E+003
9 2.083225E+003
10 2.165738E+003
11 2.334589E+003
12 2.499719E+003
13 2.331244E+003
14 2.048702E+003
15 1.859852E+003
16 1.833641E+003
17 1.840328E+003
18 1.843025E+003
19 1.743703E+003
20 1.707266E+003
21 1.669315E+003
22 1.603229E+003
23 1.574972E+003
24 1.572667E+003
25 1.563034E+003
26 1.512492E+003
27 1.571017E+003
28 1.560488E+003
29 1.474228E+003
30 1.473759E+003
31 1.482881E+003
32 1.468559E+003
33 1.427073E+003
34 1.412111E+003
35 1.406820E+003
36 1.413897E+003
37 1.403850E+003
38 1.379006E+003
39 1.388116E+003
40 1.387513E+003
41 1.338179E+003
42 1.261559E+003
43 1.307058E+003
44 1.308069E+003
45 1.340493E+003
46 1.312170E+003
47 1.272414E+003
48 1.203776E+003
49 1.245762E+003
50 1.234087E+003
51 1.335293E+003
$EndNodeData

View file

@ -0,0 +1,7 @@
Alphie grid system.
The file mesh.geo can be easily modified to generate new grids to test different configurations.
Two cases are provided:
Clasical case in which ions and electrons are input from the ionization chamber.
Ionization case in which ion-electron pairs are generated due to the influx of electrons using the Ionization boundary condition.

View file

@ -8,7 +8,8 @@
"folder": "base_case"
},
"geometry": {
"type": "2DCyl",
"dimension": 2,
"type": "Cyl",
"meshType": "gmsh2",
"meshFile": "mesh.msh"
},
@ -45,14 +46,15 @@
"boundaryEM": [
{"name": "Extraction Grid", "type": "dirichlet", "potential": -150.0, "physicalSurface": 4},
{"name": "Acceleration Grid", "type": "dirichlet", "potential": -600.0, "physicalSurface": 5},
{"name": "Ionization Chamber", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1}
{"name": "Ionization Chamber", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1},
{"name": "Infinite", "type": "dirichlet", "potential": -600.0, "physicalSurface": 2}
],
"inject": [
{"name": "Ionization Argon+", "species": "Argon+", "flow": 27.0e-6, "units": "A", "v": 322.0, "T": [ 500.0, 500.0, 500.0],
{"name": "Ionization Argon+", "species": "Argon+", "flow": 1.0e-5, "units": "A", "v": 2500.0, "T": [ 500.0, 500.0, 500.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1},
{"name": "Ionization Electron", "species": "Electron", "flow": 27.0e-6, "units": "A", "v": 87000.0, "T": [ 500.0, 500.0, 500.0],
{"name": "Ionization Electron", "species": "Electron", "flow": 1.0e-5, "units": "A", "v": 87000.0, "T": [30000.0, 30000.0, 30000.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1},
{"name": "Cathode Electron", "species": "Electron", "flow": 9.0e-5, "units": "A", "v": 87000.0, "T": [2500.0, 2500.0, 2500.0],
{"name": "Cathode Electron", "species": "Electron", "flow": 1.0e-4, "units": "A", "v": 87000.0, "T": [30000.0, 30000.0, 30000.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2}
],
"reference": {
@ -61,10 +63,10 @@
"temperature": 2500.0,
"radius": 1.88e-10
},
"case": {
"solver": {
"tau": [1.0e-9, 1.0e-11],
"time": 1.0e-6,
"pusher": ["2DCylCharged", "2DCylCharged"],
"finalTime": 1.0e-6,
"pusher": ["Electrostatic", "Electrostatic"],
"WeightingScheme": "Volume",
"EMSolver": "Electrostatic"
},

View file

@ -5,10 +5,11 @@
"cpuTime": false,
"numColl": false,
"EMField": true,
"folder": "ionization_0.10mA"
"folder": "ionization_0.10mA"
},
"geometry": {
"type": "2DCyl",
"dimension": 2,
"type": "Cyl",
"meshType": "gmsh2",
"meshFile": "mesh.msh"
},
@ -19,7 +20,7 @@
"boundary": [
{"name": "Ionization Chanber", "physicalSurface": 1, "bTypes": [
{"type": "transparent"},
{"type": "ionization", "neutral": {"ion": "Argon+", "mass": 6.633e-26, "density": 1.0e17, "velocity": [323, 0, 0], "temperature": 300},
{"type": "ionization", "neutral": {"ion": "Argon+", "mass": 6.633e-26, "density": 5.0e16, "velocity": [2500, 0, 0], "temperature": 300},
"effectiveTime": 5.0e-6,"energyThreshold": 15.76, "crossSection": "./data/collisions/IO_e-Ar.dat"}
]},
{"name": "Vacuum Chamber", "physicalSurface": 2, "bTypes": [
@ -46,10 +47,11 @@
"boundaryEM": [
{"name": "Extraction Grid", "type": "dirichlet", "potential": -150.0, "physicalSurface": 4},
{"name": "Acceleration Grid", "type": "dirichlet", "potential": -600.0, "physicalSurface": 5},
{"name": "Ionization Chamber", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1}
{"name": "Ionization Chamber", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1},
{"name": "Infinite", "type": "dirichlet", "potential": -600.0, "physicalSurface": 2}
],
"inject": [
{"name": "Cathode Electron", "species": "Electron", "flow": 1.0e-4, "units": "A", "v": 87000.0, "T": [2500.0, 2500.0, 2500.0],
{"name": "Cathode Electron", "species": "Electron", "flow": 1.0e-4, "units": "A", "v": 87000.0, "T": [30000.0, 30000.0, 30000.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2}
],
"reference": {
@ -58,10 +60,10 @@
"temperature": 2500.0,
"radius": 1.88e-10
},
"case": {
"solver": {
"tau": [1.0e-9, 1.0e-11],
"time": 1.0e-6,
"pusher": ["2DCylCharged", "2DCylCharged"],
"finalTime": 1.0e-6,
"pusher": ["Electrostatic", "Electrostatic"],
"WeightingScheme": "Volume",
"EMSolver": "Electrostatic"
},

View file

@ -1,11 +1,11 @@
zg1 = 0.0025;
tg1 = 0.0004;
zg1 = 0.0020;
tg1 = 0.0003;
rg1 = 0.0005;
dg = 0.0025;
zg2 = zg1+tg1+dg;
dg = 0.0020;
zg2 = zg1 + tg1 + dg;
tg2 = tg1;
rg2 = rg1;
zEnd = 0.0042;
zEnd = 0.0050;
Lz = zg2 + tg2 + zEnd;
Lr = rg1 + 0.0001;

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

View file

@ -7,7 +7,8 @@
"folder": "CX_case"
},
"geometry": {
"type": "2DCyl",
"dimension": 2,
"type": "Cyl",
"meshType": "gmsh2",
"meshFile": "mesh.msh"
},
@ -41,10 +42,10 @@
"temperature": 300.0,
"radius": 1.88e-10
},
"case": {
"solver": {
"tau": [1.0e-6, 1.0e-6],
"time": 4.0e-3,
"pusher": ["2DCylNeutral", "2DCylNeutral"],
"finalTime": 4.0e-3,
"pusher": ["Neutral", "Neutral"],
"WeightingScheme": "Volume"
},
"interactions": {

View file

@ -7,7 +7,8 @@
"folder": "Elastic_case"
},
"geometry": {
"type": "2DCyl",
"dimension": 2,
"type": "Cyl",
"meshType": "gmsh2",
"meshFile": "mesh.msh"
},
@ -41,10 +42,10 @@
"temperature": 300.0,
"radius": 1.88e-10
},
"case": {
"solver": {
"tau": [1.0e-6, 1.0e-6],
"time": 4.0e-3,
"pusher": ["2DCylNeutral", "2DCylNeutral"],
"finalTime": 4.0e-3,
"pusher": ["Neutral", "Neutral"],
"WeightingScheme": "Volume"
},
"interactions": {

View file

@ -7,7 +7,8 @@
"folder": "Nocoll_case"
},
"geometry": {
"type": "2DCyl",
"dimension": 2,
"type": "Cyl",
"meshType": "gmsh2",
"meshFile": "mesh.msh"
},
@ -41,10 +42,10 @@
"temperature": 300.0,
"radius": 1.88e-10
},
"case": {
"solver": {
"tau": [1.0e-6, 1.0e-6],
"time": 4.0e-3,
"pusher": ["2DCylNeutral", "2DCylNeutral"],
"finalTime": 4.0e-3,
"pusher": ["Neutral", "Neutral"],
"WeightingScheme": "Volume"
},
"parallel": {

7
runs/cylFlow/README.txt Normal file
View file

@ -0,0 +1,7 @@
Neutral flow around a cylinder with elastic collisions.
Case with the same mesh for particle weighting and collisions and another case using the dual mesh mode, with a more refined mesh around the cylinder.
This demostrates the capacity of fpakc to work with independent meshes for particles and collisions.
CPU time is outputed.
Gnuplot scripts for plotting the CPU time are provided.

View file

@ -6,7 +6,8 @@
"numColl": true
},
"geometry": {
"type": "2DCyl",
"dimension": 2,
"type": "Cyl",
"meshType": "gmsh2",
"meshFile": "meshSingle.msh"
},
@ -40,11 +41,10 @@
"temperature": 300.0,
"radius": 1.88e-10
},
"case": {
"solver": {
"tau": [5.0e-7],
"time": 1.0e-3,
"pusher": ["2DCylNeutral"],
"WeightingScheme": "Volume"
"finalTime": 5.0e-4,
"pusher": ["Neutral"],
},
"interactions": {
"folderCollisions": "./data/collisions/",

View file

@ -6,9 +6,11 @@
"numColl": true
},
"geometry": {
"type": "2DCyl",
"dimension": 2,
"type": "Cyl",
"meshType": "gmsh2",
"meshFile": "mesh.msh"
"meshFile": "mesh.msh",
"meshCollisions": "meshColl.msh"
},
"species": [
{"name": "Argon", "type": "neutral", "mass": 6.633e-26, "weight": 5.0e8}
@ -40,15 +42,14 @@
"temperature": 300.0,
"radius": 1.88e-10
},
"case": {
"solver": {
"tau": [5.0e-7],
"time": 1.0e-3,
"pusher": ["2DCylNeutral"],
"finalTime": 5.0e-4,
"pusher": ["Neutral"],
"WeightingScheme": "Volume"
},
"interactions": {
"folderCollisions": "./data/collisions/",
"meshCollisions": "meshColl.msh",
"collisions": [
{"species_i": "Argon", "species_j": "Argon",
"cTypes": [

View file

@ -1,629 +0,0 @@
General.AxesFormatX = "%.3g";
General.AxesFormatY = "%.3g";
General.AxesFormatZ = "%.3g";
General.AxesLabelX = "z (m)";
General.AxesLabelY = "r (m)";
General.AxesLabelZ = "";
General.BackgroundImageFileName = "";
General.BuildOptions = " 64Bit ALGLIB Bamg Blas Blossom DIntegration Dlopen DomHex Fltk GMP Gmm[system] Hxt Hxt3D Jpeg Kbipack Lapack LinuxJoystick MathEx Mesh Metis[system] Mmg3d Mpeg NativeFileChooser Netgen ONELAB ONELABMetamodel OpenCASCADE OpenCASCADE-CAF OpenGL OpenMP OptHom Parser Plugins Png Post QuadTri Solver TetGen/BR Zlib";
General.DefaultFileName = "untitled.geo";
General.Display = "";
General.ErrorFileName = ".gmsh-errors";
General.ExecutableFileName = "/usr/bin/gmsh";
General.FileName = "/home/jorge/PPartiC/runs/cylFlow/mesh.msh";
General.FltkTheme = "";
General.GraphicsFont = "Helvetica";
General.GraphicsFontEngine = "Native";
General.GraphicsFontTitle = "Helvetica";
General.OptionsFileName = ".gmsh-options";
General.RecentFile0 = "/home/jorge/PPartiC/runs/cylFlow/mesh.msh";
General.RecentFile1 = "untitled.geo";
General.RecentFile2 = "mesh/Neutral_Expansion_Div.geo";
General.RecentFile3 = "mesh/Neutral_Expansion_Div.msh";
General.RecentFile4 = "Neutral_Expansion_Div.msh";
General.RecentFile5 = "Neutral_Expansion.geo";
General.RecentFile6 = "Neutral_Expansion_Div.geo";
General.RecentFile7 = "/home/jorge/Dropbox/UPMPlasmaLab/Post-Doc/Codes/PICCIL2D/PIC-FEM/Neutral_Expansion_Div.msh";
General.RecentFile8 = "Neutral_Expansion.msh";
General.RecentFile9 = "/home/jorge/Dropbox/UPMPlasmaLab/Post-Doc/Codes/PICCIL2D/PIC-FEM/Neutral_Expansion.msh";
General.TextEditor = "gedit '%s'";
General.TmpFileName = ".gmsh-tmp";
General.Version = "4.4.1";
General.WatchFilePattern = "";
General.AlphaBlending = 1;
General.Antialiasing = 0;
General.ArrowHeadRadius = 0.12;
General.ArrowStemLength = 0.5600000000000001;
General.ArrowStemRadius = 0.02;
General.Axes = 1;
General.AxesMikado = 0;
General.AxesAutoPosition = 1;
General.AxesForceValue = 0;
General.AxesMaxX = 1;
General.AxesMaxY = 1;
General.AxesMaxZ = 1;
General.AxesMinX = 0;
General.AxesMinY = 0;
General.AxesMinZ = 0;
General.AxesTicsX = 8;
General.AxesTicsY = 4;
General.AxesTicsZ = 5;
General.AxesValueMaxX = 1;
General.AxesValueMaxY = 1;
General.AxesValueMaxZ = 1;
General.AxesValueMinX = 0;
General.AxesValueMinY = 0;
General.AxesValueMinZ = 0;
General.BackgroundGradient = 1;
General.BackgroundImage3D = 0;
General.BackgroundImagePage = 0;
General.BackgroundImagePositionX = 0;
General.BackgroundImagePositionY = 0;
General.BackgroundImageWidth = -1;
General.BackgroundImageHeight = -1;
General.BoundingBoxSize = 0.07615773105863909;
General.Camera = 0;
General.CameraAperture = 40;
General.CameraEyeSeparationRatio = 1.5;
General.CameraFocalLengthRatio = 1;
General.Clip0A = 1;
General.Clip0B = 0;
General.Clip0C = 0;
General.Clip0D = 0;
General.Clip1A = 0;
General.Clip1B = 1;
General.Clip1C = 0;
General.Clip1D = 0;
General.Clip2A = 0;
General.Clip2B = 0;
General.Clip2C = 1;
General.Clip2D = 0;
General.Clip3A = -1;
General.Clip3B = 0;
General.Clip3C = 0;
General.Clip3D = 1;
General.Clip4A = 0;
General.Clip4B = -1;
General.Clip4C = 0;
General.Clip4D = 1;
General.Clip5A = 0;
General.Clip5B = 0;
General.Clip5C = -1;
General.Clip5D = 1;
General.ClipFactor = 5;
General.ClipOnlyDrawIntersectingVolume = 0;
General.ClipOnlyVolume = 0;
General.ClipPositionX = 650;
General.ClipPositionY = 150;
General.ClipWholeElements = 0;
General.ColorScheme = 1;
General.ConfirmOverwrite = 1;
General.ContextPositionX = 235;
General.ContextPositionY = 962;
General.DetachedMenu = 0;
General.DisplayBorderFactor = 0.2;
General.DoubleBuffer = 1;
General.DrawBoundingBoxes = 0;
General.ExpertMode = 0;
General.ExtraPositionX = 650;
General.ExtraPositionY = 350;
General.ExtraHeight = 100;
General.ExtraWidth = 100;
General.FastRedraw = 0;
General.FieldPositionX = 650;
General.FieldPositionY = 550;
General.FieldHeight = 488;
General.FieldWidth = 651;
General.FileChooserPositionX = 200;
General.FileChooserPositionY = 200;
General.FltkColorScheme = 0;
General.FontSize = -1;
General.GraphicsFontSize = 15;
General.GraphicsFontSizeTitle = 18;
General.GraphicsHeight = 1003;
General.GraphicsPositionX = 274;
General.GraphicsPositionY = 263;
General.GraphicsWidth = 1920;
General.HighOrderToolsPositionX = 650;
General.HighOrderToolsPositionY = 150;
General.HighResolutionGraphics = 1;
General.HighResolutionPointSizeFactor = 2;
General.InitialModule = 0;
General.InputScrolling = 1;
General.Light0 = 1;
General.Light0X = 0.65;
General.Light0Y = 0.65;
General.Light0Z = 1;
General.Light0W = 0;
General.Light1 = 0;
General.Light1X = 0.5;
General.Light1Y = 0.3;
General.Light1Z = 1;
General.Light1W = 0;
General.Light2 = 0;
General.Light2X = 0.5;
General.Light2Y = 0.3;
General.Light2Z = 1;
General.Light2W = 0;
General.Light3 = 0;
General.Light3X = 0.5;
General.Light3Y = 0.3;
General.Light3Z = 1;
General.Light3W = 0;
General.Light4 = 0;
General.Light4X = 0.5;
General.Light4Y = 0.3;
General.Light4Z = 1;
General.Light4W = 0;
General.Light5 = 0;
General.Light5X = 0.5;
General.Light5Y = 0.3;
General.Light5Z = 1;
General.Light5W = 0;
General.LineWidth = 1;
General.ManipulatorPositionX = 650;
General.ManipulatorPositionY = 150;
General.MaxX = 0.07000000000000001;
General.MaxY = 0.03;
General.MaxZ = 0;
General.MenuWidth = 219;
General.MenuHeight = 200;
General.MenuPositionX = 400;
General.MenuPositionY = 400;
General.MessageFontSize = -1;
General.MessageHeight = 300;
General.MinX = 0;
General.MinY = 0;
General.MinZ = 0;
General.MouseHoverMeshes = 0;
General.MouseSelection = 1;
General.MouseInvertZoom = 0;
General.NonModalWindows = 1;
General.NoPopup = 0;
General.NumThreads = 1;
General.OptionsPositionX = 827;
General.OptionsPositionY = 541;
General.Orthographic = 1;
General.PluginPositionX = 58;
General.PluginPositionY = 658;
General.PluginHeight = 488;
General.PluginWidth = 708;
General.PointSize = 3;
General.PolygonOffsetAlwaysOn = 0;
General.PolygonOffsetFactor = 1;
General.PolygonOffsetUnits = 1;
General.ProgressMeterStep = 20;
General.QuadricSubdivisions = 6;
General.RotationX = -0;
General.RotationY = 0;
General.RotationZ = -0;
General.RotationCenterGravity = 1;
General.RotationCenterX = 0;
General.RotationCenterY = 0;
General.RotationCenterZ = 0;
General.SaveOptions = 0;
General.SaveSession = 1;
General.ScaleX = 1;
General.ScaleY = 1;
General.ScaleZ = 1;
General.Shininess = 0.4;
General.ShininessExponent = 40;
General.ShowModuleMenu = 1;
General.ShowOptionsOnStartup = 0;
General.ShowMessagesOnStartup = 0;
General.SmallAxes = 1;
General.SmallAxesPositionX = -60;
General.SmallAxesPositionY = -40;
General.SmallAxesSize = 30;
General.StatisticsPositionX = 650;
General.StatisticsPositionY = 150;
General.Stereo = 0;
General.SystemMenuBar = 1;
General.Terminal = 0;
General.Tooltips = 1;
General.Trackball = 1;
General.TrackballHyperbolicSheet = 1;
General.TrackballQuaternion0 = 0;
General.TrackballQuaternion1 = 0;
General.TrackballQuaternion2 = 0;
General.TrackballQuaternion3 = 1;
General.TranslationX = 0;
General.TranslationY = 0;
General.TranslationZ = 0;
General.VectorType = 4;
General.Verbosity = 5;
General.VisibilityPositionX = 1118;
General.VisibilityPositionY = 464;
General.ZoomFactor = 4;
General.Color.Background = {255,255,255};
General.Color.BackgroundGradient = {208,215,255};
General.Color.Foreground = {85,85,85};
General.Color.Text = {0,0,0};
General.Color.Axes = {0,0,0};
General.Color.SmallAxes = {0,0,0};
General.Color.AmbientLight = {25,25,25};
General.Color.DiffuseLight = {255,255,255};
General.Color.SpecularLight = {255,255,255};
Geometry.DoubleClickedPointCommand = "";
Geometry.DoubleClickedLineCommand = "";
Geometry.DoubleClickedSurfaceCommand = "";
Geometry.DoubleClickedVolumeCommand = "";
Geometry.OCCTargetUnit = "";
Geometry.AutoCoherence = 1;
Geometry.Clip = 0;
Geometry.CopyMeshingMethod = 0;
Geometry.DoubleClickedEntityTag = 0;
Geometry.ExactExtrusion = 1;
Geometry.ExtrudeReturnLateralEntities = 1;
Geometry.ExtrudeSplinePoints = 5;
Geometry.HighlightOrphans = 0;
Geometry.LabelType = 0;
Geometry.Light = 1;
Geometry.LightTwoSide = 1;
Geometry.Lines = 1;
Geometry.LineNumbers = 0;
Geometry.LineSelectWidth = 3;
Geometry.LineType = 0;
Geometry.LineWidth = 2;
Geometry.MatchGeomAndMesh = 0;
Geometry.MatchMeshScaleFactor = 1;
Geometry.MatchMeshTolerance = 1e-06;
Geometry.Normals = 0;
Geometry.NumSubEdges = 40;
Geometry.OCCAutoFix = 1;
Geometry.OCCBooleanPreserveNumbering = 1;
Geometry.OCCDisableSTL = 0;
Geometry.OCCFixDegenerated = 0;
Geometry.OCCFixSmallEdges = 0;
Geometry.OCCFixSmallFaces = 0;
Geometry.OCCImportLabels = 1;
Geometry.OCCParallel = 0;
Geometry.OCCScaling = 1;
Geometry.OCCSewFaces = 0;
Geometry.OffsetX = 0;
Geometry.OffsetY = 0;
Geometry.OffsetZ = 0;
Geometry.OldCircle = 0;
Geometry.OldRuledSurface = 0;
Geometry.OldNewReg = 1;
Geometry.Points = 1;
Geometry.PointNumbers = 0;
Geometry.PointSelectSize = 6;
Geometry.PointSize = 4;
Geometry.PointType = 0;
Geometry.ReparamOnFaceRobust = 0;
Geometry.ScalingFactor = 1;
Geometry.OrientedPhysicals = 1;
Geometry.SnapX = 0.1;
Geometry.SnapY = 0.1;
Geometry.SnapZ = 0.1;
Geometry.Surfaces = 0;
Geometry.SurfaceNumbers = 0;
Geometry.SurfaceType = 0;
Geometry.Tangents = 0;
Geometry.Tolerance = 1e-08;
Geometry.ToleranceBoolean = 0;
Geometry.Transform = 0;
Geometry.TransformXX = 1;
Geometry.TransformXY = 0;
Geometry.TransformXZ = 0;
Geometry.TransformYX = 0;
Geometry.TransformYY = 1;
Geometry.TransformYZ = 0;
Geometry.TransformZX = 0;
Geometry.TransformZY = 0;
Geometry.TransformZZ = 1;
Geometry.Volumes = 0;
Geometry.VolumeNumbers = 0;
Geometry.Color.Points = {90,90,90};
Geometry.Color.Lines = {0,0,255};
Geometry.Color.Surfaces = {128,128,128};
Geometry.Color.Volumes = {255,255,0};
Geometry.Color.Selection = {255,0,0};
Geometry.Color.HighlightZero = {255,0,0};
Geometry.Color.HighlightOne = {255,150,0};
Geometry.Color.HighlightTwo = {255,255,0};
Geometry.Color.Tangents = {255,255,0};
Geometry.Color.Normals = {255,0,0};
Geometry.Color.Projection = {0,255,0};
Mesh.Algorithm = 2;
Mesh.Algorithm3D = 1;
Mesh.AngleSmoothNormals = 30;
Mesh.AngleToleranceFacetOverlap = 0.1;
Mesh.AnisoMax = 9.999999999999999e+32;
Mesh.AllowSwapAngle = 10;
Mesh.BdfFieldFormat = 1;
Mesh.Binary = 0;
Mesh.BoundaryLayerFanPoints = 5;
Mesh.CgnsImportOrder = 1;
Mesh.CgnsConstructTopology = 0;
Mesh.CharacteristicLengthExtendFromBoundary = 1;
Mesh.CharacteristicLengthFactor = 1;
Mesh.CharacteristicLengthMin = 0;
Mesh.CharacteristicLengthMax = 1e+22;
Mesh.CharacteristicLengthFromCurvature = 0;
Mesh.CharacteristicLengthFromPoints = 1;
Mesh.Clip = 0;
Mesh.ColorCarousel = 0;
Mesh.CompoundClassify = 1;
Mesh.CompoundCharacteristicLengthFactor = 0.5;
Mesh.CpuTime = 0;
Mesh.DrawSkinOnly = 0;
Mesh.Dual = 0;
Mesh.ElementOrder = 1;
Mesh.Explode = 1;
Mesh.FlexibleTransfinite = 0;
Mesh.NewtonConvergenceTestXYZ = 0;
Mesh.Format = 10;
Mesh.Hexahedra = 1;
Mesh.HighOrderIterMax = 100;
Mesh.HighOrderNumLayers = 6;
Mesh.HighOrderOptimize = 0;
Mesh.HighOrderPassMax = 25;
Mesh.HighOrderPeriodic = 0;
Mesh.HighOrderPoissonRatio = 0.33;
Mesh.HighOrderPrimSurfMesh = 0;
Mesh.HighOrderDistCAD = 0;
Mesh.HighOrderThresholdMin = 0.1;
Mesh.HighOrderThresholdMax = 2;
Mesh.LabelSampling = 1;
Mesh.LabelType = 0;
Mesh.LcIntegrationPrecision = 1e-09;
Mesh.Light = 1;
Mesh.LightLines = 2;
Mesh.LightTwoSide = 1;
Mesh.Lines = 1;
Mesh.LineNumbers = 0;
Mesh.LineWidth = 1;
Mesh.MaxNumThreads1D = 0;
Mesh.MaxNumThreads2D = 0;
Mesh.MaxNumThreads3D = 0;
Mesh.MeshOnlyVisible = 0;
Mesh.MetisAlgorithm = 1;
Mesh.MetisEdgeMatching = 2;
Mesh.MetisMaxLoadImbalance = -1;
Mesh.MetisObjective = 1;
Mesh.MetisMinConn = -1;
Mesh.MetisRefinementAlgorithm = 2;
Mesh.MinimumCirclePoints = 7;
Mesh.MinimumCurvePoints = 3;
Mesh.MshFileVersion = 4.1;
Mesh.MedFileMinorVersion = -1;
Mesh.MedImportGroupsOfNodes = 0;
Mesh.MedSingleModel = 0;
Mesh.PartitionHexWeight = -1;
Mesh.PartitionLineWeight = -1;
Mesh.PartitionPrismWeight = -1;
Mesh.PartitionPyramidWeight = -1;
Mesh.PartitionQuadWeight = -1;
Mesh.PartitionTrihedronWeight = 0;
Mesh.PartitionTetWeight = -1;
Mesh.PartitionTriWeight = -1;
Mesh.PartitionCreateTopology = 1;
Mesh.PartitionCreatePhysicals = 1;
Mesh.PartitionCreateGhostCells = 0;
Mesh.PartitionSplitMeshFiles = 0;
Mesh.PartitionTopologyFile = 0;
Mesh.PartitionOldStyleMsh2 = 1;
Mesh.NbHexahedra = 0;
Mesh.NbNodes = 27248;
Mesh.NbPartitions = 0;
Mesh.NbPrisms = 0;
Mesh.NbPyramids = 0;
Mesh.NbTrihedra = 0;
Mesh.NbQuadrangles = 26877;
Mesh.NbTetrahedra = 0;
Mesh.NbTriangles = 0;
Mesh.Normals = 0;
Mesh.NumSubEdges = 2;
Mesh.Optimize = 1;
Mesh.OptimizeThreshold = 0.3;
Mesh.OptimizeNetgen = 0;
Mesh.Points = 0;
Mesh.PointNumbers = 0;
Mesh.PointSize = 4;
Mesh.PointType = 0;
Mesh.Prisms = 1;
Mesh.Pyramids = 1;
Mesh.Trihedra = 1;
Mesh.Quadrangles = 1;
Mesh.QualityInf = 0;
Mesh.QualitySup = 0;
Mesh.QualityType = 2;
Mesh.RadiusInf = 0;
Mesh.RadiusSup = 0;
Mesh.RandomFactor = 1e-09;
Mesh.RandomFactor3D = 1e-12;
Mesh.RandomSeed = 1;
Mesh.PreserveNumberingMsh2 = 0;
Mesh.IgnorePeriodicity = 0;
Mesh.RecombinationAlgorithm = 1;
Mesh.RecombineAll = 0;
Mesh.RecombineOptimizeTopology = 5;
Mesh.Recombine3DAll = 0;
Mesh.Recombine3DLevel = 0;
Mesh.Recombine3DConformity = 0;
Mesh.RefineSteps = 10;
Mesh.Renumber = 1;
Mesh.SaveAll = 0;
Mesh.SaveElementTagType = 1;
Mesh.SaveTopology = 0;
Mesh.SaveParametric = 0;
Mesh.SaveGroupsOfNodes = 0;
Mesh.ScalingFactor = 1;
Mesh.SecondOrderExperimental = 0;
Mesh.SecondOrderIncomplete = 0;
Mesh.SecondOrderLinear = 0;
Mesh.Smoothing = 1;
Mesh.SmoothCrossField = 0;
Mesh.CrossFieldClosestPoint = 1;
Mesh.SmoothNormals = 0;
Mesh.SmoothRatio = 1.8;
Mesh.StlOneSolidPerSurface = 0;
Mesh.StlRemoveDuplicateTriangles = 0;
Mesh.SubdivisionAlgorithm = 0;
Mesh.SurfaceEdges = 0;
Mesh.SurfaceFaces = 0;
Mesh.SurfaceNumbers = 0;
Mesh.SwitchElementTags = 0;
Mesh.Tangents = 0;
Mesh.Tetrahedra = 1;
Mesh.ToleranceEdgeLength = 0;
Mesh.ToleranceInitialDelaunay = 1e-08;
Mesh.Triangles = 1;
Mesh.UnvStrictFormat = 1;
Mesh.VolumeEdges = 1;
Mesh.VolumeFaces = 0;
Mesh.VolumeNumbers = 0;
Mesh.Voronoi = 0;
Mesh.ZoneDefinition = 0;
Mesh.Color.Points = {0,0,255};
Mesh.Color.PointsSup = {255,0,255};
Mesh.Color.Lines = {0,0,0};
Mesh.Color.Triangles = {160,150,255};
Mesh.Color.Quadrangles = {130,120,225};
Mesh.Color.Tetrahedra = {160,150,255};
Mesh.Color.Hexahedra = {130,120,225};
Mesh.Color.Prisms = {232,210,23};
Mesh.Color.Pyramids = {217,113,38};
Mesh.Color.Trihedra = {20,255,0};
Mesh.Color.Tangents = {255,255,0};
Mesh.Color.Normals = {255,0,0};
Mesh.Color.Zero = {255,120,0};
Mesh.Color.One = {0,255,132};
Mesh.Color.Two = {255,160,0};
Mesh.Color.Three = {0,255,192};
Mesh.Color.Four = {255,200,0};
Mesh.Color.Five = {0,216,255};
Mesh.Color.Six = {255,240,0};
Mesh.Color.Seven = {0,176,255};
Mesh.Color.Eight = {228,255,0};
Mesh.Color.Nine = {0,116,255};
Mesh.Color.Ten = {188,255,0};
Mesh.Color.Eleven = {0,76,255};
Mesh.Color.Twelve = {148,255,0};
Mesh.Color.Thirteen = {24,0,255};
Mesh.Color.Fourteen = {108,255,0};
Mesh.Color.Fifteen = {84,0,255};
Mesh.Color.Sixteen = {68,255,0};
Mesh.Color.Seventeen = {104,0,255};
Mesh.Color.Eighteen = {0,255,52};
Mesh.Color.Nineteen = {184,0,255};
Solver.Executable0 = "";
Solver.Executable1 = "";
Solver.Executable2 = "";
Solver.Executable3 = "";
Solver.Executable4 = "";
Solver.Executable5 = "";
Solver.Executable6 = "";
Solver.Executable7 = "";
Solver.Executable8 = "";
Solver.Executable9 = "";
Solver.Name0 = "GetDP";
Solver.Name1 = "";
Solver.Name2 = "";
Solver.Name3 = "";
Solver.Name4 = "";
Solver.Name5 = "";
Solver.Name6 = "";
Solver.Name7 = "";
Solver.Name8 = "";
Solver.Name9 = "";
Solver.Extension0 = ".pro";
Solver.Extension1 = "";
Solver.Extension2 = "";
Solver.Extension3 = "";
Solver.Extension4 = "";
Solver.Extension5 = "";
Solver.Extension6 = "";
Solver.Extension7 = "";
Solver.Extension8 = "";
Solver.Extension9 = "";
Solver.OctaveInterpreter = "octave";
Solver.PythonInterpreter = "python";
Solver.RemoteLogin0 = "";
Solver.RemoteLogin1 = "";
Solver.RemoteLogin2 = "";
Solver.RemoteLogin3 = "";
Solver.RemoteLogin4 = "";
Solver.RemoteLogin5 = "";
Solver.RemoteLogin6 = "";
Solver.RemoteLogin7 = "";
Solver.RemoteLogin8 = "";
Solver.RemoteLogin9 = "";
Solver.SocketName = ".gmshsock";
Solver.AlwaysListen = 0;
Solver.AutoArchiveOutputFiles = 0;
Solver.AutoCheck = 1;
Solver.AutoLoadDatabase = 0;
Solver.AutoSaveDatabase = 1;
Solver.AutoMesh = 2;
Solver.AutoMergeFile = 1;
Solver.AutoShowViews = 2;
Solver.AutoShowLastStep = 1;
Solver.Plugins = 0;
Solver.ShowInvisibleParameters = 0;
Solver.Timeout = 5;
PostProcessing.DoubleClickedGraphPointCommand = "";
PostProcessing.GraphPointCommand = "";
PostProcessing.AnimationDelay = 0.1;
PostProcessing.AnimationCycle = 0;
PostProcessing.AnimationStep = 1;
PostProcessing.CombineRemoveOriginal = 1;
PostProcessing.DoubleClickedGraphPointX = 0;
PostProcessing.DoubleClickedGraphPointY = 0;
PostProcessing.DoubleClickedView = 0;
PostProcessing.ForceElementData = 0;
PostProcessing.ForceNodeData = 0;
PostProcessing.Format = 10;
PostProcessing.GraphPointX = 0;
PostProcessing.GraphPointY = 0;
PostProcessing.HorizontalScales = 1;
PostProcessing.Link = 0;
PostProcessing.NbViews = 0;
PostProcessing.Plugins = 1;
PostProcessing.SaveInterpolationMatrices = 1;
PostProcessing.SaveMesh = 1;
PostProcessing.Smoothing = 0;
Print.ParameterCommand = "Mesh.Clip=1; View.Clip=1; General.ClipWholeElements=1; General.Clip0D=Print.Parameter; SetChanged;";
Print.Parameter = 0;
Print.ParameterFirst = -1;
Print.ParameterLast = 1;
Print.ParameterSteps = 10;
Print.Background = 0;
Print.CompositeWindows = 0;
Print.PgfTwoDim = 1;
Print.PgfExportAxis = 0;
Print.PgfHorizontalBar = 0;
Print.DeleteTemporaryFiles = 1;
Print.EpsBestRoot = 1;
Print.EpsCompress = 0;
Print.EpsLineWidthFactor = 1;
Print.EpsOcclusionCulling = 1;
Print.EpsPointSizeFactor = 1;
Print.EpsPS3Shading = 0;
Print.EpsQuality = 1;
Print.Format = 10;
Print.GeoLabels = 1;
Print.GeoOnlyPhysicals = 0;
Print.GifDither = 0;
Print.GifInterlace = 0;
Print.GifSort = 1;
Print.GifTransparent = 0;
Print.Height = -1;
Print.JpegQuality = 100;
Print.JpegSmoothing = 0;
Print.PostElementary = 1;
Print.PostElement = 0;
Print.PostGamma = 0;
Print.PostEta = 0;
Print.PostSICN = 0;
Print.PostSIGE = 0;
Print.PostDisto = 0;
Print.TexAsEquation = 0;
Print.Text = 1;
Print.X3dCompatibility = 0;
Print.X3dPrecision = 1e-09;
Print.X3dRemoveInnerBorders = 0;
Print.X3dTransparency = 0;
Print.Width = -1;

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

View file

@ -19,6 +19,7 @@ MODULE moduleMesh0D
PROCEDURE, PASS:: randPos => randPos0D
PROCEDURE, NOPASS:: fPsi => fPsi0D
PROCEDURE, PASS:: gatherEF => gatherEF0D
PROCEDURE, PASS:: gatherMF => gatherMF0D
PROCEDURE, PASS:: elemK => elemK0D
PROCEDURE, PASS:: elemF => elemF0D
PROCEDURE, PASS:: phy2log => phy2log0D
@ -62,6 +63,7 @@ MODULE moduleMesh0D
!Inits dummy 0D volume
SUBROUTINE initVol0D(self, n, p, nodes)
USE moduleRefParam
USE moduleSpecies
IMPLICIT NONE
CLASS(meshVol0D), INTENT(out):: self
@ -75,10 +77,11 @@ MODULE moduleMesh0D
self%volume = 1.D0
self%n1%v = 1.D0
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initVol0D
PURE FUNCTION getNodes0D(self) RESULT(n)
@ -123,6 +126,17 @@ MODULE moduleMesh0D
END FUNCTION gatherEF0D
PURE FUNCTION gatherMF0D(self, xi) RESULT(MF)
IMPLICIT NONE
CLASS(meshVol0D), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: MF(1:3)
MF = 0.D0
END FUNCTION gatherMF0D
PURE FUNCTION elemK0D(self) RESULT(localK)
IMPLICIT NONE

View file

@ -77,6 +77,7 @@ MODULE moduleMesh1DCart
PROCEDURE, PASS:: elemF => elemFSegm
PROCEDURE, NOPASS:: inside => insideSegm
PROCEDURE, PASS:: gatherEF => gatherEFSegm
PROCEDURE, PASS:: gatherMF => gatherMFSegm
PROCEDURE, PASS:: getNodes => getNodesSegm
PROCEDURE, PASS:: phy2log => phy2logSegm
PROCEDURE, PASS:: nextElement => nextElementSegm
@ -215,10 +216,11 @@ MODULE moduleMesh1DCart
self%n1%v = self%n1%v + self%arNodes(1)
self%n2%v = self%n2%v + self%arNodes(2)
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initVol1DCartSegm
!Calculates a random position in 1D volume
@ -388,6 +390,28 @@ MODULE moduleMesh1DCart
END FUNCTION gatherEFSegm
PURE FUNCTION gatherMFSegm(self, xi) RESULT(MF)
IMPLICIT NONE
CLASS(meshVol1DCartSegm), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: fPsi(1:2)
REAL(8):: MF_Nodes(1:2, 1:3)
REAL(8):: MF(1:3)
REAL(8):: invJ
MF_Nodes(1:2,1) = (/ self%n1%emData%B(1), &
self%n2%emData%B(1) /)
MF_Nodes(1:2,2) = (/ self%n1%emData%B(2), &
self%n2%emData%B(2) /)
MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), &
self%n2%emData%B(3) /)
fPsi = self%fPsi(xi)
MF = MATMUL(fPsi, MF_Nodes)
END FUNCTION gatherMFSegm
!Get nodes from 1D volume
PURE FUNCTION getNodesSegm(self) RESULT(n)
IMPLICIT NONE

View file

@ -78,6 +78,7 @@ MODULE moduleMesh1DRad
PROCEDURE, PASS:: elemF => elemFRad
PROCEDURE, NOPASS:: inside => insideRad
PROCEDURE, PASS:: gatherEF => gatherEFRad
PROCEDURE, PASS:: gatherMF => gatherMFRad
PROCEDURE, PASS:: getNodes => getNodesRad
PROCEDURE, PASS:: phy2log => phy2logRad
PROCEDURE, PASS:: nextElement => nextElementRad
@ -217,10 +218,11 @@ MODULE moduleMesh1DRad
self%n1%v = self%n1%v + self%arNodes(1)
self%n2%v = self%n2%v + self%arNodes(2)
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initVol1DRadSegm
!Calculates a random position in 1D volume
@ -400,6 +402,28 @@ MODULE moduleMesh1DRad
END FUNCTION gatherEFRad
PURE FUNCTION gatherMFRad(self, xi) RESULT(MF)
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: fPsi(1:2)
REAL(8):: MF_Nodes(1:2, 1:3)
REAL(8):: MF(1:3)
REAL(8):: invJ
MF_Nodes(1:2,1) = (/ self%n1%emData%B(1), &
self%n2%emData%B(1) /)
MF_Nodes(1:2,2) = (/ self%n1%emData%B(2), &
self%n2%emData%B(2) /)
MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), &
self%n2%emData%B(3) /)
fPsi = self%fPsi(xi)
MF = MATMUL(fPsi, MF_Nodes)
END FUNCTION gatherMFRad
!Get nodes from 1D volume
PURE FUNCTION getNodesRad(self) RESULT(n)
IMPLICIT NONE

View file

@ -85,6 +85,7 @@ MODULE moduleMesh2DCart
PROCEDURE, PASS:: elemF => elemFQuad
PROCEDURE, NOPASS:: inside => insideQuad
PROCEDURE, PASS:: gatherEF => gatherEFQuad
PROCEDURE, PASS:: gatherMF => gatherMFQuad
PROCEDURE, PASS:: getNodes => getNodesQuad
PROCEDURE, PASS:: phy2log => phy2logQuad
PROCEDURE, PASS:: nextElement => nextElementQuad
@ -114,6 +115,7 @@ MODULE moduleMesh2DCart
PROCEDURE, PASS:: elemF => elemFTria
PROCEDURE, NOPASS:: inside => insideTria
PROCEDURE, PASS:: gatherEF => gatherEFTria
PROCEDURE, PASS:: gatherMF => gatherMFTria
PROCEDURE, PASS:: getNodes => getNodesTria
PROCEDURE, PASS:: phy2log => phy2logTria
PROCEDURE, PASS:: nextElement => nextElementTria
@ -305,10 +307,11 @@ MODULE moduleMesh2DCart
self%n3%v = self%n3%v + self%arNodes(3)
self%n4%v = self%n4%v + self%arNodes(4)
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initVolQuad2DCart
!Computes element area
@ -505,6 +508,33 @@ MODULE moduleMesh2DCart
END FUNCTION gatherEFQuad
PURE FUNCTION gatherMFQuad(self,xi) RESULT(MF)
IMPLICIT NONE
CLASS(meshVol2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: fPsi(1:4)
REAL(8):: MF_Nodes(1:4,1:3)
REAL(8):: MF(1:3)
MF_Nodes(1:4,1) = (/self%n1%emData%B(1), &
self%n2%emData%B(1), &
self%n3%emData%B(1), &
self%n4%emData%B(1) /)
MF_Nodes(1:4,2) = (/self%n1%emData%B(2), &
self%n2%emData%B(2), &
self%n3%emData%B(2), &
self%n4%emData%B(2) /)
MF_Nodes(1:4,3) = (/self%n1%emData%B(3), &
self%n2%emData%B(3), &
self%n3%emData%B(3), &
self%n4%emData%B(3) /)
fPsi = self%fPsi(xi)
MF = MATMUL(fPsi(:), MF_Nodes)
END FUNCTION gatherMFQuad
!Gets nodes from quadrilateral element
PURE FUNCTION getNodesQuad(self) RESULT(n)
IMPLICIT NONE
@ -605,10 +635,11 @@ MODULE moduleMesh2DCart
self%n2%v = self%n2%v + self%arNodes(2)
self%n3%v = self%n3%v + self%arNodes(3)
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initVolTria2DCart
!Random position in quadrilateral volume
@ -816,6 +847,30 @@ MODULE moduleMesh2DCart
END FUNCTION gatherEFTria
PURE FUNCTION gatherMFTria(self,xi) RESULT(MF)
IMPLICIT NONE
CLASS(meshVol2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: fPsi(1:3)
REAL(8):: MF_Nodes(1:3,1:3)
REAL(8):: MF(1:3)
MF_Nodes(1:3,1) = (/self%n1%emData%B(1), &
self%n2%emData%B(1), &
self%n3%emData%B(1) /)
MF_Nodes(1:3,2) = (/self%n1%emData%B(2), &
self%n2%emData%B(2), &
self%n3%emData%B(2) /)
MF_Nodes(1:3,3) = (/self%n1%emData%B(3), &
self%n2%emData%B(3), &
self%n3%emData%B(3) /)
fPsi = self%fPsi(xi)
MF = MATMUL(fPsi, MF_Nodes)
END FUNCTION gatherMFTria
!Gets node indexes from triangular element
PURE FUNCTION getNodesTria(self) RESULT(n)
IMPLICIT NONE

View file

@ -86,6 +86,7 @@ MODULE moduleMesh2DCyl
PROCEDURE, PASS:: elemF => elemFQuad
PROCEDURE, NOPASS:: inside => insideQuad
PROCEDURE, PASS:: gatherEF => gatherEFQuad
PROCEDURE, PASS:: gatherMF => gatherMFQuad
PROCEDURE, PASS:: getNodes => getNodesQuad
PROCEDURE, PASS:: phy2log => phy2logQuad
PROCEDURE, PASS:: nextElement => nextElementQuad
@ -115,6 +116,7 @@ MODULE moduleMesh2DCyl
PROCEDURE, PASS:: elemF => elemFTria
PROCEDURE, NOPASS:: inside => insideTria
PROCEDURE, PASS:: gatherEF => gatherEFTria
PROCEDURE, PASS:: gatherMF => gatherMFTria
PROCEDURE, PASS:: getNodes => getNodesTria
PROCEDURE, PASS:: phy2log => phy2logTria
PROCEDURE, PASS:: nextElement => nextElementTria
@ -293,10 +295,11 @@ MODULE moduleMesh2DCyl
self%n3%v = self%n3%v + self%arNodes(3)
self%n4%v = self%n4%v + self%arNodes(4)
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initVolQuad2DCyl
!Computes element area
@ -526,6 +529,33 @@ MODULE moduleMesh2DCyl
END FUNCTION gatherEFQuad
PURE FUNCTION gatherMFQuad(self,xi) RESULT(MF)
IMPLICIT NONE
CLASS(meshVol2DCylQuad), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: fPsi(1:4)
REAL(8):: MF_Nodes(1:4,1:3)
REAL(8):: MF(1:3)
MF_Nodes(1:4,1) = (/self%n1%emData%B(1), &
self%n2%emData%B(1), &
self%n3%emData%B(1), &
self%n4%emData%B(1) /)
MF_Nodes(1:4,2) = (/self%n1%emData%B(2), &
self%n2%emData%B(2), &
self%n3%emData%B(2), &
self%n4%emData%B(2) /)
MF_Nodes(1:4,3) = (/self%n1%emData%B(3), &
self%n2%emData%B(3), &
self%n3%emData%B(3), &
self%n4%emData%B(3) /)
fPsi = self%fPsi(xi)
MF = MATMUL(fPsi(:), MF_Nodes)
END FUNCTION gatherMFQuad
!Gets nodes from quadrilateral element
PURE FUNCTION getNodesQuad(self) RESULT(n)
IMPLICIT NONE
@ -626,10 +656,11 @@ MODULE moduleMesh2DCyl
self%n2%v = self%n2%v + self%arNodes(2)
self%n3%v = self%n3%v + self%arNodes(3)
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initVolTria2DCyl
!Random position in quadrilateral volume
@ -845,6 +876,30 @@ MODULE moduleMesh2DCyl
END FUNCTION gatherEFTria
PURE FUNCTION gatherMFTria(self,xi) RESULT(MF)
IMPLICIT NONE
CLASS(meshVol2DCylTria), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: fPsi(1:3)
REAL(8):: MF_Nodes(1:3,1:3)
REAL(8):: MF(1:3)
MF_Nodes(1:3,1) = (/self%n1%emData%B(1), &
self%n2%emData%B(1), &
self%n3%emData%B(1) /)
MF_Nodes(1:3,2) = (/self%n1%emData%B(2), &
self%n2%emData%B(2), &
self%n3%emData%B(2) /)
MF_Nodes(1:3,3) = (/self%n1%emData%B(3), &
self%n2%emData%B(3), &
self%n3%emData%B(3) /)
fPsi = self%fPsi(xi)
MF = MATMUL(fPsi, MF_Nodes)
END FUNCTION gatherMFTria
!Gets node indexes from triangular element
PURE FUNCTION getNodesTria(self) RESULT(n)
IMPLICIT NONE

View file

@ -78,6 +78,7 @@ MODULE moduleMesh3DCart
PROCEDURE, PASS:: elemF => elemFTetra
PROCEDURE, NOPASS:: inside => insideTetra
PROCEDURE, PASS:: gatherEF => gatherEFTetra
PROCEDURE, PASS:: gatherMF => gatherMFTetra
PROCEDURE, PASS:: getNodes => getNodesTetra
PROCEDURE, PASS:: phy2log => phy2logTetra
PROCEDURE, PASS:: nextElement => nextElementTetra
@ -280,10 +281,11 @@ MODULE moduleMesh3DCart
self%n3%v = self%n3%v + volNodes(3)
self%n4%v = self%n4%v + volNodes(4)
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initVolTetra3DCart
!Random position in volume tetrahedron
@ -498,6 +500,33 @@ MODULE moduleMesh3DCart
END FUNCTION gatherEFTetra
PURE FUNCTION gatherMFTetra(self, xi) RESULT(MF)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: fPsi(1:4)
REAL(8):: MF_Nodes(1:4,1:3)
REAL(8):: MF(1:3)
MF_Nodes(1:4,1) = (/self%n1%emData%B(1), &
self%n2%emData%B(1), &
self%n3%emData%B(1), &
self%n4%emData%B(1) /)
MF_Nodes(1:4,2) = (/self%n1%emData%B(2), &
self%n2%emData%B(2), &
self%n3%emData%B(2), &
self%n4%emData%B(2) /)
MF_Nodes(1:4,3) = (/self%n1%emData%B(3), &
self%n2%emData%B(3), &
self%n3%emData%B(3), &
self%n4%emData%B(3) /)
fPsi = self%fPsi(xi)
MF = MATMUL(fPsi, MF_Nodes)
END FUNCTION gatherMFTetra
PURE FUNCTION getNodesTetra(self) RESULT(n)
IMPLICIT NONE

View file

@ -45,6 +45,7 @@ MODULE moduleMeshOutput0D
CLASS(meshGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: t
CHARACTER(:), ALLOCATABLE:: fileName
INTEGER:: k
fileName='OUTPUT_Collisions.dat'
IF (t == 0) THEN
@ -56,7 +57,7 @@ MODULE moduleMeshOutput0D
END IF
OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
WRITE(20, "(ES20.6E3, I20)") REAL(t)*tauMin*ti_ref, self%vols(1)%obj%nColl
WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%vols(1)%obj%tallyColl(k)%tally, k=1,nCollPairs)
CLOSE(20)
END SUBROUTINE printColl0D

View file

@ -67,24 +67,35 @@ MODULE moduleMeshInputGmsh2
!Read the nodes information
DO e = 1, self%numNodes
READ(10, *) n, r(1), r(2), r(3)
SELECT CASE(self%geometry)
CASE("3DCart")
SELECT CASE(self%dimen)
CASE(3)
ALLOCATE(meshNode3Dcart::self%nodes(n)%obj)
self%connectMesh => connectMesh3DCart
CASE("2DCyl")
ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj)
CASE(2)
SELECT CASE(self%geometry)
CASE("Cyl")
ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj)
self%connectMesh => connectMesh2DCyl
CASE("Cart")
ALLOCATE(meshNode2DCart:: self%nodes(n)%obj)
self%connectMesh => connectMesh2DCart
END SELECT
r(3) = 0.D0
CASE("2DCart")
ALLOCATE(meshNode2DCart:: self%nodes(n)%obj)
r(3) = 0.D0
CASE(1)
SELECT CASE(self%geometry)
CASE("Rad")
ALLOCATE(meshNode1DRad:: self%nodes(n)%obj)
self%connectMesh => connectMesh1DRad
CASE("1DRad")
ALLOCATE(meshNode1DRad:: self%nodes(n)%obj)
r(2:3) = 0.D0
CASE("Cart")
ALLOCATE(meshNode1DCart:: self%nodes(n)%obj)
self%connectMesh => connectMesh1DCart
CASE("1DCart")
ALLOCATE(meshNode1DCart:: self%nodes(n)%obj)
END SELECT
r(2:3) = 0.D0
END SELECT
@ -92,7 +103,6 @@ MODULE moduleMeshInputGmsh2
END DO
!Skip comments
READ(10, *)
READ(10, *)
@ -106,16 +116,16 @@ MODULE moduleMeshInputGmsh2
self%numEdges = 0
DO e = 1, totalNumElem
READ(10, *) eTemp, elemType
SELECT CASE(self%geometry)
CASE("3DCart")
SELECT CASE(self%dimen)
CASE(3)
!Element type 2 is triangle in gmsh
IF (elemType == 2) self%numEdges = e
CASE("2DCyl","2DCart")
CASE(2)
!Element type 1 is segment in Gmsh
IF (elemType == 1) self%numEdges = e
CASE("1DRad","1DCart")
CASE(1)
!Element type 15 is physical point in Gmsh
IF (elemType == 15) self%numEdges = e
@ -148,8 +158,8 @@ MODULE moduleMeshInputGmsh2
!Reads edges
DO e=1, self%numEdges
!Reads the edge according to the geometry
SELECT CASE(self%geometry)
CASE("3DCart")
SELECT CASE(self%dimen)
CASE(3)
READ(10, *) n, elemType, eTemp, boundaryType
BACKSPACE(10)
@ -167,41 +177,49 @@ MODULE moduleMeshInputGmsh2
END SELECT
CASE("2DCyl")
ALLOCATE(p(1:2))
CASE (2)
SELECT CASE(self%geometry)
CASE("Cyl")
ALLOCATE(p(1:2))
READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
!Associate boundary condition procedure.
bt = getBoundaryId(boundaryType)
READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
!Associate boundary condition procedure.
bt = getBoundaryId(boundaryType)
ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj)
ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj)
CASE("2DCart")
ALLOCATE(p(1:2))
CASE("Cart")
ALLOCATE(p(1:2))
READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
!Associate boundary condition procedure.
bt = getBoundaryId(boundaryType)
READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
!Associate boundary condition procedure.
bt = getBoundaryId(boundaryType)
ALLOCATE(meshEdge2DCart:: self%edges(e)%obj)
ALLOCATE(meshEdge2DCart:: self%edges(e)%obj)
CASE("1DRad")
ALLOCATE(p(1:1))
END SELECT
READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
!Associate boundary condition
bt = getBoundaryId(boundaryType)
CASE(1)
SELECT CASE(self%geometry)
CASE("Rad")
ALLOCATE(p(1:1))
ALLOCATE(meshEdge1DRad:: self%edges(e)%obj)
READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
!Associate boundary condition
bt = getBoundaryId(boundaryType)
CASE("1DCart")
ALLOCATE(p(1:1))
ALLOCATE(meshEdge1DRad:: self%edges(e)%obj)
READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
!Associate boundary condition
bt = getBoundaryId(boundaryType)
CASE("Cart")
ALLOCATE(p(1:1))
ALLOCATE(meshEdge1DCart:: self%edges(e)%obj)
READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
!Associate boundary condition
bt = getBoundaryId(boundaryType)
ALLOCATE(meshEdge1DCart:: self%edges(e)%obj)
END SELECT
END SELECT
@ -215,8 +233,8 @@ MODULE moduleMeshInputGmsh2
!Read and initialize volumes
DO e = 1, self%numVols
!Reads the volume according to the geometry
SELECT CASE(self%geometry)
CASE("3DCart")
SELECT CASE(self%dimen)
CASE(3)
READ(10, *) n, elemType
BACKSPACE(10)
@ -229,56 +247,64 @@ MODULE moduleMeshInputGmsh2
END SELECT
CASE("2DCyl")
READ(10,*) n, elemType
BACKSPACE(10)
CASE(2)
SELECT CASE(self%geometry)
CASE("Cyl")
READ(10,*) n, elemType
BACKSPACE(10)
SELECT CASE(elemType)
CASE (2)
!Triangular element
ALLOCATE(p(1:3))
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
ALLOCATE(meshVol2DCylTria:: self%vols(e)%obj)
SELECT CASE(elemType)
CASE (2)
!Triangular element
ALLOCATE(p(1:3))
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
ALLOCATE(meshVol2DCylTria:: self%vols(e)%obj)
CASE (3)
!Quadrilateral element
ALLOCATE(p(1:4))
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
ALLOCATE(meshVol2DCylQuad:: self%vols(e)%obj)
CASE (3)
!Quadrilateral element
ALLOCATE(p(1:4))
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
ALLOCATE(meshVol2DCylQuad:: self%vols(e)%obj)
END SELECT
CASE("Cart")
READ(10,*) n, elemType
BACKSPACE(10)
SELECT CASE(elemType)
CASE (2)
!Triangular element
ALLOCATE(p(1:3))
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
ALLOCATE(meshVol2DCartTria:: self%vols(e)%obj)
CASE (3)
!Quadrilateral element
ALLOCATE(p(1:4))
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
ALLOCATE(meshVol2DCartQuad:: self%vols(e)%obj)
END SELECT
END SELECT
CASE("2DCart")
READ(10,*) n, elemType
BACKSPACE(10)
CASE(1)
SELECT CASE(self%geometry)
CASE("Rad")
ALLOCATE(p(1:2))
SELECT CASE(elemType)
CASE (2)
!Triangular element
ALLOCATE(p(1:3))
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
ALLOCATE(meshVol2DCartTria:: self%vols(e)%obj)
READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2)
ALLOCATE(meshVol1DRadSegm:: self%vols(e)%obj)
CASE (3)
!Quadrilateral element
ALLOCATE(p(1:4))
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
ALLOCATE(meshVol2DCartQuad:: self%vols(e)%obj)
CASE("Cart")
ALLOCATE(p(1:2))
READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2)
ALLOCATE(meshVol1DCartSegm:: self%vols(e)%obj)
END SELECT
CASE("1DRad")
ALLOCATE(p(1:2))
READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2)
ALLOCATE(meshVol1DRadSegm:: self%vols(e)%obj)
CASE("1DCart")
ALLOCATE(p(1:2))
READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2)
ALLOCATE(meshVol1DCartSegm:: self%vols(e)%obj)
END SELECT
CALL self%vols(e)%obj%init(n - numEdges, p, self%nodes)
@ -288,6 +314,9 @@ MODULE moduleMeshInputGmsh2
CLOSE(10)
!Call mesh connectivity
CALL self%connectMesh
END SUBROUTINE readGmsh2
!Reads the initial information from an output file for an species

View file

@ -98,6 +98,7 @@ MODULE moduleMeshOutputGmsh2
CLASS(meshGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: t
INTEGER:: numEdges
INTEGER:: k, c
INTEGER:: n
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
@ -125,19 +126,23 @@ MODULE moduleMeshOutputGmsh2
WRITE(60, "(A)") '$MeshFormat'
WRITE(60, "(A)") '2.2 0 8'
WRITE(60, "(A)") '$EndMeshFormat'
WRITE(60, "(A)") '$ElementData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Collisions"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numVols
DO n=1, self%numVols
WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%nColl
DO k = 1, nCollPairs
DO c = 1, interactionMatrix(k)%amount
WRITE(60, "(A)") '$ElementData'
WRITE(60, "(A)") '1'
WRITE(60, "(5A,I2)") '"Pair ', interactionMatrix(k)%sp_i%name, '-', interactionMatrix(k)%sp_j%name, ' collision ', c
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numVols
DO n=1, self%numVols
WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%tallyColl(k)%tally(c)
END DO
WRITE(60, "(A)") '$EndElementData'
END DO
END DO
WRITE(60, "(A)") '$EndElementData'
CLOSE(60)
@ -200,6 +205,21 @@ MODULE moduleMeshOutputGmsh2
WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref
END DO
WRITE(20, "(A)") '$EndElementData'
WRITE(20, "(A)") '$NodeData'
WRITE(20, "(A)") '1'
WRITE(20, "(A)") '"Magnetic Field (T)"'
WRITE(20, *) 1
WRITE(20, *) time
WRITE(20, *) 3
WRITE(20, *) t
WRITE(20, *) 3
WRITE(20, *) self%numNodes
DO n=1, self%numNodes
WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref
END DO
WRITE(20, "(A)") '$EndNodeData'
CLOSE(20)
END IF

View file

@ -3,6 +3,7 @@ MODULE moduleMesh
USE moduleList
USE moduleOutput
USE moduleBoundary
USE moduleCollisions
IMPLICIT NONE
!Generic mesh element
@ -148,17 +149,17 @@ MODULE moduleMesh
!Parent of Volume element
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshVol
!Maximum collision rate
REAL(8):: sigmaVrelMax = 0.D0
REAL(8), ALLOCATABLE:: sigmaVrelMax(:)
!Arrays for counting number of collisions
TYPE(tallyCollisions), ALLOCATABLE:: tallyColl(:)
!Volume
REAL(8):: volume = 0.D0
!List of particles inside the volume
TYPE(listNode):: listPart_in
TYPE(listNode), ALLOCATABLE:: listPart_in(:)
!Lock indicator for listPart_in
INTEGER(KIND=OMP_LOCK_KIND):: lock
!Number of collisions per volume
INTEGER:: nColl = 0
!Total weight of particles inside cell
REAL(8):: totalWeight = 0.D0
REAL(8), ALLOCATABLE:: totalWeight(:)
CONTAINS
PROCEDURE(initVol_interface), DEFERRED, PASS:: init
PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes
@ -166,6 +167,7 @@ MODULE moduleMesh
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
PROCEDURE, PASS:: scatter
PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF
PROCEDURE(gatherMF_interface), DEFERRED, PASS:: gatherMF
PROCEDURE(elemK_interface), DEFERRED, PASS:: elemK
PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF
PROCEDURE, PASS:: findCell
@ -194,6 +196,14 @@ MODULE moduleMesh
END FUNCTION gatherEF_interface
PURE FUNCTION gatherMF_interface(self, xi) RESULT(MF)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: MF(1:3)
END FUNCTION gatherMF_interface
PURE FUNCTION getNodesVol_interface(self) RESULT(n)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
@ -262,6 +272,8 @@ MODULE moduleMesh
!Generic mesh type
TYPE, ABSTRACT:: meshGeneric
!Dimension of the mesh
INTEGER:: dimen
!Geometry of the mesh
CHARACTER(:), ALLOCATABLE:: geometry
!Number of elements
@ -505,6 +517,7 @@ MODULE moduleMesh
CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
REAL(8):: xi(1:3)
CLASS(meshElement), POINTER:: nextElement
INTEGER:: sp
xi = self%phy2log(part%r)
!Checks if particle is inside 'self' cell
@ -514,8 +527,9 @@ MODULE moduleMesh
part%n_in = .TRUE.
!Assign particle to listPart_in
CALL OMP_SET_LOCK(self%lock)
CALL self%listPart_in%add(part)
self%totalWeight = self%totalWeight + part%weight
sp = part%species%n
CALL self%listPart_in(sp)%add(part)
self%totalWeight(sp) = self%totalWeight(sp) + part%weight
CALL OMP_UNSET_LOCK(self%lock)
ELSE
@ -575,6 +589,7 @@ MODULE moduleMesh
CLASS(meshVol), POINTER:: vol
REAL(8), DIMENSION(1:3):: xii
CLASS(meshElement), POINTER:: nextElement
INTEGER:: sp
found = .FALSE.
@ -584,8 +599,9 @@ MODULE moduleMesh
IF (vol%inside(xii)) THEN
part%volColl = vol%n
CALL OMP_SET_LOCK(vol%lock)
CALL vol%listPart_in%add(part)
vol%totalWeight = vol%totalWeight + part%weight
sp = part%species%n
CALL vol%listPart_in(sp)%add(part)
vol%totalWeight(sp) = vol%totalWeight(sp) + part%weight
CALL OMP_UNSET_LOCK(vol%lock)
found = .TRUE.
@ -643,68 +659,154 @@ MODULE moduleMesh
USE moduleList
use moduleRefParam
USE moduleRandom
USE moduleOutput
USE moduleMath
IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout), TARGET:: self
INTEGER, INTENT(in):: t
INTEGER:: e
CLASS(meshVol), POINTER:: vol
INTEGER:: nPart !Number of particles inside the cell
INTEGER:: k, i, j
INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell
REAL(8):: pMax !Maximum probability of collision
INTEGER:: rnd !random index
INTEGER:: nColl
TYPE(pointerArray), ALLOCATABLE:: partTemp_i(:), partTemp_j(:)
TYPE(particle), POINTER:: part_i, part_j
INTEGER:: n !collision
INTEGER:: ij, k
REAL(8):: sigmaVrelMaxNew
TYPE(pointerArray), ALLOCATABLE:: partTemp(:)
INTEGER:: n, c
REAL(8):: vRel, rMass, eRel
REAL(8):: sigmaVrelTotal
REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:)
REAL(8):: rnd !Random number for collision
IF (MOD(t, everyColl) == 0) THEN
!Collisions need to be performed in this iteration
!$OMP DO SCHEDULE(DYNAMIC)
!$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j)
DO e=1, self%numVols
vol => self%vols(e)%obj
nPart = vol%listPart_in%amount
!Resets the number of collisions
vol%nColl = 0
!Calculates number of collisions if there is more than one particle in the cell
IF (nPart > 1) THEN
!Probability of collision
pMax = vol%totalWeight*vol%sigmaVrelMax*tauColl/vol%volume
!Number of collisions in the cell
vol%nColl = NINT(REAL(nPart)*pMax*0.5D0)
IF (vol%nColl > 0) THEN
!Converts the list of particles to an array for easy access
partTemp = vol%listPart_in%convert2Array()
!TODO: Simplify this, to many sublevels
!Iterate over the number of pairs
DO k = 1, nCollPairs
!Reset tally of collisions
IF (collOutput) THEN
vol%tallyColl(k)%tally = 0
END IF
DO n = 1, vol%nColl
!Select random numbers
rnd = random(1, nPart)
part_i => partTemp(rnd)%part
rnd = random(1, nPart)
part_j => partTemp(rnd)%part
ij = interactionIndex(part_i%species%n, part_j%species%n)
sigmaVrelMaxNew = 0.D0
DO k = 1, interactionMatrix(ij)%amount
CALL interactionMatrix(ij)%collisions(k)%obj%collide(vol%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j)
IF (interactionMatrix(k)%amount > 0) THEN
!Select the species for the collision pair
i = interactionMatrix(k)%sp_i%n
j = interactionMatrix(k)%sp_j%n
END DO
!Number of particles per species in the collision pair
nPart_i = vol%listPart_in(i)%amount
nPart_j = vol%listPart_in(j)%amount
!Update maximum cross section*v_rel per each collision
IF (sigmaVrelMaxNew > vol%sigmaVrelMax) THEN
vol%sigmaVrelMax = sigmaVrelMaxNew
IF (nPart_i > 0 .AND. nPart_j > 0) THEN
!Total number of particles for the collision pair
nPart = nPart_i + nPart_j
!Resets the number of collisions in the cell
nColl = 0
!Probability of collision for pair i-j
pMax = (vol%totalWeight(i) + vol%totalWeight(j))*vol%sigmaVrelMax(k)*tauColl/vol%volume
!Number of collisions in the cell
nColl = NINT(REAL(nPart)*pMax*0.5D0)
!Converts the list of particles to an array for easy access
IF (nColl > 0) THEN
partTemp_i = vol%listPart_in(i)%convert2Array()
partTemp_j = vol%listPart_in(j)%convert2Array()
END IF
DO n = 1, nColl
!Select random particles
part_i => NULL()
part_j => NULL()
rnd = random(1, nPart_i)
part_i => partTemp_i(rnd)%part
rnd = random(1, nPart_j)
part_j => partTemp_j(rnd)%part
!If they are the same particle, skip
!TODO: Maybe try to improve this
IF (ASSOCIATED(part_i, part_j)) THEN
CYCLE
END IF
!If particles do not belong to the species, skip collision
!This can happen, for example, if particle has been previously ionized or removed
!TODO: Try to find a way to no lose these collisions. Maybe check new 'k' and use that for the collision, maybe?
IF (part_i%species%n /= i .OR. &
part_j%species%n /= j) THEN
CYCLE
END IF
!Obtain the cross sections for the different processes
!TODO: From here it might be a procedure in interactionMatrix
vRel = NORM2(part_i%v-part_j%v)
rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
eRel = rMass*vRel**2
CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)
!Update maximum sigma*v_rel
IF (sigmaVrelTotal > vol%sigmaVrelMax(k)) THEN
vol%sigmaVrelMax(k) = sigmaVrelTotal
END IF
ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount))
probabilityColl = 0.0
DO c = 1, interactionMatrix(k)%amount
probabilityColl(c) = sigmaVrel(c)/vol%sigmaVrelMax(k) + SUM(probabilityColl(0:c-1))
END DO
!Selects random number between 0 and 1
rnd = random()
!If the random number is below the total probability of collision, collide particles
IF (rnd < sigmaVrelTotal / vol%sigmaVrelMax(k)) THEN
!Loop over collisions
DO c = 1, interactionMatrix(k)%amount
IF (rnd <= probabilityColl(c)) THEN
CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel)
!If collisions are gonna be output, count the collision
IF (collOutput) THEN
vol%tallyColl(k)%tally(c) = vol%tallyColl(k)%tally(c) + 1
END IF
!A collision has ocurred, exit the loop
EXIT
END IF
END DO
END IF
!Deallocate arrays for next collision
DEALLOCATE(sigmaVrel, probabilityColl)
!End loop collisions in cell
END DO
END IF
END DO
END IF
END IF
!End loop collision pairs
END DO
!End loop volumes
END DO
!$OMP END DO

View file

@ -110,31 +110,34 @@ MODULE moduleMeshBoundary
USE moduleMesh
USE moduleRefParam
USE moduleRandom
USE moduleMath
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass
REAL(8):: ionizationRate
INTEGER:: ionizationPair, p
INTEGER:: nIonizations !Number of ionizations based on eRel
REAL(8):: pIonization !Probability of ionization of each event
INTEGER:: p
REAL(8):: v0(1:3) !random velocity of neutral
TYPE(particle), POINTER:: newElectron
TYPE(particle), POINTER:: newIon
SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
TYPE IS(boundaryIonization)
mRel = (bound%m0*part%species%m)*(bound%m0+part%species%m)
mRel = reducedMass(bound%m0, part%species%m)
vRel = SUM(DABS(part%v-bound%v0))
eRel = mRel*vRel**2*5.D-1
IF (eRel > bound%eThreshold) THEN
ionizationRate = part%weight*bound%n0*bound%crossSection%get(eRel)*vRel
!Maximum number of possible ionizations based on relative energy
nIonizations = FLOOR(eRel/bound%eThreshold)
!Rounds the number of particles up
ionizationPair = NINT(ionizationRate*bound%effectiveTime/bound%species%weight)
DO p = 1, nIonizations
!Get probability of ionization
pIonization = 1.D0 - DEXP(-bound%n0*bound%crossSection%get(eRel)*vRel*bound%effectiveTime/REAL(nIonizations))
!Create the new pair of particles
DO p = 1, ionizationPair
!If a random number is below the probability of ionization, create new pair of ion-electron
IF (random() < pIonization) THEN
!Assign random velocity to the neutral
v0(1) = bound%v0(1) + bound%vTh*randomMaxwellian()
v0(2) = bound%v0(2) + bound%vTh*randomMaxwellian()
@ -159,14 +162,7 @@ MODULE moduleMeshBoundary
newElectron%xi = mesh%vols(part%vol)%obj%phy2log(newElectron%r)
newIon%xi = newElectron%xi
newElectron%qm = part%qm
SELECT TYPE(spe => bound%species)
TYPE IS(speciesCharged)
newIon%qm = spe%qm
END SELECT
newElectron%weight = bound%species%weight
newElectron%weight = part%weight
newIon%weight = newElectron%weight
newElectron%n_in = .TRUE.
@ -178,9 +174,13 @@ MODULE moduleMeshBoundary
CALL partSurfaces%add(newIon)
CALL OMP_UNSET_LOCK(lockSurfaces)
END DO
!Electron loses energy due to ionization
eRel = eRel - bound%eThreshold
vRel = 2.D0*DSQRT(eRel)/mRel
END IF
END IF
END DO
END SELECT

View file

@ -7,8 +7,6 @@ MODULE moduleCollisions
!Abstract type for collision between two particles
TYPE, ABSTRACT:: collisionBinary
REAL(8):: rMass !Reduced mass
REAL(8):: sMassInv !Summed mass
TYPE(table1D):: crossSec !cross section of collision
CONTAINS
PROCEDURE(collideBinary_interface), PASS, DEFERRED:: collide
@ -16,14 +14,13 @@ MODULE moduleCollisions
END TYPE collisionBinary
ABSTRACT INTERFACE
SUBROUTINE collideBinary_interface(self, sigmaVRelMax, sigmaVrelMaxNew, part_i, part_j)
SUBROUTINE collideBinary_interface(self, part_i, part_j, vRel)
USE moduleSpecies
IMPORT:: collisionBinary
CLASS(collisionBinary), INTENT(in):: self
REAL(8), INTENT(in):: sigmaVrelMax
REAL(8), INTENT(inout):: sigmaVrelMaxNew
TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
REAL(8), INTENT(in):: vRel
END SUBROUTINE
@ -70,13 +67,24 @@ MODULE moduleCollisions
!Type for interaction matrix
TYPE:: interactionsBinary
CLASS(speciesGeneric), POINTER:: sp_i
CLASS(speciesGeneric), POINTER:: sp_j
INTEGER:: amount
TYPE(collisionCont), ALLOCATABLE:: collisions(:)
CONTAINS
PROCEDURE, PASS:: init => initInteractionBinary
PROCEDURE, PASS:: getSigmaVrel => getSigmaVrelBinary
END TYPE interactionsBinary
!Type to count number of collisions
TYPE:: tallyCollisions
INTEGER, ALLOCATABLE:: tally(:)
END TYPE
!Number of collision pairs (nSpecies*(nSpecies+1)/2)
INTEGER:: nCollPairs = 0
!Collision 'Matrix'. A symmetric 2D matrix put into a 1D array to save memory
TYPE(interactionsBinary), ALLOCATABLE, TARGET:: interactionMatrix(:)
!Folder for collision cross section tables
@ -85,13 +93,14 @@ MODULE moduleCollisions
CONTAINS
!Velocity of center of mass of two particles
PURE FUNCTION velocityCM(m_i, v_i, m_j, v_j) RESULT(vCM)
IMPLICIT NONE
REAL(8), INTENT(in):: m_i, m_j
REAL(8), INTENT(in), DIMENSION(1:3):: v_i, v_j
REAL(8):: vCM(1:3)
vCM = (m_i*v_i + m_j*v_j)/(m_i + m_j)
vCM = (m_i*v_i + m_j*v_j) / (m_i + m_j)
END FUNCTION velocityCM
@ -118,10 +127,11 @@ MODULE moduleCollisions
IMPLICIT NONE
TYPE(interactionsBinary), INTENT(inout), ALLOCATABLE:: interactionMatrix(:)
INTEGER:: nInteractions
nInteractions = (nSpecies*(nSpecies+1))/2
ALLOCATE(interactionMatrix(1:nInteractions))
nCollPairs = (nSpecies*(nSpecies+1))/2
ALLOCATE(interactionMatrix(1:nCollPairs))
interactionMatrix(:)%amount = 0
END SUBROUTINE initInteractionMatrix
@ -139,21 +149,52 @@ MODULE moduleCollisions
END FUNCTION interactionIndex
!Inits the binary interaction
SUBROUTINE initInteractionBinary(self, amount)
SUBROUTINE initInteractionBinary(self, amount, i, j)
USE moduleMath
USE moduleSpecies
IMPLICIT NONE
CLASS(interactionsBinary), INTENT(inout):: self
INTEGER, INTENT(in):: amount
INTEGER, INTENT(in):: i, j
REAL(8):: mass_i, mass_j
self%sp_i => species(i)%obj
self%sp_j => species(j)%obj
self%amount = amount
mass_i = species(i)%obj%m
mass_j = species(j)%obj%m
ALLOCATE(self%collisions(1:self%amount))
END SUBROUTINE initInteractionBinary
SUBROUTINE getSigmaVrelBinary (self, vRel, eRel, sigmaVrelTotal, sigmaVrel)
IMPLICIT NONE
CLASS(interactionsBinary), INTENT(in):: self
REAL(8), INTENT(in):: vRel, eRel
REAL(8), INTENT(out):: sigmaVrelTotal
REAL(8), INTENT(out), ALLOCATABLE:: sigmaVrel(:)
INTEGER:: c
sigmaVrelTotal = 0.D0
ALLOCATE(sigmaVrel(1:self%amount))
DO c = 1, self%amount
sigmaVrel(c) = self%collisions(c)%obj%crossSec%get(eRel)*vRel
END DO
sigmaVrelTotal = SUM(sigmaVrel)
END SUBROUTINE getSigmaVrelBinary
!ELASTIC COLLISIONS
!Inits binary elastic collision
SUBROUTINE initBinaryElastic(collision, crossSectionFilename, mass_i, mass_j)
SUBROUTINE initBinaryElastic(collision, crossSectionFilename)
USE moduleTable
USE moduleRefParam
USE moduleConstParam
@ -161,7 +202,6 @@ MODULE moduleCollisions
CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
REAL(8), INTENT(in):: mass_i, mass_j
ALLOCATE(collisionBinaryElastic:: collision)
@ -171,15 +211,10 @@ MODULE moduleCollisions
!Convert to no-dimensional units
CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
!Calculates reduced mass
collision%sMassInv = 1.D0/(mass_i+mass_j)
collision%rMass = (mass_i*mass_j)*collision%sMassInv
END SUBROUTINE initBinaryElastic
!Binary elastic process
SUBROUTINE collideBinaryElastic(self, sigmaVrelMax, sigmaVrelMaxNew, &
part_i, part_j)
SUBROUTINE collideBinaryElastic(self, part_i, part_j, vRel)
USE moduleSpecies
USE moduleConstParam
USE moduleRandom
@ -187,38 +222,26 @@ MODULE moduleCollisions
IMPLICIT NONE
CLASS(collisionBinaryElastic), INTENT(in):: self
REAL(8), INTENT(in):: sigmaVrelMax
REAL(8), INTENT(inout):: sigmaVrelMaxNew
TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
REAL(8):: sigmaVrel
REAL(8):: vRel !relative velocity
REAL(8):: eRel !relative energy
REAL(8), INTENT(in):: vRel
REAL(8):: m_i, m_j
REAL(8), DIMENSION(1:3):: vCM
REAL(8):: vp(1:3)
REAL(8), DIMENSION(1:3):: vCM, vp
vRel = NORM2(part_i%v-part_j%v)
eRel = self%rMass*vRel**2
sigmaVrel = self%crossSec%get(eRel)*vRel
sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
m_i = part_i%species%m
m_j = part_j%species%m
!Applies the collision
vCM = velocityCM(m_i, part_i%v, m_j, part_j%v)
vp = vRel*randomDirectionVHS()
m_i = part_i%species%m
m_j = part_j%species%m
!Applies the collision
vCM = velocityCM(part_i%weight*m_i, part_i%v, part_j%weight*m_j, part_j%v)
vp = vRel*randomDirectionVHS()
!Assign velocities to particles
part_i%v = vCM + m_j*vp*self%sMassInv
part_j%v = vCM - m_i*vp*self%sMassInv
END IF
!Assign velocities to particles
part_i%v = vCM + m_j*vp / (m_i + m_j)
part_j%v = vCM - m_i*vp / (m_i + m_j)
END SUBROUTINE collideBinaryElastic
!ELECTRON IMPACT IONIZATION
!Inits electron impact ionization
SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, mass_i, mass_j, electron)
SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, electron)
USE moduleTable
USE moduleRefParam
USE moduleConstParam
@ -229,7 +252,6 @@ MODULE moduleCollisions
CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
REAL(8), INTENT(in):: energyThreshold
REAL(8), INTENT(in):: mass_i, mass_j
CHARACTER(:), ALLOCATABLE, INTENT(in):: electron
INTEGER:: electronIndex
@ -241,10 +263,6 @@ MODULE moduleCollisions
!Convert to no-dimensional units
CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
!Calculates reduced mass
collision%sMassInv = 1.D0/(mass_i+mass_j)
collision%rMass = (mass_i*mass_j)*collision%sMassInv
!Specific parameters for ionization collision
SELECT TYPE(collision)
TYPE IS(collisionBinaryIonization)
@ -252,23 +270,27 @@ MODULE moduleCollisions
!Input energy is in eV. Convert to J with ev2J and then to
!non-dimensional units.
collision%eThreshold = energyThreshold*eV2J/(m_ref*v_ref**2)
!species for impacting electron
electronIndex = speciesName2Index(electron)
SELECT TYPE(sp => species(electronIndex)%obj)
TYPE IS(speciesCharged)
collision%electron => sp
CLASS DEFAULT
CALL criticalError("Species " // sp%name // " chosen for ionization is not a charged species", 'initBinaryIonization')
CALL criticalError("Species " // sp%name // " chosen for " // &
"secondary electron is not a charged species", 'initBinaryIonization')
END SELECT
!momentum change per ionization process
collision%deltaV = sqrt(collision%eThreshold / collision%electron%m)
END SELECT
END SUBROUTINE initBinaryIonization
!Binary electron impact ionization process
SUBROUTINE collideBinaryIonization(self, sigmaVrelMax, sigmaVrelMaxNew, &
part_i, part_j)
SUBROUTINE collideBinaryIonization(self, part_i, part_j, vRel)
USE moduleSpecies
USE moduleErrors
USE moduleList
@ -278,80 +300,71 @@ MODULE moduleCollisions
IMPLICIT NONE
CLASS(collisionBinaryIonization), INTENT(in):: self
REAL(8), INTENT(in):: sigmaVrelMax
REAL(8), INTENT(inout):: sigmaVrelMaxNew
TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
REAL(8), INTENT(in):: vRel
REAL(8):: rMass, eRel
TYPE(particle), POINTER:: electron => NULL(), neutral => NULL()
REAL(8), DIMENSION(1:3):: vChange
TYPE(particle), POINTER:: newElectron
REAL(8):: vRel, eRel
REAL(8):: sigmaVrel
REAL(8), DIMENSION(1:3):: vp_e, vp_n
!eRel (in units of [m][L]^2[t]^-2
vRel = NORM2(part_i%v-part_j%v)
eRel = self%rMass*vRel**2
rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
eRel = rMass*vRel**2
!Relative energy must be higher than threshold
IF (eRel > self%eThreshold) THEN
sigmaVrel = self%crossSec%get(eRel)*vRel
sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
!Find which particle is the ionizing electron
IF (ASSOCIATED(part_i%species, self%electron)) THEN
electron => part_i
neutral => part_j
IF (ASSOCIATED(part_i%species, self%electron)) THEN
electron => part_i
neutral => part_j
ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN
electron => part_j
neutral => part_i
ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN
electron => part_j
neutral => part_i
ELSE
CALL criticalError("No matching between input particles and ionizing species", 'collideBinaryIonization')
END IF
!Exchange energy between
vp_e = electron%v*(1.D0 - self%deltaV/NORM2(electron%v))
vp_n = neutral%v* (1.D0 + self%deltaV/NORM2(neutral%v) )
!Changes velocity of impacting electron
electron%v = vp_e
!Creates a new electron from ionization
ALLOCATE(newElectron)
newElectron%species => electron%species
newElectron%v = vp_n
newElectron%r = neutral%r
newElectron%xi = neutral%xi
newElectron%n_in = .TRUE.
newElectron%vol = neutral%vol
newElectron%volColl = neutral%volColl
newElectron%weight = neutral%weight
newElectron%qm = electron%qm
!Ionize neutral particle
SELECT TYPE(sp => neutral%species)
TYPE IS(speciesNeutral)
CALL sp%ionize(neutral)
CLASS DEFAULT
CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization')
END SELECT
!Adds new electron to list of new particles from collisions
CALL OMP_SET_LOCK(lockCollisions)
CALL partCollisions%add(newElectron)
CALL OMP_UNSET_LOCK(lockCollisions)
ELSE
CALL criticalError("No matching between input particles and ionizing species", 'collideBinaryIonization')
END IF
!Ionize neutral particle
SELECT TYPE(sp => neutral%species)
TYPE IS(speciesNeutral)
CALL sp%ionize(neutral)
CLASS DEFAULT
! CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization')
RETURN
END SELECT
!Exchange of velocity between particles
vChange = self%deltaV*randomDirectionVHS()
!Energy is loss by the primary electron
electron%v = electron%v - vChange
!Creates a new electron from ionization
ALLOCATE(newElectron)
newElectron%species => electron%species
!Secondary electorn gains energy from ionization
newElectron%v = vChange
newElectron%r = neutral%r
newElectron%xi = neutral%xi
newElectron%n_in = .TRUE.
newElectron%vol = neutral%vol
newElectron%volColl = neutral%volColl
newElectron%weight = neutral%weight
!Adds new electron to list of new particles from collisions
CALL OMP_SET_LOCK(lockCollisions)
CALL partCollisions%add(newElectron)
CALL OMP_UNSET_LOCK(lockCollisions)
END IF
END SUBROUTINE collideBinaryIonization
!ELECTRON ION RESONANT RECOMBINATION
!Inits electron ion recombination
SUBROUTINE initBinaryRecombination(collision, crossSectionFilename, energyBinding, mass_i, mass_j, electron)
SUBROUTINE initBinaryRecombination(collision, crossSectionFilename, energyBinding, electron)
USE moduleTable
USE moduleRefParam
USE moduleConstParam
@ -362,7 +375,6 @@ MODULE moduleCollisions
CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
REAL(8), INTENT(in):: energyBinding
REAL(8), INTENT(in):: mass_i, mass_j
CHARACTER(:), ALLOCATABLE, INTENT(in):: electron
INTEGER:: electronIndex
@ -374,10 +386,6 @@ MODULE moduleCollisions
!Convert to no-dimensional units
CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
!Calculates reduced mass
collision%sMassInv = 1.D0/(mass_i+mass_j)
collision%rMass = (mass_i*mass_j)*collision%sMassInv
!Specific parameters for ionization collision
SELECT TYPE(collision)
TYPE IS(collisionBinaryRecombination)
@ -399,9 +407,8 @@ MODULE moduleCollisions
END SUBROUTINE initBinaryRecombination
!Binary electron impact ionization process
SUBROUTINE collideBinaryRecombination(self, sigmaVrelMax, sigmaVrelMaxNew, &
part_i, part_j)
!Binary recombination
SUBROUTINE collideBinaryRecombination(self, part_i, part_j, vRel)
USE moduleSpecies
USE moduleErrors
USE moduleList
@ -410,59 +417,47 @@ MODULE moduleCollisions
IMPLICIT NONE
CLASS(collisionBinaryRecombination), INTENT(in):: self
REAL(8), INTENT(in):: sigmaVrelMax
REAL(8), INTENT(inout):: sigmaVrelMaxNew
REAL(8), INTENT(in):: vRel
TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
TYPE(particle), POINTER:: electron => NULL(), ion => NULL()
REAL(8):: vRel, eRel
REAL(8):: sigmaVrel
REAL(8), DIMENSION(1:3):: vp_i
!eRel (in units of [m][L]^2[t]^-2
vRel = NORM2(part_i%v-part_j%v)
eRel = self%rMass*vRel**2
!Relative energy must be higher than threshold
sigmaVrel = self%crossSec%get(eRel)*vRel
sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
!Find which particle is the ionizing electron
IF (ASSOCIATED(part_i%species, self%electron)) THEN
electron => part_i
ion => part_j
IF (ASSOCIATED(part_i%species, self%electron)) THEN
electron => part_i
ion => part_j
ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN
electron => part_j
ion => part_i
ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN
electron => part_j
ion => part_i
ELSE
CALL criticalError("No matching between input particles and ionizing species", 'collideBinaryIonization')
END IF
!Excess energy
!TODO: This energy should be transformed into photons
vp_i = ion%v* (1.D0 - (vRel + self%deltaV)/NORM2(ion%v))
!Remove electron from simulation
electron%n_in = .FALSE.
!Neutralize ion particle
SELECT TYPE(sp => ion%species)
TYPE IS(speciesCharged)
CALL sp%neutralize(ion)
CLASS DEFAULT
CALL criticalError(sp%name // " is not a charge", 'collideBinaryRecombination')
END SELECT
ELSE
CALL criticalError("No matching between input particles and ionizing species", 'collideBinaryIonization')
END IF
!Excess energy
!TODO: This energy should be transformed into photons
vp_i = ion%v* (1.D0 - (vRel + self%deltaV)/NORM2(ion%v))
!Remove electron from simulation
electron%n_in = .FALSE.
!Neutralize ion particle
SELECT TYPE(sp => ion%species)
TYPE IS(speciesCharged)
CALL sp%neutralize(ion)
CLASS DEFAULT
CALL criticalError(sp%name // " is not a charge", 'collideBinaryRecombination')
END SELECT
END SUBROUTINE collideBinaryRecombination
!RESONANT CHARGE EXCHANGE
!Inits resonant charge exchange
SUBROUTINE initBinaryChargeExchange(collision, crossSectionFilename, mass_i, mass_j)
SUBROUTINE initBinaryChargeExchange(collision, crossSectionFilename)
USE moduleTable
USE moduleRefParam
USE moduleConstParam
@ -470,7 +465,6 @@ MODULE moduleCollisions
CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
REAL(8), INTENT(in):: mass_i, mass_j
ALLOCATE(collisionBinaryChargeExchange:: collision)
@ -480,56 +474,39 @@ MODULE moduleCollisions
!Convert to no-dimensional units
CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
!Calculates reduced mass
collision%sMassInv = 1.D0/(mass_i+mass_j)
collision%rMass = (mass_i*mass_j)/collision%sMassInv
END SUBROUTINE initBinaryChargeExchange
SUBROUTINE collideBinaryChargeExchange(self, sigmaVrelMax, sigmaVrelMaxNew, &
part_i, part_j)
SUBROUTINE collideBinaryChargeExchange(self, part_i, part_j, vRel)
USE moduleSpecies
USE moduleRandom
USE moduleMath
IMPLICIT NONE
CLASS(collisionBinaryChargeExchange), INTENT(in):: self
REAL(8), INTENT(in):: sigmaVrelMax
REAL(8), INTENT(inout):: sigmaVrelMaxNew
REAL(8), INTENT(in):: vRel
TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
REAL(8):: sigmaVrel
REAL(8):: vRel !relative velocity
REAL(8):: eRel !relative energy
!eRel (in units of [m][L]^2[t]^-2
vRel = NORM2(part_i%v-part_j%v)
eRel = self%rMass*vRel**2
sigmaVrel = self%crossSec%get(eRel)*vRel
sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
SELECT TYPE(sp => part_i%species)
TYPE IS (speciesNeutral)
!Species i is neutral, ionize particle i
CALL sp%ionize(part_i)
SELECT TYPE(sp => part_i%species)
TYPE IS (speciesNeutral)
!Species i is neutral, ionize particle i
CALL sp%ionize(part_i)
TYPE IS (speciesCharged)
!Species i is charged, neutralize particle i
CALL sp%neutralize(part_i)
TYPE IS (speciesCharged)
!Species i is charged, neutralize particle i
CALL sp%neutralize(part_i)
END SELECT
END SELECT
SELECT TYPE(sp => part_j%species)
TYPE IS (speciesNeutral)
!Species j is neutral, ionize particle j
CALL sp%ionize(part_j)
SELECT TYPE(sp => part_j%species)
TYPE IS (speciesNeutral)
!Species j is neutral, ionize particle j
CALL sp%ionize(part_j)
TYPE IS (speciesCharged)
!Species j is charged, neutralize particle j
CALL sp%neutralize(part_j)
TYPE IS (speciesCharged)
!Species j is charged, neutralize particle j
CALL sp%neutralize(part_j)
END SELECT
END IF
END SELECT
END SUBROUTINE collideBinaryChargeExchange

View file

@ -62,6 +62,23 @@ MODULE moduleEM
END FUNCTION gatherElecField
PURE FUNCTION gatherMagnField(part) RESULT(BField)
USE moduleSpecies
USE moduleMesh
IMPLICIT NONE
TYPE(particle), INTENT(in):: part
REAl(8):: xi(1:3) !Logical coordinates of particle in element
REAL(8):: BField(1:3)
BField = 0.D0
xi = part%xi
BField = mesh%vols(part%vol)%obj%gatherMF(xi)
END FUNCTION gatherMagnField
SUBROUTINE assembleSourceVector(vectorF)
USE moduleMesh
USE moduleRefParam

View file

@ -82,6 +82,7 @@ MODULE moduleInject
INTEGER, INTENT(in):: i
REAL(8), INTENT(in):: v, n(1:3), T(1:3)
INTEGER, INTENT(in):: sp, physicalSurface
REAL(8):: tauInject
REAL(8), INTENT(in):: flow
CHARACTER(:), ALLOCATABLE, INTENT(in):: units
INTEGER:: e, et
@ -93,18 +94,19 @@ MODULE moduleInject
self%n = n / NORM2(n)
self%T = T / T_ref
self%species => species(sp)%obj
tauInject = tau(self%species%n)
SELECT CASE(units)
CASE ("sccm")
!Standard cubic centimeter per minute
self%nParticles = INT(flow*sccm2atomPerS*tauMin*ti_ref/species(sp)%obj%weight)
self%nParticles = INT(flow*sccm2atomPerS*tauInject*ti_ref/species(sp)%obj%weight)
CASE ("A")
!Input current in Ampers
self%nParticles = INT(flow*tauMin*ti_ref/(qe*species(sp)%obj%weight))
self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight))
CASE ("part/s")
!Input current in Ampers
self%nParticles = INT(flow*tauMin*ti_ref/species(sp)%obj%weight)
self%nParticles = INT(flow*tauInject*ti_ref/species(sp)%obj%weight)
CASE DEFAULT
CALL criticalError("No support for units: " // units, 'initInject')
@ -154,6 +156,7 @@ MODULE moduleInject
!Calculates cumulative probability
ALLOCATE(self%cumWeight(1:self%nEdges))
et = 1
self%cumWeight(1) = mesh%edges(self%edges(et))%obj%weight
DO et = 2, self%nEdges
self%cumWeight(et) = mesh%edges(self%edges(et))%obj%weight + self%cumWeight(et-1)
@ -161,8 +164,6 @@ MODULE moduleInject
END DO
self%sumWeight = self%cumWeight(self%nEdges)
nPartInj = nPartInj + self%nParticles
END SUBROUTINE initInject
!Injection of particles
@ -174,12 +175,24 @@ MODULE moduleInject
INTEGER:: i
!$OMP SINGLE
nPartInj = 0
DO i = 1, nInject
IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN
nPartInj = nPartInj + inject(i)%nParticles
END IF
END DO
IF (ALLOCATED(partInj)) DEALLOCATE(partInj)
ALLOCATE(partInj(1:nPartInj))
!$OMP END SINGLE
DO i=1, nInject
CALL inject(i)%addParticles()
IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN
CALL inject(i)%addParticles()
END IF
END DO
END SUBROUTINE doInjects
@ -240,23 +253,27 @@ MODULE moduleInject
CLASS(injectGeneric), INTENT(in):: self
INTEGER:: randomX
INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
INTEGER:: n
INTEGER:: i
INTEGER:: n, sp
CLASS(meshEdge), POINTER:: randomEdge
!Insert particles
!$OMP SINGLE
nMin = SUM(inject(1:(self%id-1))%nParticles) + 1
nMin = 0
DO i = 1, self%id -1
IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN
nMin = nMin + inject(i)%nParticles
END IF
END DO
nMin = nMin + 1
nMax = nMin + self%nParticles - 1
!Assign weight to particle.
partInj(nMin:nMax)%weight = self%species%weight
!Particle is considered to be outside the domain
partInj(nMin:nMax)%n_in = .FALSE.
!Assign charge/mass to charged particle.
SELECT TYPE(sp => self%species)
TYPE IS(speciesCharged)
partInj(nMin:nMax)%qm = sp%qm
END SELECT
!$OMP END SINGLE
!$OMP DO
@ -278,6 +295,7 @@ MODULE moduleInject
END IF
partInj(n)%volColl = randomEdge%eColl%n
sp = self%species%n
!Assign particle type
partInj(n)%species => self%species
@ -289,7 +307,7 @@ MODULE moduleInject
!Obtain natural coordinates of particle in cell
partInj(n)%xi = mesh%vols(partInj(n)%vol)%obj%phy2log(partInj(n)%r)
!Push new particle with the minimum time step
CALL solver%pusher(self%species%n)%pushParticle(partInj(n), tauMin)
CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
!Assign cell to new particle
CALL solver%updateParticleCell(partInj(n))

View file

@ -40,16 +40,6 @@ MODULE moduleInput
CALL readSpecies(config)
CALL checkStatus(config, "readSpecies")
!Reads case parameters
CALL verboseError('Reading Case parameters...')
CALL readCase(config)
CALL checkStatus(config, "readCase")
!Read interactions between species
CALL verboseError('Reading interaction between species...')
CALL readInteractions(config)
CALL checkStatus(config, "readInteractions")
!Read boundaries
CALL verboseError('Reading boundary conditions...')
CALL readBoundary(config)
@ -60,6 +50,16 @@ MODULE moduleInput
CALL readGeometry(config)
CALL checkStatus(config, "readGeometry")
!Read solver parameters
CALL verboseError('Reading Solver parameters...')
CALL readSolver(config)
CALL checkStatus(config, "readSolver")
!Read interactions between species
CALL verboseError('Reading interaction between species...')
CALL readInteractions(config)
CALL checkStatus(config, "readInteractions")
!Read probes
CALL verboseError('Reading Probes...')
CALL readProbes(config)
@ -85,9 +85,12 @@ MODULE moduleInput
!Copies input file to output folder
CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder)
!Copies particle mesh
CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder)
IF (doubleMesh) THEN
CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder)
IF (mesh%dimen > 0) THEN
CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder)
IF (doubleMesh) THEN
CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder)
END IF
END IF
@ -129,22 +132,6 @@ MODULE moduleInput
CALL config%get(object // '.temperature', T_ref, found)
IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference')
!If a reference cross section is given, it is used
CALL config%get(object // '.crossSection', sigma_ref, found)
!If not, the reference radius is searched
IF (.NOT. found) THEN
CALL config%get(object // '.radius', r_ref, found)
IF (found) THEN
sigma_ref = PI*(r_ref+r_ref)**2 !reference cross section
ELSE
sigma_ref = 0.D0 !Assume no collisions
END IF
END IF
!Derived parameters
L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !reference length
v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity
@ -152,11 +139,28 @@ MODULE moduleInput
Vol_ref = L_ref**3 !reference volume
EF_ref = qe*n_ref*L_ref/eps_0 !reference electric field
Volt_ref = EF_ref*L_ref !reference voltage
B_ref = m_ref / (ti_ref * qe) !reference magnetic field
!If a reference cross section is given, it is used
CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found)
!If not, the reference radius is searched
IF (.NOT. found) THEN
CALL config%get(object // '.radius', r_ref, found)
IF (found) THEN
sigmaVrel_ref = PI*(r_ref+r_ref)**2*v_ref !reference cross section times velocity
ELSE
sigmaVrel_ref = 0.D0 !Assume no collisions
END IF
END IF
END SUBROUTINE readReference
!Reads the specific case parameters
SUBROUTINE readCase(config)
SUBROUTINE readSolver(config)
USE moduleRefParam
USE moduleErrors
USE moduleCaseParam
@ -164,6 +168,7 @@ MODULE moduleInput
USE moduleSpecies
USE moduleCollisions
USE moduleOutput
USE moduleMesh
USE json_module
IMPLICIT NONE
@ -172,47 +177,68 @@ MODULE moduleInput
CHARACTER(:), ALLOCATABLE:: object
!simulation final and initial times in [t]
REAL(8):: finalTime, initialTime
CHARACTER(:), ALLOCATABLE:: pusherType, WSType
CHARACTER(:), ALLOCATABLE:: pusherType, WSType, EMType
REAL(8):: B(1:3)
INTEGER:: nTau, nSolver
INTEGER:: i
INTEGER:: i, n
CHARACTER(2):: iString
CHARACTER(1):: tString
object = 'case'
object = 'solver'
!Time parameters
CALL config%info(object // '.tau', found, n_children = nTau)
IF (.NOT. found .OR. nTau == 0) CALL criticalError('Required parameter tau not found','readCase')
IF (.NOT. found .OR. nTau == 0) THEN
CALL criticalError('Required parameter tau not found','readSolver')
END IF
ALLOCATE(tau(1:nSpecies))
DO i = 1, nTau
WRITE(iString, '(I2)') i
CALL config%get(object // '.tau(' // TRIM(iString) // ')', tau(i), found)
END DO
tauMin = MINVAL(tau(1:nTau))
IF (nTau < nSpecies) THEN
CALL warningError('Using minimum time step for some species')
tau(nTau+1:nSpecies) = MINVAL(tau(1:nTau))
tau(nTau+1:nSpecies) = tauMin
END IF
tauMin = MINVAL(tau)
!Gets the simulation final time
CALL config%get(object // '.finalTime', finalTime, found)
IF (.NOT. found) CALL criticalError('Required parameter finalTime not found','readCase')
IF (.NOT. found) THEN
CALL criticalError('Required parameter finalTime not found','readSolver')
END IF
!Convert simulation time to number of iterations
tFinal = INT(finalTime / tauMin)
!Gets the simulation initial time
CALL config%get(object // '.initialTime', initialTime, found)
IF (found) tInitial = INT(initialTime / tauMin)
IF (found) THEN
tInitial = INT(initialTime / tauMin)
END IF
!Gest the pusher for each species
CALL config%info(object // '.pusher', found, n_children = nSolver)
IF (.NOT. found .OR. nSolver /= nSpecies) CALL criticalError('Required parameter pusher not found','readCase')
IF (mesh%dimen > 0) THEN
CALL config%info(object // '.pusher', found, n_children = nSolver)
IF (.NOT. found .OR. nSolver /= nSpecies) THEN
CALL criticalError('Required parameter pusher not found','readSolver')
END IF
END IF
!Allocates all the pushers for particles
ALLOCATE(solver%pusher(1:nSpecies))
!Initialize pushers
DO i = 1, nSolver
DO i = 1, nSpecies
WRITE(iString, '(I2)') i
CALL config%get(object // '.pusher(' // TRIM(iString) // ')', pusherType, found)
@ -220,23 +246,64 @@ MODULE moduleInput
CALL solver%pusher(i)%init(pusherType, tauMin, tau(i))
END DO
!Gest the non-analogue scheme
!Get the non-analogue scheme
CALL config%get(object // '.WeightingScheme', WSType, found)
CALL solver%initWS(WSType)
!Makes tau(s) non-dimensional
!Make tau(s) non-dimensional
tau = tau / ti_ref
tauMin = tauMin / ti_ref
!Sets the format of output files accordint to iteration number
!Set the format of output files accordint to iteration number
iterationDigits = INT(LOG10(REAL(tFinal))) + 1
WRITE(tString, '(I1)') iterationDigits
iterationFormat = "(I" // tString // "." // tString // ")"
END SUBROUTINE readCase
!Get EM solver
CALL config%get(object // '.EMSolver', EMType, found)
IF (found) THEN
CALL solver%initEM(EMType)
SELECT CASE(EMType)
CASE("Electrostatic")
!Read BC
CALL readEMBoundary(config)
!Reads the initial information for the species
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)
IF (.NOT. found) THEN
CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver')
END IF
END DO
!Non-dimensional magnetic field
B = B / B_ref
!Assign it to the nodes
DO n =1, mesh%numNodes
mesh%nodes(n)%obj%emData%B(1) = B(1)
mesh%nodes(n)%obj%emData%B(2) = B(2)
mesh%nodes(n)%obj%emData%B(3) = B(3)
END DO
CASE DEFAULT
CALL criticalError('EM Solver ' // EMType // ' not found', 'readSolver')
END SELECT
END IF
END SUBROUTINE readSolver
!Read the initial information for the species
SUBROUTINE readInitial(config)
USE moduleSpecies
USE moduleMesh
@ -266,17 +333,18 @@ MODULE moduleInput
!Mean velocity and temperature at particle position
REAL(8):: velocityXi(1:3), temperatureXi
INTEGER:: nNewPart = 0.D0
CLASS(meshVol), POINTER:: vol
TYPE(particle), POINTER:: partNew
REAL(8):: vTh
TYPE(lNode), POINTER:: partCurr, partNext
CALL config%info('case.initial', found, n_children = nInitial)
CALL config%info('solver.initial', found, n_children = nInitial)
IF (found) THEN
!Reads the information from initial species
DO i = 1, nInitial
WRITE(iString, '(I2)') i
object = 'case.initial(' // iString // ')'
object = 'solver.initial(' // iString // ')'
CALL config%get(object // '.species', spName, found)
sp = speciesName2Index(spName)
CALL config%get(object // '.file', spFile, found)
@ -289,7 +357,6 @@ MODULE moduleInput
!Density at centroid of cell
nodes = mesh%vols(e)%obj%getNodes()
nNodes = SIZE(nodes)
!TODO: Procedure to obtain centroid from element (also for printing Electric Field)
fPsi = mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/))
ALLOCATE(source(1:nNodes))
DO j = 1, nNodes
@ -350,19 +417,17 @@ MODULE moduleInput
END IF
partNew%n_in = .TRUE.
partNew%weight = species(sp)%obj%weight
!If charged species, add qm to particle
SELECT TYPE(sp => species(sp)%obj)
TYPE IS (speciesCharged)
partNew%qm = sp%qm
CLASS DEFAULT
partNew%qm = 0.D0
END SELECT
!Assign particle to temporal list of particles
CALL partInitial%add(partNew)
!Assign particle to list in volume
vol => meshforMCC%vols(partNew%volColl)%obj
CALL OMP_SET_LOCK(vol%lock)
CALL vol%listPart_in(sp)%add(partNew)
vol%totalWeight(sp) = vol%totalWeight(sp) + partNew%weight
CALL OMP_UNSET_LOCK(vol%lock)
END DO
DEALLOCATE(source)
@ -372,7 +437,7 @@ MODULE moduleInput
END DO
!Convert temporal list of particles into initial partOld array
!Deallocates the list of initial particles
!Deallocate the list of initial particles
nNewPart = partInitial%amount
IF (nNewPart > 0) THEN
ALLOCATE(partOld(1:nNewPart))
@ -504,13 +569,13 @@ MODULE moduleInput
END DO
!Reads relations between species
!Read relations between species
DO i = 1, nSpecies
WRITE(iString, '(I2)') i
object = 'species(' // TRIM(iString) // ')'
SELECT TYPE(sp => species(i)%obj)
TYPE IS (speciesNeutral)
!Gets species linked ion
!Get species linked ion
CALL config%get(object // '.ion', linkName, found)
IF (found) THEN
linkID = speciesName2Index(linkName)
@ -519,7 +584,7 @@ MODULE moduleInput
END IF
TYPE IS (speciesCharged)
!Gets species linked neutral
!Get species linked neutral
CALL config%get(object // '.neutral', linkName, found)
IF (found) THEN
linkID = speciesName2Index(linkName)
@ -568,11 +633,13 @@ MODULE moduleInput
CHARACTER(:), ALLOCATABLE:: crossSecFilePath
CHARACTER(:), ALLOCATABLE:: cType
LOGICAL:: found
INTEGER:: nInteractions, nCollisions
INTEGER:: nPairs, nCollisions
INTEGER:: i, k, ij
INTEGER:: pt_i, pt_j
REAL(8):: energyThreshold, energyBinding
CHARACTER(:), ALLOCATABLE:: electron
INTEGER:: e
CLASS(meshVol), POINTER:: vol
!Firstly, checks if the object 'interactions' exists
CALL config%info('interactions', found)
@ -580,19 +647,6 @@ MODULE moduleInput
!Checks if MC collisions have been defined
CALL config%info('interactions.collisions', found)
IF (found) THEN
!Checks if a mesh for collisions has been defined
!The mesh will be initialized and reader in readGeometry
CALL config%info('interactions.meshCollisions', found)
IF (found) THEN
!Points meshForMCC to the specific mesh defined
meshForMCC => meshColl
ELSE
!Points the meshForMCC pointer to the Particles Mesh
meshForMCC => mesh
END IF
!Reads collision time step
CALL config%info('interactions.timeStep', found)
IF (found) THEN
@ -620,8 +674,8 @@ MODULE moduleInput
!Inits lock for list of particles
CALL OMP_INIT_LOCK(lockCollisions)
CALL config%info('interactions.collisions', found, n_children = nInteractions)
DO i = 1, nInteractions
CALL config%info('interactions.collisions', found, n_children = nPairs)
DO i = 1, nPairs
WRITE(iString, '(I2)') i
object = 'interactions.collisions(' // TRIM(iString) // ')'
CALL config%get(object // '.species_i', species_i, found)
@ -630,8 +684,9 @@ MODULE moduleInput
pt_j = speciesName2Index(species_j)
CALL config%info(object // '.cTypes', found, n_children = nCollisions)
ij = interactionIndex(pt_i,pt_j)
!Allocates the required number of collisions per each pair of species ij
CALL interactionMatrix(ij)%init(nCollisions)
CALL interactionMatrix(ij)%init(nCollisions, pt_i, pt_j)
DO k = 1, nCollisions
WRITE (kString, '(I2)') k
@ -646,13 +701,11 @@ MODULE moduleInput
SELECT CASE(cType)
CASE ('elastic')
!Elastic collision
CALL initBinaryElastic(interactionMatrix(ij)%collisions(k)%obj, &
crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
CALL initBinaryElastic(interactionMatrix(ij)%collisions(k)%obj, crossSecFilePath)
CASE ('chargeExchange')
!Resonant charge exchange
CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, &
crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, crossSecFilePath)
CASE ('ionization')
!Electorn impact ionization
@ -661,7 +714,7 @@ MODULE moduleInput
CALL config%get(object // '.electron', electron, found)
IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
crossSecFilePath, energyThreshold, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)
crossSecFilePath, energyThreshold, electron)
CASE ('recombination')
!Electorn impact ionization
@ -670,10 +723,10 @@ MODULE moduleInput
CALL config%get(object // '.electron', electron, found)
IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
CALL initBinaryRecombination(interactionMatrix(ij)%collisions(k)%obj, &
crossSecFilePath, energyBinding, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)
crossSecFilePath, energyBinding, electron)
CASE DEFAULT
CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions')
CALL criticalError('Collision type' // cType // 'not defined', 'readInteractions')
END SELECT
@ -681,6 +734,26 @@ MODULE moduleInput
END DO
!Init the required arrays in each volume to account for MCC.
DO e = 1, meshForMCC%numVols
vol => meshForMCC%vols(e)%obj
!Allocate Maximum cross section per collision pair and assign the initial collision rate
ALLOCATE(vol%sigmaVrelMax(1:nCollPairs))
vol%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
IF (collOutput) THEN
ALLOCATE(vol%tallyColl(1:nCollPairs))
DO k = 1, nCollPairs
ALLOCATE(vol%tallyColl(k)%tally(1:interactionmatrix(k)%amount))
vol%tallyColl(k)%tally = 0
END DO
END IF
END DO
END IF
END IF
@ -792,11 +865,6 @@ MODULE moduleInput
USE moduleMesh
USE moduleMeshInputGmsh2, ONLY: initGmsh2
USE moduleMeshInput0D, ONLY: init0D
USE moduleMesh3DCart, ONLY: connectMesh3DCart
USE moduleMesh2DCyl, ONLY: connectMesh2DCyl
USE moduleMesh2DCart, ONLY: connectMesh2DCart
USE moduleMesh1DRad, ONLY: connectMesh1DRad
USE moduleMesh1DCart, ONLY: connectMesh1DCart
USE moduleErrors
USE moduleOutput
USE moduleRefParam
@ -805,86 +873,121 @@ MODULE moduleInput
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile, EMType
CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
REAL(8):: volume
!Firstly, indicates if a specific mesh for MC collisions is being use
doubleMesh = ASSOCIATED(meshForMCC, meshColl)
object = 'geometry'
!Selects the type of geometry.
CALL config%get('geometry.type', mesh%geometry, found)
IF (doubleMesh) meshColl%geometry = mesh%geometry
!Checks if a mesh for collisions has been defined
!The mesh will be initialized and reader in readGeometry
CALL config%info(object // '.meshCollisions', found)
IF (found) THEN
!Points meshForMCC to the specific mesh defined
meshForMCC => meshColl
doubleMesh = .TRUE.
!Gets the type of mesh
CALL config%get('geometry.meshType', meshFormat, found)
SELECT CASE(meshFormat)
CASE ("gmsh2")
CALL initGmsh2(mesh)
IF (doubleMesh) CALL initGmsh2(meshColl)
ELSE
CALL config%info('interactions', found)
IF (found) THEN
!Points the meshForMCC pointer to the Particles Mesh
meshForMCC => mesh
CASE ("0D")
CALL config%get('geometry.meshType', meshFormat, found)
CALL init0D(mesh)
CASE DEFAULT
CALL criticalError("Mesh format " // meshFormat // " not recogniced", "readGeometry")
END SELECT
!Reads the mesh file
CALL config%get('geometry.meshFile', meshFile, found)
pathMeshParticle = path // meshFile
CALL mesh%readMesh(pathMeshParticle)
DEALLOCATE(meshFile)
IF (doubleMesh) THEN
!Reads the mesh file for collisions
CALL config%get('interactions.meshCollisions', meshFile, found)
pathMeshColl = path // meshFile
CALL meshColl%readMesh(pathMeshColl)
END IF
doubleMesh = .FALSE.
END IF
!Gets the volume for a 0D mesh
!TODO: Try to constrain this to the inout for 0D
IF (meshFormat == "0D") THEN
CALL config%get('geometry.volume', volume, found)
!Get the dimension of the geometry
CALL config%get(object // '.dimension', mesh%dimen, found)
IF (.NOT. found) THEN
CALL criticalError('Required parameter geometry.dimension not found', 'readGeometry')
END IF
IF (doubleMesh) THEN
meshColl%dimen = mesh%dimen
END IF
SELECT CASE(mesh%dimen)
CASE (0)
CALL init0D(mesh)
!Read the 0D mesh
CALL mesh%readMesh(pathMeshParticle)
!Get the volumne
CALL config%get(object // '.volume', volume, found)
!Rescale the volumne
IF (found) THEN
mesh%vols(1)%obj%volume = mesh%vols(1)%obj%volume*volume / Vol_ref
mesh%nodes(1)%obj%v = mesh%vols(1)%obj%volume
END IF
END IF
CASE (1:3)
!Select the type of geometry
CALL config%get(object // '.type', mesh%geometry, found)
IF (doubleMesh) THEN
meshColl%geometry = mesh%geometry
!Creates the connectivity between elements
SELECT CASE(mesh%geometry)
CASE("3DCart")
mesh%connectMesh => connectMesh3DCart
END IF
!Check the consistency between geometry and dimension
SELECT CASE(mesh%dimen)
CASE(3)
IF (mesh%geometry /= 'Cart') THEN
CALL criticalError('No available geometry ' // mesh%geometry // ' in 3D', 'readGeometry')
CASE("2DCyl")
mesh%connectMesh => connectMesh2DCyl
END IF
CASE("2DCart")
mesh%connectMesh => connectMesh2DCart
CASE(2)
IF (mesh%geometry /= 'Cart' .AND. &
mesh%geometry /= 'Cyl') THEN
CALL criticalError('No available geometry ' // mesh%geometry // ' in 2D', 'readGeometry')
CASE("1DRad")
mesh%connectMesh => connectMesh1DRad
END IF
CASE("1DCart")
mesh%connectMesh => connectMesh1DCart
CASE(1)
IF (mesh%geometry /= 'Cart' .AND. &
mesh%geometry /= 'Rad') THEN
CALL criticalError('No available geometry ' // mesh%geometry // ' in 1D', 'readGeometry')
CASE("0D")
mesh%connectMesh => NULL()
END IF
END SELECT
!Get the type of mesh
CALL config%get(object // '.meshType', meshFormat, found)
SELECT CASE(meshFormat)
CASE ("gmsh2")
CALL initGmsh2(mesh)
IF (doubleMesh) THEN
CALL initGmsh2(meshColl)
END IF
END SELECT
!Reads the mesh file
CALL config%get(object // '.meshFile', meshFile, found)
pathMeshParticle = path // meshFile
CALL mesh%readMesh(pathMeshParticle)
DEALLOCATE(meshFile)
IF (doubleMesh) THEN
!Reads the mesh file for collisions
CALL config%get(object // '.meshCollisions', meshFile, found)
pathMeshColl = path // meshFile
CALL meshColl%readMesh(pathMeshColl)
END IF
CASE DEFAULT
CALL criticalError("Mesh dimension not recogniced", "readGeometry")
END SELECT
IF (ASSOCIATED(mesh%connectMesh)) CALL mesh%connectMesh
IF (doubleMesh) THEN
meshColl%connectMesh => mesh%connectMesh
CALL meshColl%connectMesh
END IF
!Builds the K matrix for the Particles mesh
CALL mesh%constructGlobalK()
@ -898,16 +1001,6 @@ MODULE moduleInput
END IF
!Gest EM solver
CALL config%get('case.EMSolver', EMType, found)
CALL solver%initEM(EMType)
SELECT CASE(EMType)
CASE("Electrostatic")
!Reads BC
CALL readEMBoundary(config)
END SELECT
END SUBROUTINE readGeometry
SUBROUTINE readProbes(config)
@ -1010,7 +1103,6 @@ MODULE moduleInput
END DO
!TODO: Improve this
IF (ALLOCATED(boundEM)) THEN
DO e = 1, mesh%numEdges
IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN
@ -1026,7 +1118,10 @@ MODULE moduleInput
!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)
IF (info /= 0) CALL criticalError('Factorization of K matrix failed', 'readEMBoundary')
IF (info /= 0) THEN
CALL criticalError('Factorization of K matrix failed', 'readEMBoundary')
END IF
END SUBROUTINE readEMBoundary

View file

@ -71,6 +71,7 @@ MODULE moduleList
DO n=1, self%amount
!Point element in array to element in list
partArray(n)%part => tempNode%part
!Go to next element
tempNode => tempNode%next

View file

@ -49,4 +49,14 @@ MODULE moduleMath
END FUNCTION norm1
PURE FUNCTION reducedMass(m_i, m_j) RESULT(rMass)
IMPLICIT NONE
REAL(8), INTENT(in):: m_i, m_j
REAL(8):: rMass
rMass = (m_i * m_j) / (m_i + m_j)
END FUNCTION
END MODULE moduleMath

View file

@ -11,6 +11,7 @@ MODULE moduleOutput
TYPE emNode
CHARACTER(:), ALLOCATABLE:: type
REAL(8):: phi
REAL(8):: B(1:3)
END TYPE emNode

View file

@ -1,9 +1,9 @@
!Reference parameters
MODULE moduleRefParam
!Parameters that define the problem (inputs)
REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigma_ref
REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigmaVrel_ref
!Reference parameters for non-dimensional problem
REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref
REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref
END MODULE moduleRefParam

View file

@ -60,53 +60,63 @@ MODULE moduleSolver
!Init Pusher
SUBROUTINE initPusher(self, pusherType, tau, tauSp)
USE moduleErrors
USE moduleMesh, ONLY: mesh
IMPLICIT NONE
CLASS(pusherGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE:: pusherType
REAL(8):: tau, tauSp
SELECT CASE(pusherType)
!3D Cartesian
CASE('3DCartNeutral')
self%pushParticle => push3DCartNeutral
CASE('3DCartCharged')
self%pushParticle => push3DCartCharged
!2D Cylindrical
CASE('2DCylNeutral')
self%pushParticle => push2DCylNeutral
CASE('2DCylCharged')
self%pushParticle => push2DCylCharged
!2D Cartesian
CASE('2DCartNeutral')
self%pushParticle => push2DCartNeutral
CASE('2DCartCharged')
self%pushParticle => push2DCartCharged
!1D Cartesian
CASE('1DCartNeutral')
self%pushParticle => push1DCartNeutral
CASE('1DCartCharged')
self%pushParticle => push1DCartCharged
!1D Radial
CASE('1DRadNeutral')
self%pushParticle => push1DRadNeutral
CASE('1DRadCharged')
self%pushParticle => push1DRadCharged
CASE('0D')
!TODO: Reorganize if Cart pushers are combined
SELECT CASE(mesh%dimen)
CASE(0)
self%pushParticle => push0D
CASE DEFAULT
CALL criticalError('Pusher ' // pusherType // ' not found','initPusher')
CASE(1:3)
SELECT CASE(mesh%geometry)
CASE ('Cart')
SELECT CASE(pusherType)
CASE('Neutral')
self%pushParticle => pushCartNeutral
CASE('Electrostatic')
self%pushParticle => pushCartElectrostatic
CASE('Electromagnetic')
self%pushParticle => pushCartElectromagnetic
CASE DEFAULT
CALL criticalError('Pusher ' // pusherType // ' not found for Cart','initPusher')
END SELECT
CASE('Cyl')
SELECT CASE(pusherType)
CASE('Neutral')
self%pushParticle => push2DCylNeutral
CASE('Electrostatic')
self%pushParticle => push2DCylElectrostatic
CASE DEFAULT
CALL criticalError('Pusher ' // pusherType // ' not found for Cyl','initPusher')
END SELECT
CASE('Rad')
SELECT CASE(pusherType)
CASE('Neutral')
self%pushParticle => push1DRadNeutral
CASE('Electrostatic')
self%pushParticle => push1DRadElectrostatic
CASE DEFAULT
CALL criticalError('Pusher ' // pusherType // ' not found for Rad','initPusher')
END SELECT
END SELECT
END SELECT
@ -122,7 +132,7 @@ MODULE moduleSolver
CHARACTER(:), ALLOCATABLE:: EMType
SELECT CASE(EMType)
CASE('Electrostatic')
CASE('Electrostatic','ConstantB')
self%solveEM => solveElecField
END SELECT
@ -155,7 +165,7 @@ MODULE moduleSolver
INTEGER:: sp
!$OMP DO
DO n=1, nPartOld
DO n = 1, nPartOld
!Select species type
sp = partOld(n)%species%n
!Checks if the species sp is update this iteration
@ -172,67 +182,76 @@ MODULE moduleSolver
END SUBROUTINE doPushes
!Push neutral particles in 3D cartesian coordinates
PURE SUBROUTINE push3DCartNeutral(part, tauIn)
!Push neutral particles in cartesian coordinates
PURE SUBROUTINE pushCartNeutral(part, tauIn)
USE moduleSPecies
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
part_temp = part
part%r = part%r + part%v*tauIn
!x
part_temp%v(1) = part%v(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
END SUBROUTINE pushCartNeutral
!y
part_temp%v(2) = part%v(2)
part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn
!z
part_temp%v(3) = part%v(3)
part_temp%r(3) = part%r(3) + part_temp%v(3)*tauIn
part_temp%n_in = .FALSE.
part = part_temp
END SUBROUTINE push3DCartNeutral
!Push charged particles in 3D cartesian coordinates
PURE SUBROUTINE push3DCartCharged(part, tauIn)
PURE SUBROUTINE pushCartElectrostatic(part, tauIn)
USE moduleSPecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
REAL(8):: qmEFt(1:3)
part_temp = part
!Get the electric field at particle position
qmEFt = part_temp%qm*gatherElecField(part_temp)*tauIn
qmEFt = part%species%qm*gatherElecField(part)*tauIn
!x
part_temp%v(1) = part%v(1) + qmEFt(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
!Update velocity
part%v = part%v + qmEFt
!y
part_temp%v(2) = part%v(2) + qmEFt(2)
part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn
!Update position
part%r = part%r + part%v*tauIn
!z
part_temp%v(3) = part%v(3) + qmEFt(3)
part_temp%r(3) = part%r(3) + part_temp%v(3)*tauIn
END SUBROUTINE pushCartElectrostatic
part_temp%n_in = .FALSE.
PURE SUBROUTINE pushCartElectromagnetic(part, tauIn)
USE moduleSPecies
USE moduleEM
USE moduleMath
IMPLICIT NONE
part = part_temp
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
REAL(8):: tauInHalf
REAL(8):: qmEFt(1:3)
REAL(8):: B(1:3), BNorm
REAL(8):: fn
REAL(8):: v_minus(1:3), v_prime(1:3), v_plus(1:3)
END SUBROUTINE push3DCartCharged
tauInHalf = tauIn *0.5D0
!Half of the force o f the electric field
qmEFt = part%species%qm*gatherElecField(part)*tauInHalf
!Half step for electrostatic
v_minus = part%v + qmEFt
!Full step rotation
B = gatherMagnField(part)
BNorm = NORM2(B)
IF (BNorm > 0.D0) THEN
fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm
v_prime = v_minus + fn * crossProduct(v_minus, B)
v_plus = v_minus + 2.D0 * fn / (1.D0 + fn**2 * B**2)*crossProduct(v_prime, B)
END IF
!Half step for electrostatic
part%v = v_plus + qmEFt
!Update position
part%r = part%r + part%v*tauIn
END SUBROUTINE pushCartElectromagnetic
!Push one particle. Boris pusher for 2D Cyl Neutral particle
PURE SUBROUTINE push2DCylNeutral(part, tauIn)
@ -265,14 +284,14 @@ MODULE moduleSolver
END IF
part_temp%v(2) = cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
part_temp%n_in = .FALSE. !Assume particle is outside until cell is found
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push2DCylNeutral
!Push one particle. Boris pusher for 2D Cyl Charged particle
PURE SUBROUTINE push2DCylCharged(part, tauIn)
!Push one particle. Boris pusher for 2D Cyl Electrostatic particle
PURE SUBROUTINE push2DCylElectrostatic(part, tauIn)
USE moduleSpecies
USE moduleEM
IMPLICIT NONE
@ -286,7 +305,7 @@ MODULE moduleSolver
part_temp = part
!Get electric field at particle position
qmEFt = part_temp%qm*gatherElecField(part_temp)*tauIn
qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauIn
!z
part_temp%v(1) = part%v(1) + qmEFt(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
@ -306,112 +325,11 @@ MODULE moduleSolver
END IF
part_temp%v(2) = cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
part_temp%n_in = .FALSE. !Assume particle is outside until cell is found
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push2DCylCharged
!Push neutral particles in 2D cartesian coordinates
PURE SUBROUTINE push2DCartNeutral(part, tauIn)
USE moduleSPecies
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
part_temp = part
!x
part_temp%v(1) = part%v(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
!y
part_temp%v(2) = part%v(2)
part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn
part_temp%n_in = .FALSE.
part = part_temp
END SUBROUTINE push2DCartNeutral
!Push charged particles in 2D cartesian coordinates
PURE SUBROUTINE push2DCartCharged(part, tauIn)
USE moduleSPecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
REAL(8):: qmEFt(1:3)
part_temp = part
!Get the electric field at particle position
qmEFt = part_temp%qm*gatherElecField(part_temp)*tauIn
!x
part_temp%v(1) = part%v(1) + qmEFt(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
!y
part_temp%v(2) = part%v(2) + qmEFt(2)
part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn
part_temp%n_in = .FALSE.
part = part_temp
END SUBROUTINE push2DCartCharged
!Push neutral particles in 1D cartesian coordinates
PURE SUBROUTINE push1DCartNeutral(part, tauIn)
USE moduleSPecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
part_temp = part
!x
part_temp%v(1) = part%v(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
part_temp%n_in = .FALSE.
part = part_temp
END SUBROUTINE push1DCartNeutral
!Push charged particles in 1D cartesian coordinates
PURE SUBROUTINE push1DCartCharged(part, tauIn)
USE moduleSPecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
REAL(8):: qmEFt(1:3)
part_temp = part
!Get the electric field at particle position
qmEFt = part_temp%qm*gatherElecField(part_temp)*tauIn
!x
part_temp%v(1) = part%v(1) + qmEFt(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
part_temp%n_in = .FALSE.
part = part_temp
END SUBROUTINE push1DCartCharged
END SUBROUTINE push2DCylElectrostatic
!Push one particle. Boris pusher for 1D Radial Neutral particle
PURE SUBROUTINE push1DRadNeutral(part, tauIn)
@ -442,14 +360,14 @@ MODULE moduleSolver
END IF
part_temp%v(1) = cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)
part_temp%n_in = .FALSE. !Assume particle is outside until cell is found
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push1DRadNeutral
!Push one particle. Boris pusher for 1D Radial Charged particle
PURE SUBROUTINE push1DRadCharged(part, tauIn)
!Push one particle. Boris pusher for 1D Radial Electrostatic particle
PURE SUBROUTINE push1DRadElectrostatic(part, tauIn)
USE moduleSpecies
USE moduleEM
IMPLICIT NONE
@ -463,7 +381,7 @@ MODULE moduleSolver
part_temp = part
!Get electric field at particle position
qmEFt = part_temp%qm*gatherElecField(part_temp)*tauMin
qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauMin
!r,theta
v_p_oh_star(1) = part%v(1) + qmEFt(1)
x_new = part%r(1) + v_p_oh_star(1)*tauIn
@ -480,11 +398,11 @@ MODULE moduleSolver
END IF
part_temp%v(1) = cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)
part_temp%n_in = .FALSE. !Assume particle is outside until cell is found
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push1DRadCharged
END SUBROUTINE push1DRadElectrostatic
!Dummy pusher for 0D geometry
PURE SUBROUTINE push0D(part, tauIn)
@ -539,6 +457,7 @@ MODULE moduleSolver
INTEGER, SAVE:: nPartNew
INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions, nSurfaces
TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:)
INTEGER:: s
!$OMP SECTIONS
!$OMP SECTION
@ -622,18 +541,30 @@ MODULE moduleSolver
END DO
!$OMP SECTION
!Erase the list of particles inside the cell
DO e = 1, mesh%numVols
mesh%vols(e)%obj%totalWeight = 0.D0
CALL mesh%vols(e)%obj%listPart_in%erase()
!Erase the list of particles inside the cell if particles have been pushed
DO s = 1, nSpecies
DO e = 1, mesh%numVols
IF (solver%pusher(s)%pushSpecies) THEN
CALL mesh%vols(e)%obj%listPart_in(s)%erase()
mesh%vols(e)%obj%totalWeight(s) = 0.D0
END IF
END DO
END DO
!$OMP SECTION
!Erase the list of particles inside the cell in coll mesh
DO e = 1, meshColl%numVols
meshColl%vols(e)%obj%totalWeight = 0.D0
CALL meshColl%vols(e)%obj%listPart_in%erase()
DO s = 1, nSpecies
DO e = 1, meshColl%numVols
IF (solver%pusher(s)%pushSpecies) THEN
CALL meshColl%vols(e)%obj%listPart_in(s)%erase()
meshColl%vols(e)%obj%totalWeight(s) = 0.D0
END IF
END DO
END DO
@ -655,7 +586,7 @@ MODULE moduleSolver
!Loops over the particles to scatter them
!$OMP DO
DO n=1, nPartOld
DO n = 1, nPartOld
CALL mesh%vols(partOld(n)%vol)%obj%scatter(partOld(n))
END DO
@ -767,6 +698,7 @@ MODULE moduleSolver
REAL(8):: newWeight
TYPE(particle), POINTER:: newPart
INTEGER:: p
INTEGER:: sp
newWeight = part%weight / nSplit
@ -779,12 +711,14 @@ MODULE moduleSolver
ALLOCATE(newPart)
!Copy data from original particle
newPart = part
!Add particle to list of new particles from weighting scheme
CALL OMP_SET_LOCK(lockWScheme)
CALL partWScheme%add(newPart)
CALL OMP_UNSET_LOCK(lockWScheme)
!Add particle to cell list
CALL OMP_SET_lock(vol%lock)
CALL vol%listPart_in%add(newPart)
sp = part%species%n
CALL vol%listPart_in(sp)%add(newPart)
CALL OMP_UNSET_lock(vol%lock)
END DO
@ -800,6 +734,9 @@ MODULE moduleSolver
TYPE(particle), INTENT(inout):: part
CLASS(meshVol), POINTER:: volOld, volNew
!Assume that particle is outside the domain
part%n_in = .FALSE.
volOld => mesh%vols(part%vol)%obj
CALL volOld%findCell(part)
CALL findCellColl(part)

View file

@ -6,7 +6,7 @@ MODULE moduleSpecies
TYPE, ABSTRACT:: speciesGeneric
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: m=0.D0, weight=0.D0
REAL(8):: m=0.D0, weight=0.D0, qm=0.D0
INTEGER:: n=0
END TYPE speciesGeneric
@ -18,7 +18,7 @@ MODULE moduleSpecies
END TYPE speciesNeutral
TYPE, EXTENDS(speciesGeneric):: speciesCharged
REAL(8):: q=0.D0, qm=0.D0
REAL(8):: q=0.D0
CLASS(speciesGeneric), POINTER:: ion => NULL(), neutral => NULL()
CONTAINS
PROCEDURE, PASS:: ionize => ionizeCharged
@ -43,7 +43,6 @@ MODULE moduleSpecies
REAL(8):: xi(1:3) !Logical coordinates of particle in element e_p.
LOGICAL:: n_in !Flag that indicates if a particle is in the domain
REAL(8):: weight=0.D0 !weight of particle
REAL(8):: qm = 0.D0 !charge over mass
END TYPE particle
@ -91,7 +90,7 @@ MODULE moduleSpecies
part%species => self%ion
ELSE
CALL criticalError('No ion defined for species' // self%name, 'ionizeNeutral')
CALL criticalError('No ion defined for species :' // self%name, 'ionizeNeutral')
END IF
@ -109,7 +108,7 @@ MODULE moduleSpecies
part%species => self%ion
ELSE
CALL criticalError('No ion defined for species' // self%name, 'ionizeCharged')
CALL criticalError('No ion defined for species :' // self%name, 'ionizeCharged')
END IF