Peano 4
Loading...
Searching...
No Matches
Density.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_end_density(Particle* localParticle) {
7
8 if (localParticle->getSmoothingLengthConverged())
9 return;
10
11 // macro expansions from swift2/legacy/kernels are cumbersome otherwise.
13
14 /* Grab static parameters */
15 const int hydroDimensions = localParticle->getHydroDimensions();
16 const double etaFactor = localParticle->getEtaFactor();
17 const double alphaAV = localParticle->getAlphaAV();
18
19 /* Retrieve some properties */
20 const double m_i = localParticle->getMass();
21 const double h_i = localParticle->getSmoothingLength();
22
23 /* Some smoothing length multiples */
24 const double h_inv = 1.0 / h_i;
25 const double h_inv_dim = swift2::kernels::legacy::hydroDimensions::
26 pow_dimension(h_inv);
27 const double h_inv_dim_plus_one = h_inv_dim * h_inv;
28
29 /* Final operation. Add the self-contributions (i.e. at q=0) */
30 localParticle->setDensity(localParticle->getDensity() + m_i * kernel_root);
31 localParticle->setRho_dh(
32 localParticle->getRho_dh() - m_i * hydroDimensions * kernel_root
33 );
34 localParticle->setWcount(localParticle->getWcount() + kernel_root);
35 localParticle->setWcount_dh(
36 localParticle->getWcount_dh() - hydroDimensions * kernel_root
37 );
38
39
40 localParticle->setDensity(localParticle->getDensity() * h_inv_dim);
41 localParticle->setRho_dh(localParticle->getRho_dh() * h_inv_dim_plus_one);
42 localParticle->setWcount(localParticle->getWcount() * h_inv_dim);
43 localParticle->setWcount_dh(
44 localParticle->getWcount_dh() * h_inv_dim_plus_one
45 );
46
47 const double rho_inv = 1.0 / localParticle->getDensity();
48 // @TODO: implement cosmological factors
49 const double a_inv2 = 1.0;
50
51 /* Finish calculation of the velocity curl components */
52 localParticle->setRot_v(
53 localParticle->getRot_v() * h_inv_dim_plus_one * a_inv2 * rho_inv
54 );
55
56 /* Finish calculation of the velocity divergence */
57 localParticle->setDiv_v(
58 localParticle->getDiv_v() * h_inv_dim_plus_one * a_inv2 * rho_inv
59 );
60
62 std::isfinite(localParticle->getDensity()),
63 localParticle->toString()
64 );
65 assertion1(localParticle->getDensity() >= 0., localParticle->toString());
67 std::isfinite(localParticle->getWcount()),
68 localParticle->toString()
69 );
70 assertion1(localParticle->getWcount() >= 0., localParticle->toString());
72 std::isfinite(localParticle->getRho_dh()),
73 localParticle->toString()
74 );
76 std::isfinite(localParticle->getWcount_dh()),
77 localParticle->toString()
78 );
79}
80
81
82template <typename Particle>
84 // temporary to get things to work
85 if (localParticle->getSmoothingLengthConverged())
86 return;
87
88 localParticle->setWcount(0.0);
89 localParticle->setWcount_dh(0.0);
90
91 localParticle->setDensity(0.0);
92 localParticle->setRho_dh(0.0);
93
94 localParticle->setDiv_v(0.0);
95 localParticle->setRot_v(0.0);
96
97#if PeanoDebug > 0
98 localParticle->setDensityNeighbourCount(0);
99#endif
100}
101
102
103template <typename Particle>
105 Particle* localParticle,
106 const Particle* activeParticle
107) {
108
109 if (localParticle->getSmoothingLengthConverged())
110 return;
111
112 // macro expansions from swift2/legacy/kernels are cumbersome otherwise.
114
115 const double hydroDimensions = localParticle->getHydroDimensions();
116
117 /* Retrieve basic properties of local and active particles */
118 const double hi = localParticle->getSmoothingLength();
119
120 /* Distance between the particles */
122 dx = localParticle->getX() - activeParticle->getX();
123 const double r = tarch::la::norm2(dx);
124
125 /* Compute the actual interaction radius */
126 const double iactR = kernel_gamma * hi;
127
128 // This will probably cause issues with vectorisation and should be outside of
129 // the loop.
130 if (
132 and tarch::la::greater(r, 0.0) /* Self-contribution is calculated elsewhere
133 */
134 // and
135 // not localParticle->getIsBoundaryParticle()
136 ) {
137
138 const double mj = activeParticle->getMass();
139 const double hi_inv = 1. / hi;
140 // Normalized distance to evaluate kernel
141 const double qi = r * hi_inv;
142
143 // Evaluate kernel and its derivative
144 double wi, dwi_dx;
145 kernel_deval(qi, wi, dwi_dx);
146 const double dwi_dh = -(hydroDimensions * wi + qi * dwi_dx);
147
148 // Increment density
149 localParticle->setDensity(localParticle->getDensity() + mj * wi);
150
151 // Increment drho/dh. Note that minus sign was absorbed into dwi_dh
152 localParticle->setRho_dh(localParticle->getRho_dh() + mj * dwi_dh);
153
154 // Increment fractional number density
155 localParticle->setWcount(localParticle->getWcount() + wi);
156
157 // Increment dWcount/dh. Note that minus sign was absorbed into dwi_dh
158 localParticle->setWcount_dh(localParticle->getWcount_dh() + dwi_dh);
159
160
161 // Compute dv dot r
162 // Velocity terms for viscosity
164
165 /* Now we need to compute the div terms */
166 const double r_inv = r ? 1.0 / r : 0.0;
167 const double fac_i = mj * dwi_dx * r_inv;
168
169 dv = localParticle->getV() - activeParticle->getV();
170 const double dvdr = tarch::la::dot(dv, dx);
171 localParticle->setDiv_v(localParticle->getDiv_v() - fac_i * dvdr);
172
173 /* Compute dv cross r */
174 /* 2D Term only has 1 component (z) */
175#if Dimensions < 3
176 const double curlvr2D = dv[0] * dx[1] - dv[1] * dx[0];
177 localParticle->setRot_v(localParticle->getRot_v() + fac_i * curlvr2D);
178#else
179 const tarch::la::Vector<3, double> curlvr = {
180 dv[1] * dx[2] - dv[2] * dx[1],
181 dv[2] * dx[0] - dv[0] * dx[2],
182 dv[0] * dx[1] - dv[1] * dx[0]};
183 localParticle->setRot_v(localParticle->getRot_v() + fac_i * curlvr);
184#endif
185
186#if PeanoDebug > 0
187 localParticle->setDensityNeighbourCount(
188 localParticle->getDensityNeighbourCount() + 1
189 );
190#endif
191
192 } /* End if r<iactR */
193}
#define assertion1(expr, param)
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
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.
#define kernel_root
InlineMethod T pow_dimension(T x)
returns x^ndim
void hydro_end_density(Particle *localParticle)
Finishes the density calculation.
Definition Density.cpph:6
void hydro_prepare_density(Particle *localParticle)
Prepares a particle for the density calculation.
Definition Density.cpph:83
void density_kernel(Particle *localParticle, const Particle *activeParticle)
The actual density kernel, which interacts a local particle and an active particle,...
Definition Density.cpph:104
Scalar dot(const Vector< Size, Scalar > &lhs, const Vector< Size, Scalar > &rhs)
Performs the dot (=inner) product of two vectors.
bool greater(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
bool smallerEquals(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
Scalar norm2(const Vector< Size, Scalar > &vector)
Computes the 2-norm of the vector, i.e.
Simple vector class.
Definition Vector.h:134