AutoPas  3.0.0
Loading...
Searching...
No Matches
SPHCalcHydroForceFunctor.h
Go to the documentation of this file.
1
7#pragma once
8
9#include "SPHKernels.h"
10#include "autopas/particles/OwnershipState.h"
11
12namespace sphLib {
18template <class Particle_T>
19class SPHCalcHydroForceFunctor : public autopas::PairwiseFunctor<Particle_T, SPHCalcHydroForceFunctor<Particle_T>> {
20 public:
22 using SoAArraysType = typename Particle_T::SoAArraysType;
23
24 SPHCalcHydroForceFunctor()
25 // the actual cutoff used is dynamic. 0 is used to pass the sanity check.
26 : autopas::PairwiseFunctor<Particle_T, SPHCalcHydroForceFunctor<Particle_T>>(0.){};
27
28 virtual std::string getName() override { return "SPHHydroForceFunctor"; }
29
30 bool isRelevantForTuning() override { return true; }
31
32 bool allowsNewton3() override { return true; }
33
34 bool allowsNonNewton3() override { return true; }
35
45 void AoSFunctor(Particle_T &i, Particle_T &j, bool newton3 = true) override {
46 using namespace autopas::utils::ArrayMath::literals;
47
48 if (i.isDummy() or j.isDummy()) {
49 return;
50 }
51 const std::array<double, 3> dr = i.getR() - j.getR();
52 // const PS::F64vec dr = ep_i[i].pos - ep_j[j].pos;
53
54 double cutoff = i.getSmoothingLength() * sphLib::SPHKernels::getKernelSupportRadius();
55
56 if (autopas::utils::ArrayMath::dot(dr, dr) >= cutoff * cutoff) {
57 return;
58 }
59
60 const std::array<double, 3> dv = i.getV() - j.getV();
61 // const PS::F64vec dv = ep_i[i].vel - ep_j[j].vel;
62
63 double dvdr = autopas::utils::ArrayMath::dot(dv, dr);
64 const double w_ij = (dvdr < 0) ? dvdr / autopas::utils::ArrayMath::L2Norm(dr) : 0;
65 // const PS::F64 w_ij = (dv * dr < 0) ? dv * dr / sqrt(dr * dr) : 0;
66
67 const double v_sig = i.getSoundSpeed() + j.getSoundSpeed() - 3.0 * w_ij;
68 // const PS::F64 v_sig = ep_i[i].snds + ep_j[j].snds - 3.0 * w_ij;
69
70 i.checkAndSetVSigMax(v_sig);
71 if (newton3) {
72 j.checkAndSetVSigMax(v_sig); // Newton 3
73 // v_sig_max = std::max(v_sig_max, v_sig);
74 }
75 const double AV = -0.5 * v_sig * w_ij / (0.5 * (i.getDensity() + j.getDensity()));
76 // const PS::F64 AV = - 0.5 * v_sig * w_ij / (0.5 * (ep_i[i].dens +
77 // ep_j[j].dens));
78
79 const std::array<double, 3> gradW_ij =
80 (SPHKernels::gradW(dr, i.getSmoothingLength()) + SPHKernels::gradW(dr, j.getSmoothingLength())) * 0.5;
81 // const PS::F64vec gradW_ij = 0.5 * (gradW(dr, ep_i[i].smth) + gradW(dr,
82 // ep_j[j].smth));
83
84 double scale =
85 i.getPressure() / (i.getDensity() * i.getDensity()) + j.getPressure() / (j.getDensity() * j.getDensity()) + AV;
86 i.subAcceleration(gradW_ij * (scale * j.getMass()));
87 // hydro[i].acc -= ep_j[j].mass * (ep_i[i].pres / (ep_i[i].dens *
88 // ep_i[i].dens) + ep_j[j].pres / (ep_j[j].dens * ep_j[j].dens) + AV) *
89 // gradW_ij;
90 if (newton3) {
91 j.addAcceleration(gradW_ij * (scale * i.getMass()));
92 // Newton3, gradW_ij = -gradW_ji
93 }
94 double scale2i = j.getMass() * (i.getPressure() / (i.getDensity() * i.getDensity()) + 0.5 * AV);
95 i.addEngDot(autopas::utils::ArrayMath::dot(gradW_ij, dv) * scale2i);
96 // hydro[i].eng_dot += ep_j[j].mass * (ep_i[i].pres / (ep_i[i].dens *
97 // ep_i[i].dens) + 0.5 * AV) * dv * gradW_ij;
98
99 if (newton3) {
100 double scale2j = i.getMass() * (j.getPressure() / (j.getDensity() * j.getDensity()) + 0.5 * AV);
101 j.addEngDot(autopas::utils::ArrayMath::dot(gradW_ij, dv) * scale2j);
102 // Newton 3
103 }
104 }
105
110 void SoAFunctorSingle(autopas::SoAView<SoAArraysType> soa, bool newton3) override {
111 using namespace autopas::utils::ArrayMath::literals;
112 if (soa.size() == 0) return;
113
114 double *const __restrict massptr = soa.template begin<Particle_T::AttributeNames::mass>();
115 double *const __restrict densityptr = soa.template begin<Particle_T::AttributeNames::density>();
116 double *const __restrict smthptr = soa.template begin<Particle_T::AttributeNames::smth>();
117 double *const __restrict soundSpeedptr = soa.template begin<Particle_T::AttributeNames::soundSpeed>();
118 double *const __restrict pressureptr = soa.template begin<Particle_T::AttributeNames::pressure>();
119 double *const __restrict vsigmaxptr = soa.template begin<Particle_T::AttributeNames::vsigmax>();
120 double *const __restrict engDotptr = soa.template begin<Particle_T::AttributeNames::engDot>();
121
122 double *const __restrict xptr = soa.template begin<Particle_T::AttributeNames::posX>();
123 double *const __restrict yptr = soa.template begin<Particle_T::AttributeNames::posY>();
124 double *const __restrict zptr = soa.template begin<Particle_T::AttributeNames::posZ>();
125 double *const __restrict velXptr = soa.template begin<Particle_T::AttributeNames::velX>();
126 double *const __restrict velYptr = soa.template begin<Particle_T::AttributeNames::velY>();
127 double *const __restrict velZptr = soa.template begin<Particle_T::AttributeNames::velZ>();
128 double *const __restrict accXptr = soa.template begin<Particle_T::AttributeNames::accX>();
129 double *const __restrict accYptr = soa.template begin<Particle_T::AttributeNames::accY>();
130 double *const __restrict accZptr = soa.template begin<Particle_T::AttributeNames::accZ>();
131
132 const auto *const __restrict ownedStatePtr = soa.template begin<Particle_T::AttributeNames::ownershipState>();
133
134 for (unsigned int indexFirst = 0; indexFirst < soa.size(); ++indexFirst) {
135 // checks whether particle i is owned.
136 if (ownedStatePtr[indexFirst] == autopas::OwnershipState::dummy) {
137 continue;
138 }
139
140 double localvsigmax = 0.;
141 double localengdotsum = 0.;
142 double localAccX = 0.;
143 double localAccY = 0.;
144 double localAccZ = 0.;
145
146 // icpc vectorizes this.
147 // g++ only with -ffast-math or -funsafe-math-optimizations
148 // #pragma omp simd reduction(+ : localengdotsum, localAccX, localAccY, localAccZ), reduction(max : localvsigmax)
149 for (unsigned int j = indexFirst + 1; j < soa.size(); ++j) {
150 using namespace autopas::utils::ArrayMath::literals;
151
152 const double drx = xptr[indexFirst] - xptr[j];
153 const double dry = yptr[indexFirst] - yptr[j];
154 const double drz = zptr[indexFirst] - zptr[j];
155
156 const double drx2 = drx * drx;
157 const double dry2 = dry * dry;
158 const double drz2 = drz * drz;
159
160 const double dr2 = drx2 + dry2 + drz2;
161 double cutoff = smthptr[indexFirst] * sphLib::SPHKernels::getKernelSupportRadius();
162 if (dr2 >= cutoff * cutoff or ownedStatePtr[j] == autopas::OwnershipState::dummy) continue;
163
164 const double dvX = velXptr[indexFirst] - velXptr[j];
165 const double dvY = velYptr[indexFirst] - velYptr[j];
166 const double dvZ = velZptr[indexFirst] - velZptr[j];
167 // const PS::F64vec dv = ep_i[i].vel - ep_j[j].vel;
168
169 double dvdr = dvX * drx + dvY * dry + dvZ * drz;
170 const double w_ij = (dvdr < 0) ? dvdr / sqrt(dr2) : 0;
171 // const PS::F64 w_ij = (dv * dr < 0) ? dv * dr / sqrt(dr * dr) : 0;
172
173 const double v_sig = soundSpeedptr[indexFirst] + soundSpeedptr[j] - 3.0 * w_ij;
174 // const PS::F64 v_sig = ep_i[i].snds + ep_j[j].snds - 3.0 * w_ij;
175
176 localvsigmax = std::max(localvsigmax, v_sig);
177 // vsigmaxptr[j] = std::max(vsigmaxptr[j], v_sig); // Newton 3
178 vsigmaxptr[j] = vsigmaxptr[j] > v_sig ? vsigmaxptr[j] : v_sig; // Newton 3
179 // v_sig_max = std::max(v_sig_max, v_sig);
180
181 const double AV = -0.5 * v_sig * w_ij / (0.5 * (densityptr[indexFirst] + densityptr[j]));
182 // const PS::F64 AV = - 0.5 * v_sig * w_ij / (0.5 * (ep_i[i].dens +
183 // ep_j[j].dens));
184
185 const std::array<double, 3> gradW_ij =
186 (SPHKernels::gradW({drx, dry, drz}, smthptr[indexFirst]) + SPHKernels::gradW({drx, dry, drz}, smthptr[j])) *
187 0.5;
188 // const PS::F64vec gradW_ij = 0.5 * (gradW(dr, ep_i[i].smth) + gradW(dr,
189 // ep_j[j].smth));
190
191 double scale = pressureptr[indexFirst] / (densityptr[indexFirst] * densityptr[indexFirst]) +
192 pressureptr[j] / (densityptr[j] * densityptr[j]) + AV;
193 const double massscale = scale * massptr[j];
194 localAccX -= gradW_ij[0] * massscale;
195 localAccY -= gradW_ij[1] * massscale;
196 localAccZ -= gradW_ij[2] * massscale;
197 // hydro[i].acc -= ep_j[j].mass * (ep_i[i].pres / (ep_i[i].dens *
198 // ep_i[i].dens) + ep_j[j].pres / (ep_j[j].dens * ep_j[j].dens) + AV) *
199 // gradW_ij;
200
201 const double massscale2 = scale * massptr[indexFirst];
202 accXptr[j] += gradW_ij[0] * massscale2;
203 accYptr[j] += gradW_ij[1] * massscale2;
204 accZptr[j] += gradW_ij[2] * massscale2;
205 // Newton3, gradW_ij = -gradW_ji
206
207 double scale2i =
208 massptr[j] * (pressureptr[indexFirst] / (densityptr[indexFirst] * densityptr[indexFirst]) + 0.5 * AV);
209 localengdotsum += (gradW_ij[0] * dvX + gradW_ij[1] * dvY + gradW_ij[2] * dvZ) * scale2i;
210 // hydro[i].eng_dot += ep_j[j].mass * (ep_i[i].pres / (ep_i[i].dens *
211 // ep_i[i].dens) + 0.5 * AV) * dv * gradW_ij;
212
213 double scale2j = massptr[indexFirst] * (pressureptr[j] / (densityptr[j] * densityptr[j]) + 0.5 * AV);
214 engDotptr[j] += (gradW_ij[0] * dvX + gradW_ij[1] * dvY + gradW_ij[2] * dvZ) * scale2j;
215 // Newton 3
216 }
217
218 engDotptr[indexFirst] += localengdotsum;
219 accXptr[indexFirst] += localAccX;
220 accYptr[indexFirst] += localAccY;
221 accZptr[indexFirst] += localAccZ;
222 vsigmaxptr[indexFirst] = std::max(localvsigmax, vsigmaxptr[indexFirst]);
223 }
224 }
225
229 void SoAFunctorPair(autopas::SoAView<SoAArraysType> soa1, autopas::SoAView<SoAArraysType> soa2,
230 bool newton3) override {
231 using namespace autopas::utils::ArrayMath::literals;
232 if (soa1.size() == 0 || soa2.size() == 0) return;
233
234 double *const __restrict massptr1 = soa1.template begin<Particle_T::AttributeNames::mass>();
235 double *const __restrict densityptr1 = soa1.template begin<Particle_T::AttributeNames::density>();
236 double *const __restrict smthptr1 = soa1.template begin<Particle_T::AttributeNames::smth>();
237 double *const __restrict soundSpeedptr1 = soa1.template begin<Particle_T::AttributeNames::soundSpeed>();
238 double *const __restrict pressureptr1 = soa1.template begin<Particle_T::AttributeNames::pressure>();
239 double *const __restrict vsigmaxptr1 = soa1.template begin<Particle_T::AttributeNames::vsigmax>();
240 double *const __restrict engDotptr1 = soa1.template begin<Particle_T::AttributeNames::engDot>();
241
242 double *const __restrict xptr1 = soa1.template begin<Particle_T::AttributeNames::posX>();
243 double *const __restrict yptr1 = soa1.template begin<Particle_T::AttributeNames::posY>();
244 double *const __restrict zptr1 = soa1.template begin<Particle_T::AttributeNames::posZ>();
245 double *const __restrict velXptr1 = soa1.template begin<Particle_T::AttributeNames::velX>();
246 double *const __restrict velYptr1 = soa1.template begin<Particle_T::AttributeNames::velY>();
247 double *const __restrict velZptr1 = soa1.template begin<Particle_T::AttributeNames::velZ>();
248 double *const __restrict accXptr1 = soa1.template begin<Particle_T::AttributeNames::accX>();
249 double *const __restrict accYptr1 = soa1.template begin<Particle_T::AttributeNames::accY>();
250 double *const __restrict accZptr1 = soa1.template begin<Particle_T::AttributeNames::accZ>();
251
252 double *const __restrict massptr2 = soa2.template begin<Particle_T::AttributeNames::mass>();
253 double *const __restrict densityptr2 = soa2.template begin<Particle_T::AttributeNames::density>();
254 double *const __restrict smthptr2 = soa2.template begin<Particle_T::AttributeNames::smth>();
255 double *const __restrict soundSpeedptr2 = soa2.template begin<Particle_T::AttributeNames::soundSpeed>();
256 double *const __restrict pressureptr2 = soa2.template begin<Particle_T::AttributeNames::pressure>();
257 double *const __restrict vsigmaxptr2 = soa2.template begin<Particle_T::AttributeNames::vsigmax>();
258 double *const __restrict engDotptr2 = soa2.template begin<Particle_T::AttributeNames::engDot>();
259
260 double *const __restrict xptr2 = soa2.template begin<Particle_T::AttributeNames::posX>();
261 double *const __restrict yptr2 = soa2.template begin<Particle_T::AttributeNames::posY>();
262 double *const __restrict zptr2 = soa2.template begin<Particle_T::AttributeNames::posZ>();
263 double *const __restrict velXptr2 = soa2.template begin<Particle_T::AttributeNames::velX>();
264 double *const __restrict velYptr2 = soa2.template begin<Particle_T::AttributeNames::velY>();
265 double *const __restrict velZptr2 = soa2.template begin<Particle_T::AttributeNames::velZ>();
266 double *const __restrict accXptr2 = soa2.template begin<Particle_T::AttributeNames::accX>();
267 double *const __restrict accYptr2 = soa2.template begin<Particle_T::AttributeNames::accY>();
268 double *const __restrict accZptr2 = soa2.template begin<Particle_T::AttributeNames::accZ>();
269
270 const auto *const __restrict ownedStatePtr1 = soa1.template begin<Particle_T::AttributeNames::ownershipState>();
271 const auto *const __restrict ownedStatePtr2 = soa2.template begin<Particle_T::AttributeNames::ownershipState>();
272
273 for (unsigned int indexFirst = 0; indexFirst < soa1.size(); ++indexFirst) {
274 // checks whether particle i is owned.
275 if (ownedStatePtr1[indexFirst] == autopas::OwnershipState::dummy) {
276 continue;
277 }
278
279 double localvsigmax = 0.;
280 double localengdotsum = 0.;
281 double localAccX = 0.;
282 double localAccY = 0.;
283 double localAccZ = 0.;
284
285 // icpc vectorizes this.
286 // g++ only with -ffast-math or -funsafe-math-optimizations
287 // #pragma omp simd reduction(+ : localengdotsum, localAccX, localAccY, localAccZ), reduction(max : localvsigmax)
288 for (unsigned int j = 0; j < soa2.size(); ++j) {
289 using namespace autopas::utils::ArrayMath::literals;
290
291 const double drx = xptr1[indexFirst] - xptr2[j];
292 const double dry = yptr1[indexFirst] - yptr2[j];
293 const double drz = zptr1[indexFirst] - zptr2[j];
294
295 const double drx2 = drx * drx;
296 const double dry2 = dry * dry;
297 const double drz2 = drz * drz;
298
299 const double dr2 = drx2 + dry2 + drz2;
300 double cutoff = smthptr1[indexFirst] * sphLib::SPHKernels::getKernelSupportRadius();
301 if (dr2 >= cutoff * cutoff or ownedStatePtr2[j] == autopas::OwnershipState::dummy) continue;
302
303 const double dvX = velXptr1[indexFirst] - velXptr2[j];
304 const double dvY = velYptr1[indexFirst] - velYptr2[j];
305 const double dvZ = velZptr1[indexFirst] - velZptr2[j];
306 // const PS::F64vec dv = ep_i[i].vel - ep_j[j].vel;
307
308 double dvdr = dvX * drx + dvY * dry + dvZ * drz;
309 const double w_ij = (dvdr < 0) ? dvdr / sqrt(dr2) : 0;
310 // const PS::F64 w_ij = (dv * dr < 0) ? dv * dr / sqrt(dr * dr) : 0;
311
312 const double v_sig = soundSpeedptr1[indexFirst] + soundSpeedptr2[j] - 3.0 * w_ij;
313 // const PS::F64 v_sig = ep_i[i].snds + ep_j[j].snds - 3.0 * w_ij;
314
315 localvsigmax = std::max(localvsigmax, v_sig);
316 if (newton3) {
317 // vsigmaxptr2[j] = std::max(vsigmaxptr2[j], v_sig); // Newton 3
318 vsigmaxptr2[j] = vsigmaxptr2[j] > v_sig ? vsigmaxptr2[j] : v_sig; // Newton 3
319 // v_sig_max = std::max(v_sig_max, v_sig);
320 }
321 const double AV = -0.5 * v_sig * w_ij / (0.5 * (densityptr1[indexFirst] + densityptr2[j]));
322 // const PS::F64 AV = - 0.5 * v_sig * w_ij / (0.5 * (ep_i[i].dens +
323 // ep_j[j].dens));
324
325 const std::array<double, 3> gradW_ij = (SPHKernels::gradW({drx, dry, drz}, smthptr1[indexFirst]) +
326 SPHKernels::gradW({drx, dry, drz}, smthptr2[j])) *
327 0.5;
328 // const PS::F64vec gradW_ij = 0.5 * (gradW(dr, ep_i[i].smth) + gradW(dr,
329 // ep_j[j].smth));
330
331 double scale = pressureptr1[indexFirst] / (densityptr1[indexFirst] * densityptr1[indexFirst]) +
332 pressureptr2[j] / (densityptr2[j] * densityptr2[j]) + AV;
333 const double massscale = scale * massptr2[j];
334 localAccX -= gradW_ij[0] * massscale;
335 localAccY -= gradW_ij[1] * massscale;
336 localAccZ -= gradW_ij[2] * massscale;
337 // hydro[i].acc -= ep_j[j].mass * (ep_i[i].pres / (ep_i[i].dens *
338 // ep_i[i].dens) + ep_j[j].pres / (ep_j[j].dens * ep_j[j].dens) + AV) *
339 // gradW_ij;
340 if (newton3) {
341 const double massscale = scale * massptr1[indexFirst];
342 accXptr2[j] += gradW_ij[0] * massscale;
343 accYptr2[j] += gradW_ij[1] * massscale;
344 accZptr2[j] += gradW_ij[2] * massscale;
345 // Newton3, gradW_ij = -gradW_ji
346 }
347 double scale2i =
348 massptr2[j] * (pressureptr1[indexFirst] / (densityptr1[indexFirst] * densityptr1[indexFirst]) + 0.5 * AV);
349 localengdotsum += (gradW_ij[0] * dvX + gradW_ij[1] * dvY + gradW_ij[2] * dvZ) * scale2i;
350 // hydro[i].eng_dot += ep_j[j].mass * (ep_i[i].pres / (ep_i[i].dens *
351 // ep_i[i].dens) + 0.5 * AV) * dv * gradW_ij;
352
353 if (newton3) {
354 double scale2j = massptr1[indexFirst] * (pressureptr2[j] / (densityptr2[j] * densityptr2[j]) + 0.5 * AV);
355 engDotptr2[j] += (gradW_ij[0] * dvX + gradW_ij[1] * dvY + gradW_ij[2] * dvZ) * scale2j;
356 // Newton 3
357 }
358 }
359
360 engDotptr1[indexFirst] += localengdotsum;
361 accXptr1[indexFirst] += localAccX;
362 accYptr1[indexFirst] += localAccY;
363 accZptr1[indexFirst] += localAccZ;
364 vsigmaxptr1[indexFirst] = std::max(localvsigmax, vsigmaxptr1[indexFirst]);
365 }
366 }
367 // clang-format off
371 // clang-format on
372 void SoAFunctorVerlet(autopas::SoAView<SoAArraysType> soa, const size_t indexFirst,
373 const std::vector<size_t, autopas::AlignedAllocator<size_t>> &neighborList,
374 bool newton3) override {
375 using namespace autopas::utils::ArrayMath::literals;
376 if (soa.size() == 0) return;
377
378 const auto *const __restrict ownedStatePtr = soa.template begin<Particle_T::AttributeNames::ownershipState>();
379
380 // checks whether particle i is owned.
381 if (ownedStatePtr[indexFirst] == autopas::OwnershipState::dummy) {
382 return;
383 }
384
385 double *const __restrict massptr = soa.template begin<Particle_T::AttributeNames::mass>();
386 double *const __restrict densityptr = soa.template begin<Particle_T::AttributeNames::density>();
387 double *const __restrict smthptr = soa.template begin<Particle_T::AttributeNames::smth>();
388 double *const __restrict soundSpeedptr = soa.template begin<Particle_T::AttributeNames::soundSpeed>();
389 double *const __restrict pressureptr = soa.template begin<Particle_T::AttributeNames::pressure>();
390 double *const __restrict vsigmaxptr = soa.template begin<Particle_T::AttributeNames::vsigmax>();
391 double *const __restrict engDotptr = soa.template begin<Particle_T::AttributeNames::engDot>();
392
393 double *const __restrict xptr = soa.template begin<Particle_T::AttributeNames::posX>();
394 double *const __restrict yptr = soa.template begin<Particle_T::AttributeNames::posY>();
395 double *const __restrict zptr = soa.template begin<Particle_T::AttributeNames::posZ>();
396 double *const __restrict velXptr = soa.template begin<Particle_T::AttributeNames::velX>();
397 double *const __restrict velYptr = soa.template begin<Particle_T::AttributeNames::velY>();
398 double *const __restrict velZptr = soa.template begin<Particle_T::AttributeNames::velZ>();
399 double *const __restrict accXptr = soa.template begin<Particle_T::AttributeNames::accX>();
400 double *const __restrict accYptr = soa.template begin<Particle_T::AttributeNames::accY>();
401 double *const __restrict accZptr = soa.template begin<Particle_T::AttributeNames::accZ>();
402
403 double localvsigmax = 0.;
404 double localengdotsum = 0.;
405 double localAccX = 0.;
406 double localAccY = 0.;
407 double localAccZ = 0.;
408
409 const auto &currentList = neighborList;
410 size_t listSize = currentList.size();
411
412 // icpc vectorizes this.
413 // g++ only with -ffast-math or -funsafe-math-optimizations
414 // #pragma omp simd reduction(+ : localengdotsum, localAccX, localAccY, localAccZ), reduction(max : localvsigmax)
415 for (unsigned int j = 0; j < listSize; ++j) {
416 using namespace autopas::utils::ArrayMath::literals;
417
418 const double drx = xptr[indexFirst] - xptr[currentList[j]];
419 const double dry = yptr[indexFirst] - yptr[currentList[j]];
420 const double drz = zptr[indexFirst] - zptr[currentList[j]];
421
422 const double drx2 = drx * drx;
423 const double dry2 = dry * dry;
424 const double drz2 = drz * drz;
425
426 const double dr2 = drx2 + dry2 + drz2;
427 double cutoff = smthptr[indexFirst] * sphLib::SPHKernels::getKernelSupportRadius();
428 if (dr2 >= cutoff * cutoff or ownedStatePtr[currentList[j]] == autopas::OwnershipState::dummy) continue;
429
430 const double dvX = velXptr[indexFirst] - velXptr[currentList[j]];
431 const double dvY = velYptr[indexFirst] - velYptr[currentList[j]];
432 const double dvZ = velZptr[indexFirst] - velZptr[currentList[j]];
433 // const PS::F64vec dv = ep_i[i].vel - ep_j[currentList[j]].vel;
434
435 double dvdr = dvX * drx + dvY * dry + dvZ * drz;
436 const double w_ij = (dvdr < 0) ? dvdr / sqrt(dr2) : 0;
437 // const PS::F64 w_ij = (dv * dr < 0) ? dv * dr / sqrt(dr * dr) : 0;
438
439 const double v_sig = soundSpeedptr[indexFirst] + soundSpeedptr[currentList[j]] - 3.0 * w_ij;
440 // const PS::F64 v_sig = ep_i[i].snds + ep_j[currentList[j]].snds - 3.0 * w_ij;
441
442 localvsigmax = std::max(localvsigmax, v_sig);
443 if (newton3) {
444 // vsigmaxptr[currentList[j]] = std::max(vsigmaxptr[currentList[j]], v_sig); // Newton 3
445 vsigmaxptr[currentList[j]] =
446 vsigmaxptr[currentList[j]] > v_sig ? vsigmaxptr[currentList[j]] : v_sig; // Newton 3
447 // v_sig_max = std::max(v_sig_max, v_sig);
448 }
449 const double AV = -0.5 * v_sig * w_ij / (0.5 * (densityptr[indexFirst] + densityptr[currentList[j]]));
450 // const PS::F64 AV = - 0.5 * v_sig * w_ij / (0.5 * (ep_i[i].dens +
451 // ep_j[currentList[j]].dens));
452
453 const std::array<double, 3> gradW_ij = (SPHKernels::gradW({drx, dry, drz}, smthptr[indexFirst]) +
454 SPHKernels::gradW({drx, dry, drz}, smthptr[currentList[j]])) *
455 0.5;
456 // const PS::F64vec gradW_ij = 0.5 * (gradW(dr, ep_i[i].smth) + gradW(dr,
457 // ep_j[currentList[j]].smth));
458
459 double scale = pressureptr[indexFirst] / (densityptr[indexFirst] * densityptr[indexFirst]) +
460 pressureptr[currentList[j]] / (densityptr[currentList[j]] * densityptr[currentList[j]]) + AV;
461 const double massscale = scale * massptr[currentList[j]];
462 localAccX -= gradW_ij[0] * massscale;
463 localAccY -= gradW_ij[1] * massscale;
464 localAccZ -= gradW_ij[2] * massscale;
465 // hydro[i].acc -= ep_j[currentList[j]].mass * (ep_i[i].pres / (ep_i[i].dens *
466 // ep_i[i].dens) + ep_j[currentList[j]].pres / (ep_j[currentList[j]].dens * ep_j[currentList[j]].dens) + AV) *
467 // gradW_ij;
468 if (newton3) {
469 const double massscale = scale * massptr[indexFirst];
470 accXptr[currentList[j]] += gradW_ij[0] * massscale;
471 accYptr[currentList[j]] += gradW_ij[1] * massscale;
472 accZptr[currentList[j]] += gradW_ij[2] * massscale;
473 // Newton3, gradW_ij = -gradW_ji
474 }
475 double scale2i = massptr[currentList[j]] *
476 (pressureptr[indexFirst] / (densityptr[indexFirst] * densityptr[indexFirst]) + 0.5 * AV);
477 localengdotsum += (gradW_ij[0] * dvX + gradW_ij[1] * dvY + gradW_ij[2] * dvZ) * scale2i;
478 // hydro[i].eng_dot += ep_j[currentList[j]].mass * (ep_i[i].pres / (ep_i[i].dens *
479 // ep_i[i].dens) + 0.5 * AV) * dv * gradW_ij;
480
481 if (newton3) {
482 double scale2j =
483 massptr[indexFirst] *
484 (pressureptr[currentList[j]] / (densityptr[currentList[j]] * densityptr[currentList[j]]) + 0.5 * AV);
485 engDotptr[currentList[j]] += (gradW_ij[0] * dvX + gradW_ij[1] * dvY + gradW_ij[2] * dvZ) * scale2j;
486 // Newton 3
487 }
488 }
489
490 engDotptr[indexFirst] += localengdotsum;
491 accXptr[indexFirst] += localAccX;
492 accYptr[indexFirst] += localAccY;
493 accZptr[indexFirst] += localAccZ;
494 vsigmaxptr[indexFirst] = std::max(localvsigmax, vsigmaxptr[indexFirst]);
495 }
496
500 constexpr static auto getNeededAttr() {
501 return std::array<typename Particle_T::AttributeNames, 17>{
502 Particle_T::AttributeNames::mass, Particle_T::AttributeNames::density,
503 Particle_T::AttributeNames::smth, Particle_T::AttributeNames::soundSpeed,
504 Particle_T::AttributeNames::pressure, Particle_T::AttributeNames::vsigmax,
505 Particle_T::AttributeNames::engDot, Particle_T::AttributeNames::posX,
506 Particle_T::AttributeNames::posY, Particle_T::AttributeNames::posZ,
507 Particle_T::AttributeNames::velX, Particle_T::AttributeNames::velY,
508 Particle_T::AttributeNames::velZ, Particle_T::AttributeNames::accX,
509 Particle_T::AttributeNames::accY, Particle_T::AttributeNames::accZ,
510 Particle_T::AttributeNames::ownershipState};
511 }
512
516 constexpr static auto getNeededAttr(std::false_type) {
517 return std::array<typename Particle_T::AttributeNames, 12>{
518 Particle_T::AttributeNames::mass, Particle_T::AttributeNames::density,
519 Particle_T::AttributeNames::smth, Particle_T::AttributeNames::soundSpeed,
520 Particle_T::AttributeNames::pressure, Particle_T::AttributeNames::posX,
521 Particle_T::AttributeNames::posY, Particle_T::AttributeNames::posZ,
522 Particle_T::AttributeNames::velX, Particle_T::AttributeNames::velY,
523 Particle_T::AttributeNames::velZ, Particle_T::AttributeNames::ownershipState};
524 }
525
529 constexpr static auto getComputedAttr() {
530 return std::array<typename Particle_T::AttributeNames, 6>{
531 Particle_T::AttributeNames::vsigmax, Particle_T::AttributeNames::engDot,
532 Particle_T::AttributeNames::accX, Particle_T::AttributeNames::accY,
533 Particle_T::AttributeNames::accZ, Particle_T::AttributeNames::ownershipState};
534 }
535};
536
537} // namespace sphLib
autopas::AlignedAllocator
AlignedAllocator class.
Definition: AlignedAllocator.h:29
autopas::PairwiseFunctor
PairwiseFunctor class.
Definition: PairwiseFunctor.h:31
autopas::PairwiseFunctor< Particle_T, SPHCalcHydroForceFunctor< Particle_T > >::PairwiseFunctor
PairwiseFunctor(double cutoff)
Constructor.
Definition: PairwiseFunctor.h:42
autopas::SoAView
View on a fixed part of a SoA between a start index and an end index.
Definition: SoAView.h:23
autopas::SoAView::size
size_t size() const
Returns the number of particles in the view.
Definition: SoAView.h:83
sphLib::SPHCalcHydroForceFunctor
Class that defines the hydrodynamic force functor.
Definition: SPHCalcHydroForceFunctor.h:19
sphLib::SPHCalcHydroForceFunctor::SoAArraysType
typename Particle_T::SoAArraysType SoAArraysType
soa arrays type
Definition: SPHCalcHydroForceFunctor.h:22
sphLib::SPHCalcHydroForceFunctor::allowsNewton3
bool allowsNewton3() override
Specifies whether the functor is capable of Newton3-like functors.
Definition: SPHCalcHydroForceFunctor.h:32
sphLib::SPHCalcHydroForceFunctor::SoAFunctorSingle
void SoAFunctorSingle(autopas::SoAView< SoAArraysType > soa, bool newton3) override
PairwiseFunctor for structure of arrays (SoA)
Definition: SPHCalcHydroForceFunctor.h:110
sphLib::SPHCalcHydroForceFunctor::getNeededAttr
static constexpr auto getNeededAttr(std::false_type)
Get attributes needed for computation without N3 optimization.
Definition: SPHCalcHydroForceFunctor.h:516
sphLib::SPHCalcHydroForceFunctor::AoSFunctor
void AoSFunctor(Particle_T &i, Particle_T &j, bool newton3=true) override
Calculates the contribution of the interaction of particle i and j to the hydrodynamic force.
Definition: SPHCalcHydroForceFunctor.h:45
sphLib::SPHCalcHydroForceFunctor::getName
virtual std::string getName() override
Returns name of functor.
Definition: SPHCalcHydroForceFunctor.h:28
sphLib::SPHCalcHydroForceFunctor::getNeededAttr
static constexpr auto getNeededAttr()
Get attributes needed for computation.
Definition: SPHCalcHydroForceFunctor.h:500
sphLib::SPHCalcHydroForceFunctor::SoAFunctorPair
void SoAFunctorPair(autopas::SoAView< SoAArraysType > soa1, autopas::SoAView< SoAArraysType > soa2, bool newton3) override
PairwiseFunctor for structure of arrays (SoA)
Definition: SPHCalcHydroForceFunctor.h:229
sphLib::SPHCalcHydroForceFunctor::getComputedAttr
static constexpr auto getComputedAttr()
Get attributes computed by this functor.
Definition: SPHCalcHydroForceFunctor.h:529
sphLib::SPHCalcHydroForceFunctor::isRelevantForTuning
bool isRelevantForTuning() override
Specifies whether the functor should be considered for the auto-tuning process.
Definition: SPHCalcHydroForceFunctor.h:30
sphLib::SPHCalcHydroForceFunctor::allowsNonNewton3
bool allowsNonNewton3() override
Specifies whether the functor is capable of non-Newton3-like functors.
Definition: SPHCalcHydroForceFunctor.h:34
sphLib::SPHCalcHydroForceFunctor::SoAFunctorVerlet
void SoAFunctorVerlet(autopas::SoAView< SoAArraysType > soa, const size_t indexFirst, const std::vector< size_t, autopas::AlignedAllocator< size_t > > &neighborList, bool newton3) override
PairwiseFunctor for structure of arrays (SoA) for neighbor lists.
Definition: SPHCalcHydroForceFunctor.h:372
sphLib::SPHKernels::gradW
static std::array< double, 3 > gradW(const std::array< double, 3 > &dr, const double h)
gradient of the kernel function W
Definition: SPHKernels.h:75
sphLib::SPHKernels::getKernelSupportRadius
static double getKernelSupportRadius()
Get the kernelSupportRadius.
Definition: SPHKernels.h:28
autopas::utils::ArrayMath::dot
constexpr T dot(const std::array< T, SIZE > &a, const std::array< T, SIZE > &b)
Generates the dot product of two arrays.
Definition: ArrayMath.h:233
autopas::utils::ArrayMath::L2Norm
constexpr T L2Norm(const std::array< T, SIZE > &a)
Computes the L2Norm / Euclidean norm.
Definition: ArrayMath.h:265
autopas
This is the main namespace of AutoPas.
Definition: AutoPasDecl.h:32
autopas::OwnershipState::dummy
@ dummy
Dummy or deleted state, a particle with this state is not an actual particle!
Generated by doxygen 1.9.7