Line | Branch | Exec | Source |
---|---|---|---|
1 | #include <fstream> | ||
2 | #include <iostream> | ||
3 | #include <limits> | ||
4 | #include <sstream> | ||
5 | |||
6 | #include "analyzers.hpp" | ||
7 | #include "interaction.hpp" | ||
8 | #include "multiindex.hpp" | ||
9 | #include "particle.hpp" | ||
10 | #include "rng.hpp" | ||
11 | #include "system.hpp" | ||
12 | |||
13 | 75 | System::System(std::vector<Particle> &particles, const YAML::Node &config, double t, | |
14 | 75 | const YAML::Node &addParticlesConfig) | |
15 |
1/2✓ Branch 2 taken 75 times.
✗ Branch 3 not taken.
|
150 | : particles(particles), t(t), L(YAML::parse<Vec>(config, "L")), |
16 |
1/2✓ Branch 2 taken 75 times.
✗ Branch 3 not taken.
|
150 | T(YAML::parse<double>(config, "T", std::numeric_limits<double>::quiet_NaN())), |
17 |
1/2✓ Branch 2 taken 75 times.
✗ Branch 3 not taken.
|
150 | mu(YAML::parse<double>(config, "mu", std::numeric_limits<double>::quiet_NaN())), E(0.), |
18 |
2/4✓ Branch 2 taken 75 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 75 times.
✗ Branch 7 not taken.
|
75 | interaction(InteractionFactory::create(this, config["interaction"])) { |
19 |
1/2✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
|
75 | std::clog << "\nSystem:"; |
20 |
1/2✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
|
75 | std::clog << "\n\tL:"; |
21 |
2/2✓ Branch 0 taken 225 times.
✓ Branch 1 taken 75 times.
|
300 | for (int i = 0; i < DIM; ++i) { |
22 |
2/4✓ Branch 1 taken 225 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 225 times.
✗ Branch 6 not taken.
|
225 | std::clog << " " << this->L[i]; |
23 | } | ||
24 |
1/2✓ Branch 0 taken 75 times.
✗ Branch 1 not taken.
|
75 | if (std::isfinite(T)) { |
25 |
2/4✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 75 times.
✗ Branch 5 not taken.
|
75 | std::clog << "\n\tT: " << this->T; |
26 | } | ||
27 |
2/2✓ Branch 0 taken 50 times.
✓ Branch 1 taken 25 times.
|
75 | if (std::isfinite(mu)) { |
28 |
2/4✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
|
50 | std::clog << "\n\tmu: " << this->mu; |
29 | } | ||
30 |
1/2✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
|
75 | std::clog << "\nParticles:"; |
31 |
3/6✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 75 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 75 times.
✗ Branch 8 not taken.
|
75 | std::clog << "\n\t" << particles.size() << " already in the system"; |
32 |
1/2✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
|
75 | initParticles(addParticlesConfig); |
33 |
5/10✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 75 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 75 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 75 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 75 times.
✗ Branch 13 not taken.
|
150 | std::clog << "\n\ttotal: " << particles.size() << " (N/V = " << particles.size() / L.prod() |
34 |
1/2✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
|
75 | << ")"; |
35 |
2/4✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 75 times.
|
75 | if (!allParticlesInside()) { |
36 | ✗ | throw std::runtime_error( | |
37 | ✗ | "Could not initialize the system because some particles would be outside of the box."); | |
38 | } | ||
39 |
2/2✓ Branch 0 taken 50 times.
✓ Branch 1 taken 25 times.
|
75 | if (particles.size() == 0) { |
40 | 50 | std::clog << "\nWARNING: Initialized system with zero particles. Ensure that this is intended, " | |
41 |
1/2✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
|
50 | "otherwise check the particle initialization in the config!"; |
42 | } | ||
43 | 75 | } | |
44 | |||
45 |
1/2✓ Branch 0 taken 73 times.
✗ Branch 1 not taken.
|
73 | System::~System() { delete interaction; } |
46 | |||
47 | ✗ | void System::ensureAnalyzer(const std::string &key) { | |
48 | ✗ | if (!analyzers[key]) { | |
49 | ✗ | analyzers[key] = AnalyzerFactory::create(this, key); | |
50 | } | ||
51 | } | ||
52 | |||
53 | ✗ | void System::updateAnalyzers() { | |
54 | ✗ | for (auto &[_, analyzer] : analyzers) { | |
55 | ✗ | analyzer->update(); | |
56 | } | ||
57 | } | ||
58 | |||
59 | 10 | int System::getNumSpecies() { | |
60 | 10 | int maxSpecies = 0; | |
61 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 10 times.
|
12 | for (auto &p : particles) { |
62 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | maxSpecies = std::max(maxSpecies, p.species); |
63 | } | ||
64 | 10 | return maxSpecies + 1; | |
65 | } | ||
66 | |||
67 | 4885 | bool System::removeLastParticleIfHighEnergy(double threshold) { | |
68 | 4885 | bool wasRemoved = false; | |
69 | 4885 | try { | |
70 | 4885 | double Eparticle = 0.; | |
71 |
1/2✓ Branch 1 taken 4885 times.
✗ Branch 2 not taken.
|
4885 | Eparticle = updateParticleState(particles.back(), ParticleUpdateType::INSERT); |
72 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4885 times.
|
4885 | if (Eparticle >= threshold) { |
73 | ✗ | std::clog << "\n\tParticle insertion skipped due to energy threshold (E=" << Eparticle | |
74 | ✗ | << " >= threshold=" << threshold << ")"; | |
75 | ✗ | updateParticleMeta(particles.back(), ParticleUpdateType::DELETE); | |
76 | ✗ | particles.pop_back(); | |
77 | ✗ | wasRemoved = true; | |
78 | } | ||
79 | ✗ | } catch (const std::runtime_error &e) { | |
80 | ✗ | std::clog << "\nWARNING: Could not check energy threshold of particle insertion because " | |
81 | ✗ | "interaction does not support single particle updates"; | |
82 | } | ||
83 | 4885 | return wasRemoved; | |
84 | } | ||
85 | |||
86 | 78 | void System::initParticles(const YAML::Node &addParticlesConfig) { | |
87 |
1/2✓ Branch 0 taken 78 times.
✗ Branch 1 not taken.
|
78 | if (addParticlesConfig) { |
88 |
2/2✓ Branch 1 taken 1 times.
✓ Branch 2 taken 77 times.
|
78 | if (addParticlesConfig.IsSequence()) { |
89 |
4/6✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✓ Branch 9 taken 1 times.
|
3 | for (auto &c : addParticlesConfig) { |
90 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | initParticlesItem(c); |
91 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
3 | } |
92 | } else { | ||
93 | 77 | initParticlesItem(addParticlesConfig); | |
94 | } | ||
95 | } | ||
96 | 78 | } | |
97 | |||
98 | 78 | void System::initParticlesItem(const YAML::Node &addParticlesConfig) { | |
99 | 78 | int nBefore = particles.size(); | |
100 | 78 | double defaultEThreshold = 1e4; | |
101 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 78 times.
|
78 | if (addParticlesConfig["file"]) { |
102 | ✗ | auto &config = addParticlesConfig["file"]; | |
103 | ✗ | std::string filename = YAML::parse<std::string>(config, "filename"); | |
104 | ✗ | std::ifstream particlesInitStream(filename, std::ios_base::in); | |
105 | ✗ | if (particlesInitStream.fail()) { | |
106 | ✗ | throw std::invalid_argument("Could not open initialization file " + filename); | |
107 | } | ||
108 | ✗ | double EThreshold = YAML::parse<double>(config, "EThreshold", defaultEThreshold); | |
109 | ✗ | initParticlesFromFile(particlesInitStream, EThreshold); | |
110 | ✗ | std::clog << "\n\tAdded " << particles.size() - nBefore << " from file " << filename; | |
111 | } | ||
112 |
2/2✓ Branch 2 taken 28 times.
✓ Branch 3 taken 50 times.
|
78 | if (addParticlesConfig["lattice"]) { |
113 | 28 | auto &config = addParticlesConfig["lattice"]; | |
114 |
3/6✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
|
56 | std::string cell = YAML::parse<std::string>(config, "cell", "sc"); |
115 |
3/6✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
|
56 | Vec spacing = YAML::parse<Vec>(config, "spacing", 1.); |
116 |
3/6✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
|
56 | Vec margin = YAML::parse<Vec>(config, "margin", 0.5); |
117 |
3/6✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
|
56 | Vec offset = YAML::parse<Vec>(config, "offset", 0.); |
118 |
3/6✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 28 times.
✗ Branch 9 not taken.
|
56 | int NTarget = YAML::parse<int>(config, "N", 0); |
119 |
3/6✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 28 times.
✗ Branch 9 not taken.
|
56 | double EThreshold = YAML::parse<double>(config, "EThreshold", defaultEThreshold); |
120 |
1/2✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
|
28 | YAML::Node particleProperties; |
121 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
|
28 | if (config["particleProperties"]) { |
122 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
|
28 | particleProperties = config["particleProperties"]; |
123 | } | ||
124 |
1/2✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
|
28 | initParticlesOnLattice(cell, spacing, margin, offset, particleProperties, NTarget, EThreshold); |
125 |
5/10✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 28 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 28 times.
✗ Branch 14 not taken.
|
28 | std::clog << "\n\tAdded " << particles.size() - nBefore << " on " << cell << " lattice"; |
126 | 56 | } | |
127 | 78 | } | |
128 | |||
129 | 16 | void System::initParticlesFromFile(std::istream &stream, double EThreshold) { | |
130 | // Not very robust, I know | ||
131 | 16 | std::string line; | |
132 |
3/4✓ Branch 1 taken 171 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 155 times.
✓ Branch 4 taken 16 times.
|
171 | while (std::getline(stream, line)) { |
133 |
4/4✓ Branch 0 taken 153 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 152 times.
|
155 | if (line.empty() || line.front() == '#') { |
134 | 3 | continue; | |
135 | } | ||
136 |
1/2✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
|
152 | std::istringstream iss(line); |
137 |
1/2✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
|
152 | iss.exceptions(std::ios::failbit); |
138 | 152 | Vec r; | |
139 | 152 | Vec v = Vec::Constant(0.); | |
140 | 152 | Vec u = Vec::Constant(0.); | |
141 | 152 | Real sigma = 1; | |
142 | 152 | Real m = 1; | |
143 | 152 | Real gamma = 1; | |
144 | 152 | int species = 0; | |
145 |
2/2✓ Branch 0 taken 456 times.
✓ Branch 1 taken 152 times.
|
608 | for (int i = 0; i < DIM; ++i) { |
146 | 456 | Real ri; | |
147 |
1/2✓ Branch 1 taken 456 times.
✗ Branch 2 not taken.
|
456 | iss >> ri; |
148 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 456 times.
|
456 | if (ri >= L[i]) { |
149 | ✗ | throw std::invalid_argument("Particle " + std::to_string(i) + | |
150 | ✗ | " is not inside simulation box. Aborting"); | |
151 | } | ||
152 | 456 | r[i] = ri; | |
153 | } | ||
154 | try { | ||
155 |
2/2✓ Branch 0 taken 456 times.
✓ Branch 1 taken 152 times.
|
608 | for (int i = 0; i < DIM; ++i) { |
156 |
1/2✓ Branch 2 taken 456 times.
✗ Branch 3 not taken.
|
456 | iss >> v[i]; |
157 | } | ||
158 |
1/2✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
|
152 | iss >> sigma; |
159 |
1/2✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
|
152 | iss >> m; |
160 |
1/2✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
|
152 | iss >> gamma; |
161 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 152 times.
|
152 | iss >> species; |
162 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 152 times.
|
152 | } catch (const std::istream::failure &e) { |
163 | // ignore optional properties | ||
164 | 152 | } | |
165 |
1/2✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
|
152 | particles.emplace_back(r, v, u, sigma, m, gamma, species); |
166 |
1/2✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
|
152 | removeLastParticleIfHighEnergy(EThreshold); |
167 | 152 | } | |
168 | 16 | } | |
169 | |||
170 | 28 | static void removeAverageVelocity(std::vector<Particle> &particles) { | |
171 | 28 | Vec vAvg = Vec::Constant(0.); | |
172 |
2/2✓ Branch 0 taken 4733 times.
✓ Branch 1 taken 28 times.
|
4761 | for (const Particle &p : particles) { |
173 | 4733 | vAvg += p.v; | |
174 | } | ||
175 | 28 | vAvg /= particles.size(); | |
176 |
2/2✓ Branch 0 taken 4733 times.
✓ Branch 1 taken 28 times.
|
4761 | for (Particle &p : particles) { |
177 | 4733 | p.v -= vAvg; | |
178 | } | ||
179 | 28 | } | |
180 | |||
181 | 28 | void System::initParticlesOnLattice(const std::string &cell, Vec spacing, const Vec &margin, | |
182 | const Vec &offset, const YAML::Node &particleProperties, | ||
183 | int NTarget, double EThreshold) { | ||
184 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 28 times.
|
28 | if ((spacing.array() <= 0.).any() || (spacing.array() == 1.).all()) { |
185 | ✗ | spacing = Vec::Constant(1. + 1e-10); // Safety margin due to floating point errors | |
186 | } | ||
187 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 27 times.
|
28 | Vec a = spacing; |
188 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 27 times.
|
28 | if (cell == "bcc") { |
189 | 1 | a *= 2. / std::sqrt(3.); | |
190 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 26 times.
|
27 | } else if (cell == "fcc") { |
191 | 1 | a *= std::sqrt(2.); | |
192 | } | ||
193 |
2/2✓ Branch 0 taken 26 times.
✓ Branch 1 taken 2 times.
|
28 | std::vector<Vec> bravais; |
194 | 28 | if constexpr (DIM == 2) { | |
195 | if (cell == "sc") { | ||
196 | bravais.push_back(Vec::Constant(0.)); | ||
197 | } else { | ||
198 | throw std::invalid_argument("Grid type " + cell + " not implemented!"); | ||
199 | } | ||
200 | } | ||
201 | 28 | if constexpr (DIM == 3) { | |
202 |
2/2✓ Branch 0 taken 26 times.
✓ Branch 1 taken 2 times.
|
28 | if (cell == "sc") { |
203 |
1/2✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
|
26 | bravais.push_back(Vec::Constant(0.)); |
204 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
|
2 | } else if (cell == "bcc") { |
205 |
1/2✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
1 | bravais.push_back(Vec::Constant(0.)); |
206 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | bravais.push_back(Vec{0.5, 0.5, 0.5}); |
207 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
1 | } else if (cell == "fcc") { |
208 |
1/2✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
1 | bravais.push_back(Vec::Constant(0.)); |
209 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | bravais.push_back(Vec{0., 0.5, 0.5}); |
210 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | bravais.push_back(Vec{0.5, 0., 0.5}); |
211 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | bravais.push_back(Vec{0.5, 0.5, 0.}); |
212 | } else { | ||
213 | ✗ | throw std::invalid_argument("Grid type " + cell + " not implemented!"); | |
214 | } | ||
215 | } | ||
216 | |||
217 |
3/6✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 28 times.
✗ Branch 9 not taken.
|
56 | Real sigma = YAML::parse<Real>(particleProperties, "sigma", 1.); |
218 |
3/6✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 28 times.
✗ Branch 9 not taken.
|
56 | Real m = YAML::parse<Real>(particleProperties, "m", 1.); |
219 |
3/6✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 28 times.
✗ Branch 9 not taken.
|
56 | Real gamma = YAML::parse<Real>(particleProperties, "gamma", 1.); |
220 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
|
56 | int species = YAML::parse<int>(particleProperties, "species", 0); |
221 | 28 | bool vMaxwell = false; | |
222 | 28 | bool uRandom = false; | |
223 | 28 | Vec v, u = Vec::Zero(); | |
224 | 28 | try { | |
225 |
3/6✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 28 times.
|
84 | if (YAML::parse<std::string>(particleProperties, "v", "") == "maxwell") { |
226 | ✗ | vMaxwell = true; | |
227 | } | ||
228 | ✗ | } catch (const YAML::BadConversion &e) { | |
229 | ✗ | v = YAML::parse<Vec>(particleProperties, "v", 0.); | |
230 | } | ||
231 | 28 | try { | |
232 |
3/6✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 28 times.
|
84 | if (YAML::parse<std::string>(particleProperties, "u", "") == "random") { |
233 | ✗ | uRandom = true; | |
234 | } | ||
235 | ✗ | } catch (const YAML::BadConversion &e) { | |
236 | ✗ | u = YAML::parse<Vec>(particleProperties, "u", 0.); | |
237 | } | ||
238 | |||
239 |
1/2✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
|
28 | RNG rng; |
240 | 28 | std::gamma_distribution<double> maxwellDist(1.5, T); | |
241 | |||
242 | 28 | size_t NSpecies = 0; | |
243 | 28 | ArrayI bounds = (L.array() / a.array() + 1).cast<int>(); | |
244 |
2/2✓ Branch 1 taken 10174 times.
✓ Branch 2 taken 28 times.
|
10202 | for (ArrayI multi : Multiindex<DIM>(bounds)) { |
245 | 10174 | Vec rc = multi.cast<Real>() * a.array(); | |
246 |
2/2✓ Branch 0 taken 10342 times.
✓ Branch 1 taken 10174 times.
|
20516 | for (auto &b : bravais) { |
247 | 10342 | Vec r = rc + (b.array() * a.array()).matrix() + offset; | |
248 |
2/2✓ Branch 0 taken 4733 times.
✓ Branch 1 taken 5609 times.
|
10342 | if ((r.array() >= margin.array() && r.array() < L.array() - margin.array()).all()) { |
249 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4733 times.
|
4733 | if (vMaxwell) { |
250 | ✗ | v = rng.vecNormal().normalized(); | |
251 | ✗ | v *= maxwellDist(rng.engine); | |
252 | } | ||
253 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4733 times.
|
4733 | if (uRandom) { |
254 | ✗ | do { | |
255 | ✗ | u = rng.vecNormal(); | |
256 | ✗ | } while (u.squaredNorm() - 1. > 1e-5); | |
257 | } | ||
258 |
1/2✓ Branch 1 taken 4733 times.
✗ Branch 2 not taken.
|
4733 | u.normalize(); |
259 |
1/2✓ Branch 1 taken 4733 times.
✗ Branch 2 not taken.
|
4733 | particles.emplace_back(r, v, u, sigma, m, gamma, species); |
260 |
2/4✓ Branch 1 taken 4733 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4733 times.
✗ Branch 4 not taken.
|
4733 | if (!removeLastParticleIfHighEnergy(EThreshold)) { |
261 | 4733 | ++NSpecies; | |
262 | } | ||
263 |
1/3✗ Branch 0 not taken.
✓ Branch 1 taken 4733 times.
✗ Branch 3 not taken.
|
4733 | if (NTarget && NSpecies == (size_t)NTarget) { |
264 | ✗ | removeAverageVelocity(particles); | |
265 | ✗ | return; | |
266 | } | ||
267 | } | ||
268 | } | ||
269 | } | ||
270 | 28 | removeAverageVelocity(particles); | |
271 |
1/3✗ Branch 0 not taken.
✓ Branch 1 taken 28 times.
✗ Branch 3 not taken.
|
28 | if (NTarget && NSpecies < (size_t)NTarget) { |
272 | ✗ | throw std::runtime_error("Could not fit " + std::to_string(NTarget) + " particles of species" + | |
273 | ✗ | std::to_string(species) + " into the box (only got to " + | |
274 | ✗ | std::to_string(NSpecies) + ")."); | |
275 | } | ||
276 | 28 | } | |
277 | |||
278 | 5 | void System::updateParticleMeta(Particle &p, ParticleUpdateType particleUpdateType) { | |
279 | 5 | interaction->updateParticleMeta(p, particleUpdateType); | |
280 | 5 | } | |
281 | |||
282 | 5094 | double System::updateParticleState(Particle &p, ParticleUpdateType particleUpdateType) { | |
283 | 5094 | p.Fext.setZero(); | |
284 | 5094 | p.Eext = 0.; | |
285 | 5094 | p.Fint.setZero(); | |
286 | 5094 | return interaction->updateParticleState(p, particleUpdateType); | |
287 | } | ||
288 | |||
289 | 4495 | void System::updateState() { | |
290 | 4495 | E = 0.; | |
291 |
2/2✓ Branch 0 taken 18057 times.
✓ Branch 1 taken 4495 times.
|
22552 | for (auto &p : particles) { |
292 | 18057 | p.Fext.setZero(); | |
293 | 18057 | p.Eext = 0.; | |
294 | 18057 | p.Fint.setZero(); | |
295 | } | ||
296 | 4495 | interaction->updateState(); | |
297 | 4495 | } | |
298 | |||
299 | 13873 | bool System::particleInside(const Particle &p) { | |
300 |
2/2✓ Branch 1 taken 3133 times.
✓ Branch 2 taken 10740 times.
|
13873 | if ((p.r.array() >= L.array() || p.r.array() < 0).any()) { |
301 | 3133 | return false; | |
302 | } | ||
303 | return true; | ||
304 | } | ||
305 | |||
306 | 78 | bool System::allParticlesInside() { | |
307 |
2/2✓ Branch 0 taken 4733 times.
✓ Branch 1 taken 78 times.
|
4811 | for (auto &p : particles) { |
308 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 4733 times.
|
4733 | if (!particleInside(p)) { |
309 | ✗ | return false; | |
310 | } | ||
311 | } | ||
312 | 78 | return true; | |
313 | } | ||
314 | |||
315 | 6007 | void System::applyPeriodic(Particle &p) { | |
316 |
2/2✓ Branch 1 taken 3133 times.
✓ Branch 2 taken 6007 times.
|
9140 | while (!particleInside(p)) { |
317 | 3133 | Vec shift = (p.r.array() / L.array()).floor(); | |
318 | 3133 | p.box += shift.cast<int>().array(); | |
319 | 6266 | p.r.array() -= shift.array() * L.array(); | |
320 | } | ||
321 | 6007 | } | |
322 | |||
323 | 24 | void System::applyPeriodic() { | |
324 | 24 | #pragma omp parallel for | |
325 | for (auto &p : particles) { | ||
326 | applyPeriodic(p); | ||
327 | } | ||
328 | 24 | } | |
329 |