PeriDyno 1.0.0
Loading...
Searching...
No Matches
DivergenceFreeSphSolver.cu
Go to the documentation of this file.
1#include "DivergenceFreeSphSolver.h"
2
3#include "SummationDensity.h"
4
5namespace dyno
6{
7 template<typename TDataType>
8 DivergenceFreeSphSolver<TDataType>::DivergenceFreeSphSolver()
9 : ParticleApproximation<TDataType>()
10 {
11 this->varRestDensity()->setValue(Real(1000));
12
13 mSummation = std::make_shared<SummationDensity<TDataType>>();
14
15 this->inSmoothingLength()->connect(mSummation->inSmoothingLength());
16 this->inSamplingDistance()->connect(mSummation->inSamplingDistance());
17 this->inPosition()->connect(mSummation->inPosition());
18 this->inNeighborIds()->connect(mSummation->inNeighborIds());
19
20 mSummation->outDensity()->connect(this->outDensity());
21 }
22
23 template<typename TDataType>
24 DivergenceFreeSphSolver<TDataType>::~DivergenceFreeSphSolver()
25 {
26 mKappa_r.clear();
27 mKappa_v.clear();
28 mAlpha.clear();
29 mPredictDensity.clear();
30 mDivergence.clear();
31 }
32
33
34
35 template <typename Real, typename Coord, typename Kernel>
36 __global__ void DFSPH_AlphaCompute(
37 DArray<Real> alpha,
38 DArray<Coord> posArr,
39 DArray<Real> rhoArr,
40 DArrayList<int> neighbors,
41 Real rest_density,
42 Real mass,
43 Real smoothingLength,
44 Kernel gradient,
45 Real scale
46 ) {
47 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
48 if (pId >= alpha.size()) return;
49
50 Coord pos_i = posArr[pId];
51
52 Real inv_alpha_i = Real(0);
53 Coord grad_ci(0);
54
55 List<int>& list_i = neighbors[pId];
56 int nbSize = list_i.size();
57
58 for (int ne = 0; ne < nbSize; ne++)
59 {
60 int j = list_i[ne];
61 Real r = (pos_i - posArr[j]).norm();
62
63 if (r > EPSILON)
64 {
65 Coord g = mass * gradient(r, smoothingLength, scale) * (pos_i - posArr[j]) * (1.0f / r);
66 grad_ci += g;
67 inv_alpha_i += g.dot(g);
68 }
69 }
70
71 inv_alpha_i += grad_ci.dot(grad_ci);
72
73 Real rho_i = rhoArr[pId];
74
75 if (inv_alpha_i < EPSILON) inv_alpha_i = EPSILON;
76
77 alpha[pId] = rho_i / (inv_alpha_i);
78
79 if ((nbSize < 5)) alpha[pId] = 0.0;
80
81 //if (alpha[pId] > 1.0)
82 //{
83 // printf("a %f, inv_a %e, rho %f, nb %d \r\n", alpha[pId], inv_alpha_i, rho_i, nbSize);
84 //}
85 }
86
87
88 template <typename Real, typename Coord, typename Kernel>
89 __global__ void DFSPH_DensityPredict(
90 DArray<Real> p_rhoArr,
91 DArray<Coord> posArr,
92 DArray<Coord> veloArr,
93 DArray<Real> rhoArr,
94 DArrayList<int> neighbors,
95 Real dt,
96 Real rest_density,
97 Real mass,
98 Real smoothingLength,
99 Kernel gradient,
100 Real scale
101 ) {
102
103 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
104 if (pId >= p_rhoArr.size()) return;
105
106 Coord pos_i = posArr[pId];
107
108 Real div_i = 0.0f;
109
110 List<int>& list_i = neighbors[pId];
111 int nbSize = list_i.size();
112
113 for (int ne = 0; ne < nbSize; ne++)
114 {
115 int j = list_i[ne];
116 Real r = (pos_i - posArr[j]).norm();
117
118 if (r > EPSILON)
119 {
120 Coord g = gradient(r, smoothingLength, scale) * (pos_i - posArr[j]) * (1.0f / r);
121 div_i += mass * g.dot(veloArr[pId] - veloArr[j]);
122 }
123 }
124 Real rho_i = rhoArr[pId] > rest_density ? rhoArr[pId] : rest_density;
125
126 //if (nbSize < 10) div_i = 0;
127
128 p_rhoArr[pId] = rho_i + dt * div_i;
129
130 if (p_rhoArr[pId] < rest_density) p_rhoArr[pId] = rest_density;
131
132 }
133
134 template <typename Real, typename Coord, typename Kernel>
135 __global__ void DFSPH_DivergencePredict(
136 DArray<Real> divergenceArr,
137 DArray<Coord> posArr,
138 DArray<Coord> veloArr,
139 DArray<Real> rhoArr,
140 DArrayList<int> neighbors,
141 Real dt,
142 Real rest_density,
143 Real mass,
144 Real smoothingLength,
145 Kernel gradient,
146 Real scale
147 ) {
148
149 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
150 if (pId >= divergenceArr.size()) return;
151
152 Coord pos_i = posArr[pId];
153
154 Real div_i = 0.0f;
155
156 List<int>& list_i = neighbors[pId];
157 int nbSize = list_i.size();
158
159 for (int ne = 0; ne < nbSize; ne++)
160 {
161 int j = list_i[ne];
162 Real r = (pos_i - posArr[j]).norm();
163
164 if (r > 100 * EPSILON)
165 {
166 Coord g = gradient(r, smoothingLength, scale) * (pos_i - posArr[j]) * (1.0f / r);
167 div_i += mass * g.dot(veloArr[pId] - veloArr[j]);
168 }
169 }
170
171 if ((div_i < 0) || (nbSize < 10)) div_i = 0;
172
173 divergenceArr[pId] = div_i;
174 }
175
176
177 template <typename Real>
178 __global__ void DFSPH_DensityKappa(
179 DArray<Real> KappaArr,
180 DArray<Real> alpahArr,
181 DArray<Real> predict_RhoArr,
182 DArray<Real> rhoArr,
183 Real dt,
184 Real rho_0
185 ) {
186 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
187 if (pId >= KappaArr.size()) return;
188 KappaArr[pId] = alpahArr[pId] * (predict_RhoArr[pId] - rho_0) / (dt * dt);
189 }
190
191
192
193 template <typename Real>
194 __global__ void DFSPH_DivergenceKappa(
195 DArray<Real> KappaArr,
196 DArray<Real> alpahArr,
197 DArray<Real> divergence,
198 DArray<Real> rhoArr,
199 Real dt,
200 Real rho_0
201 ) {
202 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
203 if (pId >= KappaArr.size()) return;
204
205 KappaArr[pId] = 1.0 * alpahArr[pId] * divergence[pId] / (dt);
206 }
207
208 template <typename Real, typename Coord, typename Kernel>
209 __global__ void DFSPH_VelocityUpdatedByKappa(
210 DArray<Coord> veloArr,
211 DArray<Real> KappaArr,
212 DArray<Real> alpahArr,
213 DArray<Coord> posArr,
214 DArray<Real> predict_RhoArr,
215 DArray<Real> rhoArr,
216 DArrayList<int> neighbors,
217 Real dt,
218 Real rest_density,
219 Real mass,
220 Real smoothingLength,
221 Kernel gradient,
222 Real scale
223 ) {
224 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
225 if (pId >= veloArr.size()) return;
226
227 Coord pos_i = posArr[pId];
228
229 List<int>& list_i = neighbors[pId];
230 int nbSize = list_i.size();
231
232 Coord d_velo(0.0f);
233
234 for (int ne = 0; ne < nbSize; ne++)
235 {
236 int j = list_i[ne];
237 Real r = (pos_i - posArr[j]).norm();
238 if (r > EPSILON)
239 {
240 Coord g = gradient(r, smoothingLength, scale) * (pos_i - posArr[j]) * (1.0f / r);
241 {
242 //d_velo += mass * (KappaArr[pId] / rhoArr[pId] + KappaArr[j] / rhoArr[j]) * g;
243 d_velo += mass * (KappaArr[pId] / (rhoArr[pId]) + KappaArr[j] / (rhoArr[j]) ) * g;
244 }
245 }
246 }
247 veloArr[pId] -= dt * d_velo;
248 }
249
250 template<typename TDataType>
251 void DivergenceFreeSphSolver<TDataType>::compute()
252 {
253 int num = this->inPosition()->size();
254
255 if (this->outDensity()->size() != this->inPosition()->size())
256 this->outDensity()->resize(this->inPosition()->size());
257
258 if (mKappa_r.size() != this->inPosition()->size())
259 mKappa_r.resize(this->inPosition()->size());
260
261 if (mKappa_v.size() != this->inPosition()->size())
262 mKappa_v.resize(this->inPosition()->size());
263
264 if (mAlpha.size() != this->inPosition()->size())
265 mAlpha.resize(this->inPosition()->size());
266
267 if (mPredictDensity.size() != this->inPosition()->size())
268 mPredictDensity.resize(this->inPosition()->size());
269
270 if (mDivergence.size() != this->inPosition()->size())
271 mDivergence.resize(this->inPosition()->size());
272
273 mSummation->varRestDensity()->setValue(this->varRestDensity()->getValue());
274 mSummation->varKernelType()->setCurrentKey(this->varKernelType()->currentKey());
275 mSummation->update();
276
277 int MaxItNum = this->varMaxIterationNumber()->getValue();
278
279 this->computeAlpha();
280
281 int it = 0;
282 if (this->varDivergenceSolverDisabled()->getValue() == false)
283 {
284 Real v_err = 10000.0f;
285 while ((v_err > this->varDivergenceErrorThreshold()->getValue()) && (it < MaxItNum))
286 {
287 v_err = abs(this->takeOneDivergenIteration()) / this->varRestDensity()->getValue();
288 it++;
289 }
290 }
291
292 it = 0;
293 if (this->varDensitySolverDisabled()->getValue() == false)
294 {
295 Real d_err = 10000.0f;
296 while ((it < MaxItNum))
297 {
298 d_err = abs(this->takeOneDensityIteration());
299 it++;
300 }
301 }
302 }
303
304 template<typename TDataType>
305 void DivergenceFreeSphSolver<TDataType>::computeAlpha()
306 {
307
308 int num = this->inPosition()->size();
309 Real rho_0 = this->varRestDensity()->getValue();
310
311 cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
312 DFSPH_AlphaCompute,
313 mAlpha,
314 this->inPosition()->getData(),
315 mSummation->outDensity()->getData(),
316 this->inNeighborIds()->getData(),
317 rho_0,
318 mSummation->getParticleMass(),
319 this->inSmoothingLength()->getValue()
320 );
321 }
322
323 template<typename TDataType>
324 typename TDataType::Real DivergenceFreeSphSolver<TDataType>::takeOneDensityIteration()
325 {
326 Real dt = this->inTimeStep()->getData();
327 int num = this->inPosition()->size();
328 Real rho_0 = this->varRestDensity()->getValue();
329
330 mSummation->update();
331
332 cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
333 DFSPH_DensityPredict,
334 mPredictDensity,
335 this->inPosition()->getData(),
336 this->inVelocity()->getData(),
337 mSummation->outDensity()->getData(),
338 this->inNeighborIds()->getData(),
339 this->inTimeStep()->getValue(),
340 rho_0,
341 mSummation->getParticleMass(),
342 this->inSmoothingLength()->getValue()
343 );
344
345 cuExecute(num, DFSPH_DensityKappa,
346 mKappa_r,
347 mAlpha,
348 mPredictDensity,
349 mSummation->outDensity()->getData(),
350 this->inTimeStep()->getValue(),
351 rho_0
352 );
353
354 cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
355 DFSPH_VelocityUpdatedByKappa,
356 this->inVelocity()->getData(),
357 mKappa_r,
358 mAlpha,
359 this->inPosition()->getData(),
360 mPredictDensity,
361 mSummation->outDensity()->getData(),
362 this->inNeighborIds()->getData(),
363 this->inTimeStep()->getValue(),
364 rho_0,
365 mSummation->getParticleMass(),
366 this->inSmoothingLength()->getValue()
367 );
368
369// auto m_reduce = Reduction<Real>::Create(num);
370// Real avr_pdensity = m_reduce->average(mPredictDensity.begin(), num);
371// Real err = (avr_pdensity - rho_0) / rho_0;
372// delete m_reduce;
373 return 0;
374 }
375
376 template<typename TDataType>
377 typename TDataType::Real DivergenceFreeSphSolver<TDataType>::takeOneDivergenIteration()
378 {
379 Real err = 0.0f;
380
381 Real dt = this->inTimeStep()->getData();
382 int num = this->inPosition()->size();
383 Real rho_0 = this->varRestDensity()->getValue();
384
385 cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
386 DFSPH_DivergencePredict,
387 mDivergence,
388 this->inPosition()->getData(),
389 this->inVelocity()->getData(),
390 mSummation->outDensity()->getData(),
391 this->inNeighborIds()->getData(),
392 this->inTimeStep()->getValue(),
393 rho_0,
394 mSummation->getParticleMass(),
395 this->inSmoothingLength()->getValue()
396 );
397
398 cuExecute(num, DFSPH_DivergenceKappa,
399 mKappa_v,
400 mAlpha,
401 mDivergence,
402 mSummation->outDensity()->getData(),
403 this->inTimeStep()->getValue(),
404 rho_0
405 );
406
407 cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
408 DFSPH_VelocityUpdatedByKappa,
409 this->inVelocity()->getData(),
410 mKappa_v,
411 mAlpha,
412 this->inPosition()->getData(),
413 mPredictDensity,
414 mSummation->outDensity()->getData(),
415 this->inNeighborIds()->getData(),
416 this->inTimeStep()->getValue(),
417 rho_0,
418 mSummation->getParticleMass(),
419 this->inSmoothingLength()->getValue()
420 );
421
422 auto m_reduce2 = Reduction<Real>::Create(num);
423 Real avr_div = m_reduce2->average(mDivergence.begin(), num);
424 err = avr_div;
425 delete m_reduce2;
426 return err;
427 }
428
429 DEFINE_CLASS(DivergenceFreeSphSolver);
430}
Generated by doxygen 1.13.2