Peano 4
Loading...
Searching...
No Matches
HydroForce.cpph
Go to the documentation of this file.
1#include "equation_of_state.h"
2#include "hydro_dimensions.h"
3#include "kernel_hydro.h"
4
5template <typename Particle>
6void swift2::kernels::legacy::hydro_prepare_force(Particle* localParticle) {
7
8 static tarch::logging::Log _log("swift2::kernels::legacy");
9
10 // Grab static parameters
11 const int hydroDimensions = localParticle->getHydroDimensions();
12 const double alphaAV = localParticle->getAlphaAV();
13 double hydro_h_max = std::min(
14 localParticle->getSmlMax(),
15 localParticle->getSearchRadius() / kernel_gamma
16 );
17
18 // Retrieve basic properties of local and active particles
19 const double rho_i = localParticle->getDensity();
20 const double m_i = localParticle->getMass();
21 const double h_i = localParticle->getSmoothingLength();
22 const double u_i = localParticle->getU();
23
24 const double h_inv = 1.0 / h_i;
25
26 // This is needed for cosmological integration.
27 const double fac_balsara_eps = 1.;
28
29 // Compute the norm of the curl
30#if Dimensions < 3
31 const double curl_v = std::sqrt(
32 localParticle->getRot_v() * localParticle->getRot_v()
33 );
34#else
35 const double curl_v = tarch::la::norm2(localParticle->getRot_v());
36#endif
37
38 /* Compute the norm of div v */
39 const double abs_div_v = tarch::la::abs(localParticle->getDiv_v());
40
41 /* Compute the pressure */
42 const double P_i = eos::gas_pressure_from_internal_energy(rho_i, u_i);
43 assertion1(std::isfinite(P_i), localParticle->toString());
44 assertion1(P_i > 0., localParticle->toString());
45
46 /* Compute the sound speed */
47 const double soundspeed = legacy::eos::gas_soundspeed_from_pressure(
48 rho_i,
49 P_i
50 );
51 assertion1(std::isfinite(soundspeed), localParticle->toString());
52
53 /* Compute the "grad h" term - Note here that we have \tilde{x}
54 * as 1 as we use the local number density to find neighbours. This
55 * introduces a j-component that is considered in the force loop,
56 * meaning that this cached grad_h_term gives:
57 *
58 * f_ij = 1.f - grad_h_term_i / m_j */
59 const double common_factor = h_i / hydroDimensions
60 / localParticle->getWcount();
61
62 float grad_h_term;
63 // Ignore changing-kernel effects when h ~= h_max
64 if (h_i > 0.9999 * hydro_h_max) {
65 grad_h_term = 0.;
67 "hydro_prepare_force",
68 "h ~ h_max for particle with ID " << localParticle->getPartid(
69 ) << " h = " << h_i
70 );
71 } else {
72 const double grad_W_term = common_factor * localParticle->getWcount_dh();
73 if (grad_W_term < -0.9999) {
74 /* if we get here, we either had very small neighbour contributions
75 (which should be treated as a no neighbour case in the ghost) or
76 a very weird particle distribution (e.g. particles sitting on
77 top of each other). Either way, we cannot use the normal
78 expression, since that would lead to overflow or excessive round
79 off and cause excessively high accelerations in the force loop */
80 grad_h_term = 0.;
82 "hydro_prepare_force",
83 "grad_W_term very small for particle with ID "
84 << localParticle->getPartid() << " h:" << h_i
85 << " wcount:" << localParticle->getWcount() << " wcount_dh:"
86 << localParticle->getWcount_dh() << " grad_W_term:" << grad_W_term
87 << " common_factor:" << common_factor
88 );
89 } else {
90 grad_h_term = common_factor * localParticle->getRho_dh()
91 / (1. + grad_W_term);
92 }
93 }
94
95
96 /* Compute the Balsara Switch */
97 const double balsara
98 = alphaAV * abs_div_v
99 / (abs_div_v + curl_v + 0.0001 * fac_balsara_eps * soundspeed * h_inv);
100
102 std::isfinite(balsara),
103 localParticle->toString(),
104 balsara,
105 abs_div_v,
106 curl_v,
107 soundspeed,
108 alphaAV
109 );
110
111 /* Update variables */
112 localParticle->setF(grad_h_term);
113 localParticle->setPressure(P_i);
114 localParticle->setSoundSpeed(soundspeed);
115 localParticle->setBalsara(balsara);
116
117#if PeanoDebug > 0
118 localParticle->setForceNeighbourCount(0);
119#endif
120}
121
122
123template <typename Particle>
125) {
126
127 /* Reset the acceleration*/
128 localParticle->setA(0.0);
129
130 /* Reset the time derivatives */
131 localParticle->setUDot(0.0);
132 localParticle->setHDot(0.0);
133 localParticle->setV_sig_AV(2.0 * localParticle->getSoundSpeed());
134}
135
136
137template <typename Particle>
138void swift2::kernels::legacy::hydro_end_force(Particle* localParticle) {
139
140 const int hydroDimensions = localParticle->getHydroDimensions();
141 const double hydroDimensionsInv = 1. / (double)hydroDimensions;
142
143 localParticle->setHDot(
144 localParticle->getHDot() * localParticle->getSmoothingLength()
145 * hydroDimensionsInv
146 );
148 std::isfinite(localParticle->getHDot()),
149 localParticle->toString(),
150 localParticle->getSmoothingLength()
151 );
152
153#if PeanoDebug > 0
155 localParticle->getDensityNeighbourCount(
156 ) <= localParticle->getForceNeighbourCount(),
157 "Interacted with fewer particles in Force than in Density loop",
158 localParticle->toString()
159 );
160#endif
161}
162
163
164template <typename Particle>
166 Particle* localParticle,
167 const Particle* activeParticle
168) {
169
170 // macro expansions from swift2/legacy/kernels are cumbersome otherwise.
171 using namespace kernelHydro;
172
173 // Grab static parameters
174 const int hydroDimensions = localParticle->getHydroDimensions();
175 const double alphaAV = localParticle->getAlphaAV();
176 const double betaAV = localParticle->getBetaAV();
177
178 // No cosmology yet: a^{(3 gamma - 5 )/ 2} = 1.
179 const double fac_mu = 1.;
180 const double a2_Hubble = 0.; // a**2 * H, where H = da/dt / a. = 0 without
181 // cosmo
182
183 // Distance between the two particles
185 dx = localParticle->getX() - activeParticle->getX();
186
187 const double r = tarch::la::norm2(dx);
188 const double r_inv = r ? 1.0 / r : 0.0;
189
190 // Get actual interaction radius
191 const double h_i = localParticle->getSmoothingLength();
192 const double iactRi = kernel_gamma * h_i;
193 const double h_j = activeParticle->getSmoothingLength();
194 const double iactRj = kernel_gamma * h_j;
195
196 if ((tarch::la::smaller(r, iactRi) or (tarch::la::smaller(r, iactRj))) and tarch::la::greater(r, 0.0)) {
197
198#if PeanoDebug > 0
199 localParticle->setForceNeighbourCount(
200 localParticle->getForceNeighbourCount() + 1
201 );
202#endif
203
204 // Retrieve basic properties of local and active particles
205 const double m_i = localParticle->getMass();
206 const double m_j = activeParticle->getMass();
207
208 const double rho_i = localParticle->getDensity();
209 const double rho_j = activeParticle->getDensity();
210
211 const double P_i = localParticle->getPressure();
212 const double P_j = activeParticle->getPressure();
213
214 // Get kernel for hi
215 const double h_i_inv = 1. / h_i;
216 const double h_i_inv_dp1 = swift2::kernels::legacy::hydroDimensions::
217 pow_dimension_plus_one(h_i_inv); // 1/h^{dim + 1}
218 const double qi = r * h_i_inv;
219 double wi, dwi_dqi;
220 // Evaluate kernel and its derivative for particle a
221 kernel_deval(qi, wi, dwi_dqi);
222 const double dwi_dr = dwi_dqi * h_i_inv_dp1;
223
224 // Get kernel for hj
225 const double h_j_inv = 1. / h_j;
226 const double h_j_inv_dp1 = swift2::kernels::legacy::hydroDimensions::
227 pow_dimension_plus_one(h_j_inv); // 1/h^{dim + 1}
228 const double qj = r * h_j_inv;
229 double wj, dwj_dqj;
230 // Evaluate kernel and its derivative for particle b
231 kernel_deval(qj, wj, dwj_dqj);
232 const double dwj_dr = dwj_dqj * h_j_inv_dp1;
233
234 // Variable smoothing length term
235
236 // SWIFT "f"
237 const double f_ij = 1.0 - localParticle->getF() / m_j;
238 const double f_ji = 1.0 - activeParticle->getF() / m_i;
239
240 // Compute gradient terms
241 const double P_over_rho2_i = P_i / (rho_i * rho_i) * f_ij;
242 const double P_over_rho2_j = P_j / (rho_j * rho_j) * f_ji;
243
244 // Compute dv dot r.
245 double dvdr = tarch::la::dot(
246 (localParticle->getV() - activeParticle->getV()),
247 dx
248 );
249
250 // Add Hubble flow: Skipped for now
251 const double dvdr_Hubble = dvdr + a2_Hubble * r * r;
252
253 /* Are the particles moving towards each others ? */
254 const double omega_ij = std::min(dvdr_Hubble, 0.);
255 const double mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
256
257 /* Compute sound speeds and signal velocity */
258 const double v_sig = localParticle->getSoundSpeed()
259 + activeParticle->getSoundSpeed() - betaAV * mu_ij;
260
261 /* Balsara term */
262 const double balsara_i = localParticle->getBalsara();
263 const double balsara_j = activeParticle->getBalsara();
264
265 /* Construct the full viscosity term */
266 const double rho_ij = 0.5 * (rho_i + rho_j);
267 const double visc = -0.25 * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
268
269 /* Convolve with the kernel */
270 const double visc_acc_term = 0.5 * visc * (dwi_dr * f_ij + dwj_dr * f_ji)
271 * r_inv;
272
273 /* SPH acceleration term */
274 const double sph_acc_term = (P_over_rho2_i * dwi_dr + P_over_rho2_j * dwj_dr)
275 * r_inv;
276
277 // Assemble the acceleration
278 const double acc = sph_acc_term + visc_acc_term;
279
281 std::isfinite(visc_acc_term),
282 visc_acc_term,
283 localParticle->toString(),
284 activeParticle->toString(),
285 dwi_dr,
286 f_ij,
287 dwj_dr,
288 f_ji,
289 r_inv,
290 visc
291 );
293 std::isfinite(sph_acc_term),
294 sph_acc_term,
295 localParticle->toString(),
296 activeParticle->toString()
297 );
298
299 // Use the force
300 // NOTE: This is a_hydro in SWIFT
301 localParticle->setA(localParticle->getA() - m_j * acc * dx);
302
303 /* Get the time derivative for u. */
304 const double sph_du_term_i = P_over_rho2_i * dvdr * r_inv * dwi_dr;
305
306 /* Viscosity term */
307 const double visc_du_term = 0.5 * visc_acc_term * dvdr_Hubble;
308
309 // Assemble the energy equation term
310 const double du_dt_i = sph_du_term_i + visc_du_term;
311
312 // Internal energy time derivative
313 localParticle->setUDot(localParticle->getUDot() + m_j * du_dt_i);
314
315 /* Get the time derivative for h. */
316 localParticle->setHDot(
317 localParticle->getHDot() - m_j * dvdr * r_inv / rho_j * dwi_dr * f_ij
318 );
319
320 // Update the signal velocity
321 localParticle->setV_sig_AV(std::max(localParticle->getV_sig_AV(), v_sig));
322
324 not std::isnan(localParticle->getHDot()),
325 localParticle->toString(),
326 dvdr,
327 r_inv,
328 dwi_dr,
329 f_ij,
330 f_ji
331 );
333 not tarch::la::isEntryNan(localParticle->getA()),
334 localParticle->toString(),
335 sph_acc_term,
336 visc_acc_term,
337 f_ij,
338 f_ji,
339 activeParticle->getF(),
340 activeParticle->getIsBoundaryParticle()
341 );
342 }
343}
344
345
346template <typename Particle>
349 Particle& localParticle,
350 const Particle& activeParticle
351) {
352
353 // macro expansions from swift2/legacy/kernels are cumbersome otherwise.
355
356 // Grab static parameters
357 const int hydroDimensions = localParticle.getHydroDimensions();
358 const double alphaAV = localParticle.getAlphaAV();
359 const double betaAV = localParticle.getBetaAV();
360
361 // Retrieve basic properties of local and active particles
362 const double rho_i = localParticle.getDensity();
363 const double rho_j = activeParticle.getDensity();
364
365 const double P_i = localParticle.getPressure();
366 const double P_j = activeParticle.getPressure();
367
368 const double m_i = localParticle.getMass();
369 const double m_j = activeParticle.getMass();
370
371 const double h_i = localParticle.getSmoothingLength();
372 const double h_j = activeParticle.getSmoothingLength();
373
374 // Distance between the two particles
376 dx = localParticle.getX() - activeParticle.getX();
377 const double r = tarch::la::norm2(dx);
378 const double r_inv = r ? 1.0 / r : 0.0;
379
380 // Dimensionless distance to evaluate kernel
381 const double h_i_inv = 1. / h_i;
382 const double h_j_inv = 1. / h_j;
383 const double q_i = r * h_i_inv;
384 const double q_j = r * h_j_inv;
385
386 // Get actual interaction radius
387 const double iactRi = kernel_gamma * h_i;
388 const double iactRj = kernel_gamma * h_j;
389
390 // @Pawel I don't know if this vectorises anymore, but I think the second
391 // check is important
392 // I also don't know if a strict or non-strict boolean optimisation
393 // make a difference. If they do, please document it in the SWIFT
394 // optimisation pages.
395 double mask = ((tarch::la::smaller(r, iactRi) | tarch::la::smaller(r, iactRj))
396 & tarch::la::greater(r, 0.0))
397 ? 1.0
398 : 0.0;
399
400 // Get actual interaction radius
401 double wi, wj;
402 double wi_dr, wj_dr;
403 double dw_dqi, dw_dqj;
404
405 // Evaluate kernel and its derivative for particle a
406 kernel_deval(q_i, wi, dw_dqi);
407
408 // Evaluate kernel and its derivative for particle b
409 kernel_deval(q_j, wj, dw_dqj);
410
411 wi_dr = dw_dqi
414 wj_dr = dw_dqj
417
418 // Variable smoothing length term
419
420 // SWIFT "f"
421 const double f_ij = 1.0 - localParticle.getF() / m_j;
422 const double f_ji = 1.0 - activeParticle.getF() / m_i;
423
424 // Compute gradient terms
425 const double P_over_rho2_i = P_i / (rho_i * rho_i) * f_ij;
426 const double P_over_rho2_j = P_j / (rho_j * rho_j) * f_ji;
427
428 // Compute dv dot r.
429 double dvdr = tarch::la::dot(
430 (localParticle.getV() - activeParticle.getV()),
431 dx
432 );
433
434
435 // Are the particles moving towards each others ?
436 const double omega_ij = std::min(dvdr, 0.0);
437 const double mu_ij = r_inv * omega_ij;
438
439 // Compute sound speeds and signal velocity
440 const double v_sig = localParticle.getSoundSpeed()
441 + activeParticle.getSoundSpeed() - betaAV * mu_ij;
442
443 // Balsara term
444 const double balsara_i = localParticle.getBalsara();
445 const double balsara_j = activeParticle.getBalsara();
446
447 // Construct the full viscosity term
448 const double rho_ij = 0.5 * (rho_i + rho_j);
449 const double visc = -0.25 * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
450
451 // Convolve with the kernel
452 const double visc_acc_term = 0.5 * visc * (wi_dr * f_ij + wj_dr * f_ji)
453 * r_inv;
454
455 // SPH acceleration term
456 const double sph_acc_term = (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr)
457 * r_inv;
458
459 // Assemble the acceleration
460 const double acc = sph_acc_term + visc_acc_term;
461
463 std::isfinite(visc_acc_term),
464 visc_acc_term,
465 localParticle.toString(),
466 activeParticle.toString(),
467 wi_dr,
468 f_ij,
469 wj_dr,
470 f_ji,
471 r_inv,
472 visc
473 );
475 std::isfinite(sph_acc_term),
476 sph_acc_term,
477 localParticle.toString(),
478 activeParticle.toString()
479 );
480
481 // Use the force
482 localParticle.setA(localParticle.getA() - mask * m_j * acc * dx);
483
484 // Get the time derivative for u.
485 const double sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
486
487 // Viscosity term
488 const double visc_du_term = 0.5 * visc_acc_term * dvdr;
489
490 // Assemble the energy equation term
491 const double du_dt_i = sph_du_term_i + visc_du_term;
492
493 const double udot = localParticle.getUDot();
494
495 // Internal energy time derivative
496 localParticle.setUDot(localParticle.getUDot() + mask * m_j * du_dt_i);
497
498 // Get the time derivative for h.
499 localParticle.setHDot(
500 localParticle.getHDot() - mask * m_j * dvdr * r_inv / rho_j * wi_dr * f_ij
501 );
502
503 // Update the signal velocity
504 // localParticle.setV_sig_AV(std::max(localParticle.getV_sig_AV(), mask *
505 // v_sig));
506 const double currentAV = localParticle.getV_sig_AV();
507 const double biggerAV = mask * v_sig;
508 const double newAV = biggerAV > currentAV ? biggerAV : currentAV;
509 localParticle.setV_sig_AV(newAV);
510
512 not std::isnan(localParticle.getHDot()),
513 localParticle.toString(),
514 dvdr,
515 r_inv,
516 wi_dr,
517 f_ij,
518 f_ji
519 );
521 not tarch::la::isEntryNan(localParticle.getA()),
522 localParticle.toString(),
523 sph_acc_term,
524 visc_acc_term,
525 f_ij,
526 f_ji,
527 activeParticle.getF(),
528 activeParticle.getIsBoundaryParticle()
529 );
530}
#define assertion2(expr, param0, param1)
#define assertion3(expr, param0, param1, param2)
#define assertion7(expr, param0, param1, param2, param3, param4, param5, param6)
#define assertion1(expr, param)
#define assertion9(expr, param0, param1, param2, param3, param4, param5, param6, param7, param8)
#define assertion6(expr, param0, param1, param2, param3, param4, param5)
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ so we are integrating with f$ phi_k phi_k dx
#define logWarning(methodName, logMacroMessageStream)
Wrapper macro around tarch::tarch::logging::Log to improve logging.
Definition Log.h:440
Log Device.
Definition Log.h:516
contains functions and definitions related to the equation of state of ideal gases.
File containing macro-riddled functions related to dimensions.
Kernel functions for SPH.
InlineMethod T pow_dimension_plus_one(T x)
returns x^{ndim+1}
void hydro_reset_acceleration(Particle *localParticle)
Reset acceleration fields of a particle.
void force_kernel(Particle *localParticle, const Particle *activeParticle)
The actual force kernel, which interacts a local particle and an active particle, and updates the loc...
void hydro_end_force(Particle *localParticle)
Finishes the force calculation.
void forceKernelWithMasking(const peano4::datamanagement::CellMarker &marker, Particle &localParticle, const Particle &activeParticle)
Vectorised alternative implementation of forceKernel()
void hydro_prepare_force(Particle *localParticle)
Prepare a particle for the force calculation.
Scalar dot(const Vector< Size, Scalar > &lhs, const Vector< Size, Scalar > &rhs)
Performs the dot (=inner) product of two vectors.
static bool smaller(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE) InlineMethod
Smaller operator for floating point values.
bool greater(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
double abs(double value)
Returns the absolute value of a type by redirecting to std::abs.
int isEntryNan(const Vector< Size, Scalar > &vector)
Scalar norm2(const Vector< Size, Scalar > &vector)
Computes the 2-norm of the vector, i.e.
tarch::logging::Log _log("examples::unittests")
Simple vector class.
Definition Vector.h:134